|Home | About | Journals | Submit | Contact Us | Français|
The accuracy of molecular dynamics (MD) simulations is limited by the availability of parameters for the molecular system of interest. In most force fields parameters of common chemical groups are already present. With the development of novel small organic molecules as probes to study biological system more chemical groups require parameterization. An azide group is often used in studies of biological systems yet computational studies are impeded by the lack of parameters. In this paper we present a set of molecular mechanics (MM) parameters for aromatic and aliphatic azido groups and their application in MD simulations of a photoaffinity probe currently used in our laboratory for mapping the binding modes available in the active site of histone deacetylases. The parameters were developed for the Generalized Amber Force Field (GAFF) using density functional theory calculations (DFT) at B3LYP 6-311G(d) level. The parameters were validated by geometry optimization and MD simulations.
Molecular dynamic (MD) simulations are the method of choice to study small and macromolecular systems as they provide valuable information on their dynamic behavior at the atomic level. Combining Molecular Mechanics (MM) force field parameters for macro- and small molecules in MM/MD calculations is usually a challenging task as historically MM force fields were designed for either one of these two systems. Until recently, simulating the interactions between small chemical compounds and proteins, a typical task in studying drug-receptor interactions, meant having to re-parameterize the original force field parameters developed for small organic compounds for the force field used with the protein. To address this issue a general Amber force field (GAFF) for common chemical groups was designed to be compatible with the popular series of Amber force field for proteins and nucleic acids, which significantly accelerated the setup of ligand-protein systems for MM/MD calculations. The number of chemical groups parameterized for GAFF is gradually increasing since more and more small molecule ligands bound to their macromolecular targets require MM/MD simulations. We recently developed a series of photolabeling ligands containing an N3-aryl (N3-Ar) and N3-alkyl (N3-Alk) groups to probe the binding site of histone deacetylases , a family of enzymes that are targeted in design of therapeutics for cancer and neurological diseases [2–7]. Although the Amber parameters for an N3-Alk group have been previously reported as a part of the parameterization of zidovudine (AZT) , there are no parameters for an N3-Ar, which are required to model Photoaffinity Probe 1 (Fig. 1) that is currently used in our laboratory .
In this study we have determined a set of MM parameters for the N3-Ar group and validated them using simulations with Compound 1. For consistency with the new aromatic N3-Ar parameters, we also re-calculated parameters for the aliphatic N3-Alk group. The results of these studies are presented in this paper.
High transferability of the MM parameters in GAFF is achieved by describing the most common chemical groups using a limited set of atom types that are devised to be as general as possible and an energy function (Eq. 1) that does not include cross-terms to describe the molecular system [9, 10]. The MM parameters were calculated by fitting the QM PES using standard procedures. The regression parameters correspond to the force constants Kr and Kθ of the bonds and angles PES, respectively, and to the height of the potential Vn in the case of the dihedrals PES (Table 1).
The model for Compound 1 (Fig. 1) used for the quantum mechanical (QM) calculations was built using GaussView . All QM calculations were performed using Gaussian 03 . The geometry was optimized using a two-step procedure, first at the HF/6-31G(d) and then at the B3LYP/6-31G(d) levels of theory. The restrained electrostatic potential (RESP) charge approach was used to derive atomic charges. First the high-density ESP calculation was performed at the B3LYP/6-31G(d) using the Merz-Kollman scheme. Then the charges were derived using Antechamber, which performed a two-stage charge fitting with default hyperbolic restraint parameters.
The other type of non-bonded parameters, the Van der Waals parameters, were taken from GAFF directly since they are known to be transferable as they are mainly dependant on the number of electrons of an atom rather than their chemical environment.
Two smaller systems derived from the optimized structure were used to calculate the missing parameters. An azido group attached to benzene and an azido group attached to the methyl group of toluene were used as model systems for the aromatic and the aliphatic azido groups and are shown in Figure 2a and Figure 2b, respectively. Potential energy surface scans (PES) of the missing bonds, angles and dihedrals parameters were performed along the internal coordinates of the parameter of interest. All other parameters were kept fixed since the GAFF energy function that we want to fit does not contain cross-term potentials. The initial geometry optimization and subsequent single point energy calculations performed during the PES scans were carried out with B3LYP using the 6-311+G(d) basis set. Bond scans were carried out using an increment of 0.002 Å for a total scan length of 0.04 Å centered on the equilibrium value. For bond angles, a 1° or 2° increment was used with a total scan range varying from 15° to 30°.
Energy minimizations and MD simulations were performed using the Gromacs software package (version 3.3.3) [13–16]. Compound 1 was parameterized using GAFF [10, 17] as ported to Gromacs  and the newly developed parameters calculated for the N3-Ar and N3-Alk groups. TIP3P explicit solvent  was used to solvate the system. The simulations were run in a dodecahedral box containing 747 water molecules and using periodic boundaries conditions. The simulation box was created by placing the edge of the box at a minimum distance of 8 Å of the solute. The bonds in Compound 1 were constrained using the LINCS algorithm  while the bonds and the angles in the water molecules were constrained using the SETTLE algorithm . A time step of 2 fs was used for integrating the equations of motion. The system was simulated in an NTP ensemble at a fixed temperature of 300 K and a fixed pressure of 1 bar. Temperature was kept constant using a Berendsen thermostat  with a coupling time of 0.1 ps and pressure was controlled by a weak coupling to a reference pressure  with a coupling time of 1 ps and an isothermal compressibility of 4.6 10−5 bar−1. Non-bonded interactions were calculated using a cut-off of 9 Å. Long-range Coulombic interactions were evaluated using Particle Mesh Ewald (PME) algorithm  with an interpolation order of 4 and Fourier spacing of 1.2 Å. The system was energy minimized and equilibrated for 200 ps after initial starting velocities were randomly assigned from a Maxwell distribution centered at 300 K. Finally the production simulation was run for 20 ns.
The parameters resulting from the fitting of the QM PES are given in Table 1 and the statistical details of the calculations, such as the number of points in each data set, correlation coefficient r and F-values, are given in the supplemental material. The F-value calculated for each model (bonds, angles and torsions) is much higher than its respective critical F-value at the 1% significance level. F-values range from 477 to 19827 whereas the critical F-values are all below 12 indicating that all the models obtained from the fitting can be used to explain the data. The quality of the model, i.e. how well it correlates with the data, is given by the correlation coefficient r, which was higher than 0.990 for all the parameters (including dihedrals) with the exception of the torsion ca-c3-Ni-Nd for which r was equal to 0.983. These results indicate overall excellent correlation between the fitted MM force field potentials and the DFT-based results.
Experimental studies on aliphatic azides [24–27] indicate that the N-N-N angle is close to 180° suggesting that the middle nitrogen atom is best described as a mix between sp and to a much lesser extent sp2 hybridization states. Because of this relatively unique electron distribution in the azido group none of the nitrogen atoms can be correctly described by the existing GAFF nitrogen atom types.
We started the development of the parameters by deciding whether the aryl and alkyl N3 groups would require two different atom types for the nitrogen atom directly attached to the aromatic or aliphatic carbon atom. Interactions between the azide nitrogen atoms N1, N2 and N3, where N1, N2 and N3 correspond to the first, middle and last nitrogen atoms (Fig. 2a,b) of the azido groups respectively, were examined using B3LYP and the 6-311+G(d) basis set. The results obtained for the bonds N1-N2, N2-N3, and the angle N1-N2-N3 (Table 1) indicate that the force constants, the equilibrium bond length and the angle values are almost identical for both the aliphatic and aromatic substituents with differences smaller than 1%. Thus, the same sets of atom types can be used for aromatic and aliphatic azido groups. For comparison with the previous parameterization publication we kept the same atom type names used in , where Ni, Nd and Ne correspond to N1, N2 and N3 respectively.
With the azido group geometry and force characteristics independent of its substituent, the difference between the parameters covering the interface between the azido group and the substituent in N3-Ar and N3-Alk systems still required further evaluation. The differences observed in the bonds properties between ca-Ni (aromatic) and c3-Ni (aliphatic), in angle properties between ca-Ni-Nd and c3-Ni-Nd, ca-ca-Ni and ca-c3-Ni and in dihedral properties clearly indicate that different parameters should be used to model N3-Ar and N3-Alk groups. This is notably the case for the bonds ca-Ni and c3-Ni where the bond with the aromatic carbon is much shorter (1.421 Å) and stronger (380.9 kcal mol−1 Å−2) than its aliphatic equivalent with 1.491 Å and 293.5 kcal mol−1 Å−2, respectively. The bond angles values for x-N1-N2 and x-x-N1 are similar in N3-Ar and N3-Alk groups, although, as expected, the force constants are greater for the parameters in the aromatic N3-Ar group due to electron conjugation between its N3 and aryl portions.
It should noted that for the most part the parameters for the aliphatic azido group are similar to those described in  with the main difference in the length of the bond Ni-Nd being shorter in our study, 1.230 Å vs. 1.340 Å. The bond length obtained in our calculations is in agreement with the experimental values reported to be 1.243 Å [24–26] for hydrazoic acid, 1.216 Å  and 1.240 Å  for azidomethane, and 1.229 Å  for azidoethane. When comparing these values with the equivalent parameters in GAFF, the N1-N2 bond characteristics appears as a combination of the characteristics of the bonds represented by GAFF atom types n2-n1, n3-n1 and n2-n2, where n3 corresponds to sp3 hybridized N atom, n2 to sp2, and n1 to sp. This observation confirms that the middle nitrogen atom N2 is a mix of sp and sp2 hybridization states.
Finally, since the main goal of this work was to obtain MM parameters to study dynamics of bis-azido Probe 1, we also needed to determine the ligand atom charges. To calculate a reliable charge distribution in 1 finding the correct protonation state of the hydroxamic acid is critical. With a pKa of 9.4 the hydroxamic acid is protonated in water at pH 7.4. However there is evidence that it might not be the case for hydroxamate-based ligands when bound to HDAC active site. It was shown that the pKa of hydroxamic acids decreases by ~3.3 log units upon forming a complex with the Zinc atom in the HDAC catalytic site . Therefore, the charges for the deprotonated state of Probe 1 were calculated. The schematic representation of Compound 1 is shown in Figure 1 and the corresponding atom types and partial charges are given in Table 2.
The high internal charges observed for the nitrogen atoms of the azido group are consistent with the main resonance structures of the azido group (Fig. 3) where the middle nitrogen carries a positive charge and either of the terminal nitrogen atoms carries a negative charge. The charges obtained for the deprotonated Probe 1 were found to be virtually identical to those obtained for azidobenzene and azidoethane demonstrating that the charge distribution does not depend on the overall charge of the molecules but instead is specific to the azido group. The detailed results of these calculations are provided in the supplemental material.
To check the accuracy of the calculated parameters, we compared the structures obtained from the geometry optimization at the DFT level with the geometry obtained after energy minimization at the MM level using a tight convergence criterion of 1 kJ.mol−1.nm−1 for the maximum force. The root mean square deviation between the heavy atoms of the two structures after least square fit was only 0.38 Å, indicating that the QM geometry could accurately be reproduced at the MM level. An overlay of Compound 1 optimized at the MM and QM levels is shown in Figure 4.
The parameters were further checked by performing the MD simulations with the Gromacs software package of Ligand 1 in explicit solvent during a 20 ns production run. Figure 5 shows that the RMSD of each individual azide group is low throughout the course of the simulation with fluctuation around 0.02 Å, indicating that bonds and angles remain close to their equilibrium values.
Figure 7 shows that the distribution of the dihedral angles of the aromatic azido group is consistent with the dihedral potentials previously calculated for the model systems and shown in Figure 6. In the case of the dihedral CA-CA-Ni-Nd the two minima observed at −90° and 90° in the dihedral distribution shown in Figure 7a correspond to the maxima of the dihedral potential shown in Figure 6a. Similarly the two maxima at 0° and 180° in the dihedral distribution match the minima in the dihedral potential. In the case of the dihedral CA-Ni-Nd-Ne there is only one maximum at 180° and one minimum at 0° in the dihedral distribution (Fig. 7b) and both are consistent with the minimum and maximum observed in the dihedral potential (Fig. 6b). The dihedral distribution for the dihedral CA-CA-Ni-Nd-Ne shows that the most dominant conformations of the aromatic azido group correspond to a coplanar arrangement with the aromatic ring. This effect is expected to be due to π electron delocalization between the π orbitals of the phenyl ring and of the azide group and is represented in the MM model by the potential well. There is a strong energy penalty for leaving the planar arrangement however the potential still allows interconversion at 300 K between the two dihedral minima since they both appear equally populated.
Overall, the results suggest that the parameters are well suited for reproducing the geometry of the azido groups during energy minimizations. An analysis of the MD data for Ligand 1 shows that the potential energy and its components oscillate around their average values and remain stable during the extensive MD simulation in water suggesting that the parameters are suitable for MD simulation studies.
The GAFF parameters for N3-Ar and N3-Alk were calculated. The quality of parameters was confirmed by accurately reproducing the QM PES. Our study indicates that the 3 new nitrogen atom types are likely to be sufficient to describe both the N3-Ar and N3-Alk groups in small organic molecules in MM force fields. Further analysis has shown that two different sets of bond, angle, and torsion parameters for the group of atoms consisting of the nitrogen atoms in the azido group and the adjacent carbon atoms should be used to describe these systems. As expected, the force constants for the bonds and angles of the N3-Ar azido group are greater than those for N3-Alk group because of the effect of delocalization in the former group. The parameters were tested in an MD simulation of the hydroxamate HDAC Probe 1 in water. The system was found to be stable indicating that the parameters are suitable for MD simulations of and should be useful for computational studies of ligands containing either or both N3-Ar and N3-Alk groups.
We would like to thank Dr. Michael Brünsteiner for critical analysis of the manuscript. This study was in part funded by the Breast Cancer Congressionally Directed Research Program of Department of Defense Idea Award BC051554 and by the National Cancer Institute/NIH grant 1R01 CA131970-01A1. We also thank OpenEye Scientific Software for providing academic license for modeling software.