Development of AMBER Parameters for Molecular Simulations of Selected Boron-Based Covalent Ligands

Boron containing compounds (BCCs) aroused increasing interest in the scientific community due to their wide application as drugs in various fields. In order to design new compounds hopefully endowed with pharmacological activity and also investigate their conformational behavior, the support of computational studies is crucial. Nevertheless, the suitable molecular mechanics parameterization and the force fields needed to perform these simulations are not completely available for this class of molecules. In this paper, Amber force field parameters for phenyl-, benzyl-, benzylamino-, and methylamino-boronates, a group of boron-containing compounds involved in different branches of the medicinal chemistry, were created. The robustness of the obtained data was confirmed through molecular dynamics simulations on ligand/β-lactamases covalent complexes. The ligand torsional angles, populated over the trajectory frames, were confirmed by values found in the ligand geometries, located through optimizations at the DFT/B3LYP/6-31g(d) level, using water as a solvent. In summary, this study successfully provided a library of parameters, opening the possibility to perform molecular dynamics simulations of this class of boron-containing compounds.


Introduction
Boron is an element widely distributed in nature. It plays essential functions, being involved in the growth, development, and metabolism of plants. Moreover, it is able to regulate vitamin D levels, help brain function, reduce the risk of developing cancer, and promote bone health in mammals [1]. Concerning this last aspect, it was demonstrated that diets deficient in boron hinder bone formation with respect to control diets that were supplemented with boron. Nevertheless, the exact mechanism behind this effect remains unclear [2].
The first isolated natural product, containing trace amounts of boron, was Boromycin. This is a macrolide that was isolated in the African soil from a strain of Streptomyces antibioticus. This antibiotic has been extensively studied for its therapeutic properties. In fact, it showed nanomolar potency against several HIV-infected cell lines [3] and against various bacterial strains, including Mycobacterium tuberculosis [4]. A compound, structurally related to Boromycin, is Aplasmomycin. It was isolated from Streptomyces griseus and has anti-plasmodium activity. In both the abovementioned natural products, the boron atom plays a structural role, causing the polyols to fold into compact structures. All these considerations suggest the increasing role that boron-containing compounds (BCCs) have assumed in scientific research.
Indeed, in the last years, BCCs have attracted growing attention because their clinical applications span from the treatment of fungal infection (i.e., tavaborole, Kerydin ® , Pfizer, New York, NY, USA) [5] to the treatment of cancer (i.e., velcade, Bortezomib ® , Janssen-Cilag International N.V., Olen, Belgium) [6]. Moreover, as evidence of their significant pharmacological activity, they also act as β-lactamase inhibitors (i.e., vaborbactam, in association with meropenem, Vaborem ® , Menarini International Operation Luxembourg S.A., Luxembourg) [7], preventing the antibiotic-cleavage activity of antibiotic-resistant bacterial strains [8][9][10][11]. Additionally, they are able to form covalent complexes with other serine-proteases, such as chymotrypsin, trypsin, and thrombin [12]. Recently, some researchers have reported on the activity of BCCs in the inhibition of SARS-CoV-2 M pro , paving the way for the development of new therapeutic options to fight viral infections [13].
The warhead of BCCs can also react with organic compounds containing the hydroxide group, giving rise to a wide range of molecules endowed with biological activity [14]. In particular, they can irreversibly bound to structures containing cis-hydroxyl groups, such as riboflavin, pyridoxine, adenosine monophosphate, pyrimidine nucleotides, ascorbic acid, ribose, and polysaccharides [15].
One of the primary factors that drives the increasing use of BCCs in research and development was also their ability to switch between neutral trigonal planar sp2 and tetrahedral sp3 hybridization states. This property enables them to adopt various binding modes during the target recognition process and makes them attractive as high-affinity ligands with low molecular mass. Between BCCs, boronic acids can form different types of covalent adducts with nucleophiles in target proteins, including trigonal covalent, tetragonal covalent, or bidentate covalent adducts. These products are reversible, so unplanned covalent modifications of non-target proteins are minimized. The two hydroxyl groups, present in the chemical structure of boronic acids, offer six opportunities to form hydrogen bond contacts with amino acid residues, increasing the ligand's affinity with the target protein. The ability of boron to alter its hybridization states enables boronate-based inhibitors to imitate the sp 2 state of β-lactamase substrates, allowing them to efficiently bind to them. They can subsequently react with the nucleophilic serine of serine-β-lactamases or metal-β-lactamases to form the sp 3 state, mimicking a high-energy intermediate. In the last case, the boron moiety can interact with metal ions in enzymes, which is useful in designing candidate drugs that can withstand selection pressures for drug resistance.
Despite the widely demonstrated pharmacological importance of BCCs, the appropriate molecular mechanics parameterization and the force fields, which allow us to perform computational studies essential for designing new ones and evaluate their conformational behavior, are not completely available to the researchers. Therefore, molecular dynamics (MD) simulations of BCCs, covalently bound to targets, can be performed only using computationally demanding QM-MM methods [16] or, as an alternative, by using the OPLS4 force fields (property of Schrödinger, LLC, New York, NY, USA) [17] available after the release of the commercial license.
In recent years, the parametrization of a limited number of aryl-, alkyl-boronic acids [12,18], and boronate esters [19] appeared in the literature. Tafi et al. [12] reported that the bonded, non-bonded, and point charges MacroModel/Amber force field was retrieved, in addition to GB/SA solvation parameters, for modeling boronic acids as tetrahedral adducts that are formed after coordination of the protease's serine Oγ. The new force field was validated through flexible docking studies conducted on three crystallographic complexes of β-lactamases with boronic acids. The output of these studies matched up with the crystallographic conformation of the complexes as the global minimum energy structure. In another paper, Kurt et al. [18] generated the Amber force field parameters for benzodioxaboroles, a group of boron compounds with aromatic structures. Their study produced the necessary parameter library for performing molecular dynamics simulations of this class of compounds. In fact, the root mean square deviation (RMSD) value between the minimized geometries and the x-ray structures was found to closely match the values obtained by quantum mechanical calculations. Using an anti-cancer BCC as a ligand, the molecular dynamics (MD) simulation of the DNA-ligand complex was then 3 of 14 successfully conducted, and the experimentally reported anti-cancer effect was confirmed by the simulations.
Despite these significant forward steps, to the best of our knowledge, the AMBER force field of phenyl-, benzyl-, benzylamino-, and methylamino-boronates are not available, even if these compounds are significant in numerous therapeutic fields. In fact, phenyl-and methyl-amino boronates are β-lactamases inhibitors [20] and anti-HIV agents [21]. Phenylboronates are also amazing prodrugs, like ZB483, able to produce a 40-fold increase in the endoxifen concentration peak in plasma when used in breast cancer therapy [22]. Based on their importance in medicinal chemistry, our attention has been focused on compounds 1-4, reported in Figure 1, to potentially open the way to the search of new drugs with the fundamental support of molecular modeling. simulations of this class of compounds. In fact, the root mean square deviation (RMSD) value between the minimized geometries and the x-ray structures was found to closely match the values obtained by quantum mechanical calculations. Using an anti-cancer BCC as a ligand, the molecular dynamics (MD) simulation of the DNA-ligand complex was then successfully conducted, and the experimentally reported anti-cancer effect was confirmed by the simulations. Despite these significant forward steps, to the best of our knowledge, the AMBER force field of phenyl-, benzyl-, benzylamino-, and methylamino-boronates are not available, even if these compounds are significant in numerous therapeutic fields. In fact, phenyl-and methyl-amino boronates are -lactamases inhibitors [20] and anti-HIV agents [21]. Phenyl-boronates are also amazing prodrugs, like ZB483, able to produce a 40-fold increase in the endoxifen concentration peak in plasma when used in breast cancer therapy [22]. Based on their importance in medicinal chemistry, our attention has been focused on compounds 1-4, reported in Figure 1, to potentially open the way to the search of new drugs with the fundamental support of molecular modeling. With this aim and to cover the above reported computational gap, in this paper we calculated the molecular mechanics parameters of the representative BCCs 1-4, using the Paramfit procedure, as suggested by AMBER developers [23]. The warheads of 1-4 were covalently bound to the side chain of a serine residue (in blue in Figure 1) to simulate the covalent bond (in red in Figure 1) to a putative serine-proteases. The accuracy of the retrieved parameters was established by performing MD simulations of the four BCCs in complex with AmpC -lactamases, produced by E. coli [24]. Finally, to verify the new force field performance, the BCC torsional angles values, populated during the MD simulations, were compared with those found in the geometries located by a DFT conformational analysis [25].

Results and Discussion
The Paramfit procedure [23] was applied to compounds 1-4 to generate the dihedral parameters shown in Table 1, while the bond and angle parameters, reported in Table 2, were retrieved from the literature data [12,19]. In particular, the Paramfit procedure was applied to predict the parameters useful to simulate the covalent bond connecting the serine residue and the BCCs. In fact, in a standard AMBER simulation, in which a non- With this aim and to cover the above reported computational gap, in this paper we calculated the molecular mechanics parameters of the representative BCCs 1-4, using the Paramfit procedure, as suggested by AMBER developers [23]. The warheads of 1-4 were covalently bound to the side chain of a serine residue (in blue in Figure 1) to simulate the covalent bond (in red in Figure 1) to a putative serine-proteases. The accuracy of the retrieved parameters was established by performing MD simulations of the four BCCs in complex with AmpC β-lactamases, produced by E. coli [24]. Finally, to verify the new force field performance, the BCC torsional angles values, populated during the MD simulations, were compared with those found in the geometries located by a DFT conformational analysis [25].

Results and Discussion
The Paramfit procedure [23] was applied to compounds 1-4 to generate the dihedral parameters shown in Table 1, while the bond and angle parameters, reported in Table 2, were retrieved from the literature data [12,19]. In particular, the Paramfit procedure was applied to predict the parameters useful to simulate the covalent bond connecting the serine residue and the BCCs. In fact, in a standard AMBER simulation, in which a non-covalent ligand was simulated in complex with a serine-protease, the ligand is parameterized by GAFF [26] or the new version GAFF2 [27], while the biological counterpart is simulated by applying one of the force fields available in AMBER for simulating proteins, such as the ff14SB [28]. Conversely, when the bond between the ligand and the protein is covalent, the standard procedure, which uses the tleap module of AMBER to assign the force field parameters, fails. In fact, in this last case, the molecular mechanics parameters of the ligand-enzyme covalent bond (in red in Figure 1) are missing. Table 1. Generated dihedral Amber parameters for compounds 1-4, using the Paramfit procedure. For other dihedrals of the molecules, General Amber Force Field (GAFF) wildcard parameters were used. The atom names are indicated as in the GAFF/GAFF2 and ff14SB force field.

Dihedral
Divider Moreover, improper dihedrals are also described in the same way [26,29], and the improper parameter, calculated through paramfit for our molecules, is reported in Table 1.
MD simulations. The above AMBER parameters, reported in Tables 1 and 2, were used to accomplish the MD simulations of BCCs 1-4 in complex with the AmpC βlactamase [24] in order to evaluate the conformational behavior of the single bonds around the Boron atom. Before starting MD simulations, BCCs 1-4, deprived of the serine portion (colored in blue in Figure 1), were covalently docked by GOLD [30] into the active site of E. coli AmpC β-lactamase [24], adopting a computational procedure previously reported by us [9]. As an example, the binding mode of the BCC 2 is reported in Figure 2. used to accomplish the MD simulations of BCCs 1-4 in complex with the AmpC -lactamase [24] in order to evaluate the conformational behavior of the single bonds around the Boron atom. Before starting MD simulations, BCCs 1-4, deprived of the serine portion (colored in blue in Figure 1), were covalently docked by GOLD [30] into the active site of E. coli AmpC -lactamase [24], adopting a computational procedure previously reported by us [9]. As an example, the binding mode of the BCC 2 is reported in Figure 2.   Table 3.  Table 3.  Conformational analysis of BCCs 1-4. With the aim of evaluating the accuracy of the new AMBER parameters developed by us, the full conformational analysis of BCCs 1-4 ( Figure 1) was performed through DFT calculations at the B3LYP/6-31G(d) level, using water as a solvent [31,32]. Successively, the torsional angles values τ1-τ4 in the located conformations of each compound (Table 4) were compared to those obtained by MD simulations ( Table 3). The 3D plots of the most significant conformations of BCCs 1-4 are shown in Figure 4.   Table 3. Average values of dihedral angles τ 1 -τ 4 for BCCs 1-4, which resulted from MD simulations.
Compound Conformational analysis of BCCs 1-4. With the aim of evaluating the accuracy of the new AMBER parameters developed by us, the full conformational analysis of BCCs 1-4 ( Figure 1) was performed through DFT calculations at the B3LYP/6-31G(d) level, using water as a solvent [31,32]. Successively, the torsional angles values τ 1 -τ 4 in the located conformations of each compound (Table 4) were compared to those obtained by MD simulations ( Table 3). The 3D plots of the most significant conformations of BCCs 1-4 are shown in Figure 4. Table 4. Geometrical features, relative energies, and equilibrium percentages in water of the located conformations of compounds 1-4. The attained results suggested that in the case of compound 1, the minimum energy conformation 1A shows geometrical preferences visited during the MD simulations. In fact, τ 1 , τ 2 , τ 3 , and τ 4 Table 4 and Figure 4).

∆E
For compound 2, the preferred geometry 2A, found during the DFT analysis, is the only significantly populated one (>90%). Once again, in the case of compound 3, the preferences determined through the MD simulations correspond to the conformational ones found in the minimum energy conformer 3A of the DFT analysis, which accounts for more than 50% of the population. In fact, τ 1 , τ 2 , τ 3 , and τ 4   The attained results suggested that in the case of compound 1, the minimum energy conformation 1A shows geometrical preferences visited during the MD simulations. In fact, τ1, τ2, τ3, and τ4 have values of 171°, 142°, −68°, 176° in 1A, fitting well with regions showing average values of 150°, 175°, −65°, 152°, for τ1, τ2, τ3, and τ4, respectively (see Table  4 and Figure 4).
Despite the good correspondence between the data attained by MD simulations and those of the DFT study, as explained above, all the free energy minima found by the conformational analysis studies were not fully visited during the MD simulations, considering the geometries populated by more than 10%. This fact might be due to the plausible steric clash between the ligand atoms and the residue atoms that can be found in the AmpC β-lactamases. Therefore, trying to solve this concern, MD simulations were accomplished on simplified systems constituted by the four investigated BCCs 1-4, in which dihedral angle fluctuations over the MD simulations time (25 ns) were observed and reported in  As expected, the conformational mobility of the boronic moieties significantly increased and new torsional angles were visited over the MD simulations, especially for dihedrals τ 3 and τ 4 , describing the orientation of hydroxyl groups. Moreover, despite the absence of AmpC residues and the potential creation of steric hindrance or hydrogen bonds limiting the conformational freedom of the moiety, in the case of dihedral τ 2 , only the trans orientation resulted to be visited for compounds 3 and 4.
These results represent a further confirmation of the good quality of the new determined parameters.
formational analysis studies were not fully visited during the MD simulations, considering the geometries populated by more than 10%. This fact might be due to the plausible steric clash between the ligand atoms and the residue atoms that can be found in the AmpC -lactamases. Therefore, trying to solve this concern, MD simulations were accomplished on simplified systems constituted by the four investigated BCCs 1-4, in which dihedral angle fluctuations over the MD simulations time (25 ns) were observed and reported in Figures 5-8.       As expected, the conformational mobility of the boronic moieties significantly creased and new torsional angles were visited over the MD simulations, especially dihedrals τ3 and τ4, describing the orientation of hydroxyl groups. Moreover, despite absence of AmpC residues and the potential creation of steric hindrance or hydro bonds limiting the conformational freedom of the moiety, in the case of dihedral τ2, the trans orientation resulted to be visited for compounds 3 and 4.
These results represent a further confirmation of the good quality of the new d mined parameters.

Materials and Methods
Paramfit procedure. Paramfit [23] is a program of the Amber package able to ge

Materials and Methods
Paramfit procedure. Paramfit [23] is a program of the Amber package able to generate or improve force field parameters. Preliminarily, RESP atomic charges were calculated for all studied compounds, using the Gaussian16 program package [33] at the B3LYP level of calculation with a 6-31G(d) basis set [34]. Four different steps were followed to determine the requested parameters.
Step 1. Firstly, the system setup was formulated. A topology file, with the parameters to be fit, was created. Moreover, a frcmod file (or force field modification) was prepared containing all parameters to be fit.
Step 2. At this point, different compound conformations were generated. In particular, a variety of molecular structures that sampled the conformational space, involving parameters that had to be fit, were systematically generated. The quality of parameter set was controlled to verify that structures adequately sampled the space. The structure set quality was then evaluated. Quantum calculations were conducted to determine the energies of the different geometries. In this step, either the energy or forces of each conformer at the quantum level of theory B3LYP/6-31G(d) [34], already used above, were calculated, finding a list of ab initio quantum energies. The final energy values were extracted into a quantum energy data file, showing the energy of each structure in the same order as the coordinate file.
Step 3. Quantum output files were processed. Using paramfit and the energies file, parameters for force field equations were determined. When fitting to calculated energies, Paramfit derives parameters performing the least squares. In this way, the program is able to regulate the parameters. This procedure allows us to minimize the least squares difference between the quantum energy of the starting geometries and the determined AMBER energy. To this aim, the following equation is applied: In this equation the different terms are: N, which is the number of geometries located for the molecule; E QM , which is the quantum energy evaluated through single-point calculations on the different conformations; and E MM , which is the calculated AMBER energy value for the same geometries. The constant term K, which depends on the molecule and set of input structures, was thus calculated. Contextually, R 2 was determined.
Step 4. Finally, parameters to fit were defined and a new frcmod file was prepared, before starting the MD simulations.
Docking and MD simulations. Initially, the BCCs 1-4 were deprived of the structural moieties mimicking the serine residue. Covalent docking calculations were then performed by GOLD (version 2021.3, CCDC), using as a target the AmpC β-lactamases of E. coli (PDB accession code 1KE3) [24]. The ligand binding site was defined selecting residues included in a sphere with a radius of 12.0 Å from the side chain oxygen atom of Ser64 [5], and all residues were kept rigid in these calculations. The water molecules found in the X-ray were removed to avoid any steric clashes during this preliminary step of simulations. The pose acquiring the highest GoldScore was selected to generate the ligand-enzyme complexes for MD simulations. The RESP atom charges were assigned to the ligands by calculations at QM level using Gaussian16 [33] at the DFT/B3LYP/6-31G(d) level of theory [25,34]. The molecular mechanics parameters needed for MD simulations were assigned to the ligand-protein covalent complexes by the antechamber module of AMBER18 [29]: the ff14SB force fields [35] and TIP3P model [36] were used to represent the enzyme and the water solvent, respectively. In this step, a solvent box with a minimum distance of 15 Å from the AmpC surface was built for each complex. Prior to starting the MD simulations production runs, a minimization of the bulk solvent molecules was performed by applying a gradient criterion convergence of 0.2 kcal mol −1 Å −1 . The entire systems were then optimized by a convergence criterion of 0.0001 kcal mol −1 Å −1 . Successively, the whole systems were equilibrated, gradually increasing the temperature from 0 to 300 K over 60 ps of MD simulations in isocore conditions (NVT). Finally, production runs of 50 ns were accomplished for each ligand/enzyme complex in isothermal-isobaric ensemble at 300 K, with a 1 fs time-step (NPT). In these simulations, the systems were performed in periodic boundary conditions, while the Van der Waals and short-range electrostatic interactions were estimated within an 8 Å cutoff. The attained trajectory frames were visually inspected by VMD software [37]. MD simulations on compounds 1-4 (without the protein environment) were accomplished adopting the protocol here reported for the ligand/AmpC complexes.
DFT Calculations. All the calculations were carried out by using the GAUSSIAN16 program package [33]. The conformational space of compounds 1-4 was explored through the optimization of all the possible starting geometries using the DFT approach at the B3LYP level with the 6-31G(d) basis set [38,39]. To take into account the influence of the solvent, optimizations were conducted in water, using the polarizable continuum model (PCM) [40]. All the degrees of conformational freedom were considered. The conformational preferences of the single bonds of the molecules were determined through the evaluation of the three different gauche, mgauche and anti orientations, combining them and locating different geometries. The percentage contribution of each optimized conformation to the overall population was determined at 298 K through the Boltzmann equation. Vibrational frequencies were computed at the same level of theory to verify that the optimized structures were minima.

Conclusions
In this paper, we have developed the Amber force field parameters for phenyl-, benzyl-, benzylamino-, and methylamino-boronates, a group of BCCs capable of creating covalent complexes with the reactive amino acids (like serine, threonine, or cysteine) located in the active sites of enzymes like β-lactamases. In fact, the computational design of new compounds of this class required an accurate molecular mechanics parametrization that is only available for a limited number of aryl-, alkyl-boronic acids, and boronate esters. By our calculations, we have developed the force field parameters for selected BCCs 1-4, testing them by performing MD simulations on ligand/β-lactamases covalent complexes (obtained by docking studies) and also by simulating them covalently bound with a putative serine residue. The robustness of our results was highlighted by the comparison of the values populated by torsional angles over the trajectory frames with those determined through a complete DFT conformational analysis. The good correspondence achieved the goal to create a library of parameters to perform MD simulations on new, selected BCCs that have important clinical applications and attracted the interest of medicinal and organic chemists.