16590
J. Phys. Chem. B 2009, 113, 16590–16595
Reoptimization of the AMBER Force Field Parameters for Peptide Bond (Omega) Torsions Using Accelerated Molecular Dynamics Urmi Doshi and Donald Hamelberg* Department of Chemistry and The Center for Biotechnology and Drug Design, Georgia State UniVersity, Atlanta, Georgia 30302-4098 ReceiVed: July 31, 2009; ReVised Manuscript ReceiVed: October 22, 2009
Improving the accuracy of molecular mechanics force field parameters for atomistic simulations of proteins and nucleic acids has been an ongoing effort. The availability of computer power and improved methodologies for conformational sampling has allowed the assessment of these parameters by comparing the free energies calculated from molecular dynamic (MD) simulations and those measured from thermodynamic experiments. Here, we focus on testing and optimizing the AMBER force field parameters for the ω dihedral, which represents rotation around the peptide bond of proteins. Due to the very slow isomerization rate of the peptide bond, it is not possible to sample the phase space with standard MD simulations. We therefore employed an accelerated MD method in explicit water in which the original Hamiltonian is modified to speed up conformational sampling and the correct canonical distribution is recaptured. Using well-studied model systems for the peptide and peptidyl prolyl bonds, we discovered that the AMBER ω dihedral parameters underestimated experimentally measured activation free energy barriers for cis/trans conversion as well as failed to reproduce the free energy difference between the two isomers. We reoptimized the original AMBER ω dihedral parameters and further validated their transferability on several experimentally studied dipeptides. The revised set of parameters successfully reproduced the cis/trans equilibria and free energy barriers within experimental and simulation errors. We also investigated the structures of the transition state and cis/trans isomers of prolyl peptide bonds in terms of pyramidality, a measure of the puckering of the prolyl ring. We observed, as expected from quantum mechanical studies, significant bidirectional, out-of-plane motions of prolyl nitrogen in the transition state. Introduction Computational simulations have become an absolutely essential tool complementing experimental studies for the investigation of structures, thermodynamics, and dynamics of biomolecular systems. Due to large sizes of biological macromolecules, there are limitations in using a quantum mechanical approach for the evaluation of the potential energy function. Therefore, physics-based molecular mechanics force fields, which have a simple potential energy function and empirically derived parameters, are generally employed for MD simulations of biomolecules. Since the reliability of the results obtained from such theoretical studies depends on how accurately structures are represented, energies are calculated, and conformational states are sampled, it is not surprising that such empirical force fields constantly undergo fine-tuning in their functional form and parameters. The AMBER force field, which was originally developed more than two decades ago, was parametrized essentially from gas-phase quantum mechanical and empirical data on short, model peptides and small, organic compounds.1,2 Since then, it has received several updates based on the availability of higher level ab initio theories, more accurate solvent models, and experimental data.3–11 The potential energy functional of AMBER includes terms for bond stretching, angle bending, torsion rotating, van der Waals forces, and electrostatic interactions. The total torsional energy is expressed as the Pitzer potential,12 a Fourier series term given by ∑dihedrals(Vn/2)[1 + cos(nφ - γ)], * Corresponding author. Phone: 404-413-5564. Fax: 404-413-5551. E-mail:
[email protected].
where Vn, n, φ, and γ are the dihedral force constants, periodicity, torsional angle, and phase angle, respectively. Due to the pseudo-double bond character of the peptide bond, its rotation is restricted such that the ω dihedral occurs predominantly in the trans (ω centered around 180°) and cis (ω centered around 0°) conformations. The first two terms in the Fourier series are sufficient to represent the peptide bond torsional energy profile. The one-fold rotational term (n ) 1, γ ) 0°) with a maximum at 0° reproduces the energy difference between the trans and cis isomers, whereas the second term (n ) 2, γ ) 180°) reproduces the 2-fold symmetry of the peptide bond torsional potential with minima at 180° (and -180° due to periodicity) and 0°. Although the V2 constant essentially determines the barriers to rotation about the C-N bond, the V1 constant is responsible for the correct cis/trans equilibria. In force fields, the torsional parameters V1 and V2 should be specified for each set of four atoms (i.e., X-Ci-Ni+1-X where X can be C, O, or H) that define the ω dihedral around the C-N bond. For non-proline-containing peptide bonds (i.e., Xaa-non-Pro) both V1 and V2 are required to obtain the right free energy difference between the cis and the trans isomers as well as the barrier heights for rotation around the C-N bond. However, V1 is not used in defining the ω dihedral of prolinecontaining peptide bonds (i.e., Xaa-Pro) because the cis and the trans isomers are almost iso-energetic in terms of free energy; that is, ∆Gc-t is less than 1 kcal/mol and can be reproduced by the 2-fold rotational term, without the need for a V1 term.
10.1021/jp907388m 2009 American Chemical Society Published on Web 11/25/2009
AMBER Parameters for Peptide Bond Torsions
Figure 1. Structures of N-acetyl-N′-methylprolineamide (Ace-ProNHMe) and N-methylacetamide (NMA). The ω dihedral is defined by atoms C′-C-N-CR in Ace-Pro-NHMe and C′-C-N-C′ in NMA where C′ is methyl carbon. Cis-trans isomerization of the ω dihedral is simulated in this study, whereas ω′, another peptide bond dihedral in Ace-Pro-NHMe is not considered. The free energy profiles in Figure 2 are projections along the ω dihedral degree of freedom. For AcePro-NHMe, the four paths for the ω dihedral are (i) C′-C-N-CR, (ii) O-C-N-CR, (iii) O-C-N-Cδ, and (iv) C′-C-N-Cδ. In the case of NMA, the four paths along the ω dihedral are (i) C′-C-N-C′, (ii) H-N-C-O, (iii) O-C-N-C′, and (iv) H-N-C-C′. In AMBER, the dihedral constant V2 is for the general X-C-N-X torsion, whereas V1 and V2′ are parameters specific for the H-N-C-O torsion.
The degree of steric conflict of the two CR atoms flanking the peptide bond is much more in the cis conformation, resulting in only 5-6% occurrence of cis peptide bonds in protein structures.13–15 From spectroscopic studies on N-methylacetamide (NMA) (Figure 1), the nearest analog of peptide bonds in proteins, the trans form is found to be almost 30-45 times more favorable than the cis form.16,17 Due to the small population of the cis form, it has been difficult to characterize the kinetics of cis-trans isomerization experimentally, as a result of which very limited data on the free energy barriers separating the two isomers is available. The force constant V2 ) 20 kcal/mol of AMBER1 was derived from rotational barriers estimated from approximate molecular orbital calculations18,19 and normal mode analysis of vibrational spectra obtained from infrared studies20,21 on the prototypic NMA. To represent an effective torsional barrier of 20 kcal/mol, each of the four possible combinations X-Ci-Ni+1-X (Figure 1) associated with the C-N bond was given a value of 20/4 ) 5 kcal/mol for V2 in the original AMBER force field.1 A 2-fold term with V2′/2 ) 2.5 kcal/mol and a one-fold rotational term with V1/2 ) 0.65 kcal/mol were included in the later Weiner et al. force field2 specifically for the H-N-C-O torsion. It is not clear from the literature on AMBER force fields how the dihedral constants for the H-N-C-O torsion were obtained. Since other nonbonded van der Waals and electrostatic parameters were modified in the ff94 version of AMBER,4 V1/2 was corrected to 2.0 kcal/mol so that the potential energy difference between the cis and the trans isomers of NMA in vacuo remained at ∼2.3 kcal/mol. Although the parameters of the potential energy functional are derived from empirical data, assessment of their accuracy can be performed only by direct comparison of both the population of conformers and the free energy barriers for interconversion between conformers obtained from experiments and MD simulations. In the past, such exercises have revealed the limitations of AMBER backbone dihedral parameters in estimating the correct stability of R-helical and β-strand conformations.6 However, this requires calculation of free energy profiles from MD simulations, which in the case of the ω dihedral has not been attempted. The reason is the slow isomerization rate of the ω dihedral (on the order of several seconds) that would result in rare or no transitions at all between the cis and trans wells in normal MD simulations and lead to
J. Phys. Chem. B, Vol. 113, No. 52, 2009 16591 poor sampling of cis/trans conformers. There would be loss of ergodicity, and even with very long simulations, it would no longer be safe to assume that a Boltzmann distribution is achieved; that is, convergence would not be reached. Hence, in the present study, we implement the accelerated MD approach to generate free energy profiles for the rotation of the ω dihedral and reoptimize peptide bond torsional parameters. The accelerated MD method was developed earlier to increase the escape rates from potential energy minima22,23 and freely explore the energy landscape without any conformational bias. It allows for simulations to be carried out in atomistic detail on a modified landscape and recovery of equilibrium properties on the original landscape by reweighting the distribution. We employ NMA and N-acetyl-N′-methylprolineamide (Ace-Pro-NHMe) as models of a generic peptide bond and a peptidyl-prolyl bond, respectively (Figure 1). Since the N-methyl groups of NMA are magnetically unequivalent in the cis and trans isomers, it was possible to obtain rate constants and thereby free energy barriers of cis/trans conversion from NMR line shape analysis of N-methyl signals.16 From these studies, a barrier height of 21.3 kcal/mol for the trans f cis conversion and a free energy difference of 2.5 kcal/mol between the cis and trans isomers were obtained. In proline-containing peptides, the trans and the cis isomers are energetically almost similar (i.e., the trans isomer is favored only by less than 1 kcal/mol over the cis form24–27). Due to the steric hindrance between the CR of the preceding residue and Cδ of the proline ring, the trans isomer is relatively destabilized. This drastically increases the occurrence of the cis isomer up to 40% in Xaa-Pro peptides, as compared to ∼1.5-2% in Xaa-non-Pro peptides.16,17 For this reason, a considerable number of experimental and computational studies have been carried out on the cis-trans isomerization of several Xaa-Pro peptides28,29 and their analog Ace-Pro-NHMe.24,26,30–35 In addition, cis-trans isomerization of Xaa-Pro motifs manifests as the slow kinetic phase and, hence, the rate-determining step for the folding of many proteins.36 Its study is also of interest while investigating peptidyl-prolyl cis-trans isomerases (PPIases), a class of enzymes important for regulating various biological processes.37,38 Therefore, there is a need to have force field parameters that can reproduce experimental results. Since the ω torsional energy of Ace-Pro-NHMe is determined only by the second rotational term (Figure 1) (the one-fold term is applicable only for the H-N-C-O torsion), we first optimize the dihedral constant V2 by matching the cis/trans free energy barriers estimated from experiments and accelerated MD simulations in explicit water. We then use the modified V2 to calibrate the dihedral constants (V1 and V2′) for the H-N-C-O torsion by carrying out accelerated MD of NMA, also in explicit water. Furthermore, we also test the transferability of these reoptimized parameters on Xaa-Pro and Xaa-non-Pro dipeptides. Computational Methods NMA and Ace-Pro-NHMe were used as models for a generic non-proline peptide bond and a Xaa-Pro peptide bond, respectively. The parameters governing the energy barrier to the rotation of the peptide bond and the free energy difference between the trans and the cis forms were revised using the above model systems. We then tested the transferability of the modified parameters on dipeptides such as alanyl-phenylalanine (AlaPhe), representing Xaa-non-Pro, and alanyl-proline (Ala-Pro) and phenylalanyl-proline (Phe-Pro), representing Xaa-Pro, in their zwitterionic forms. In addition, we also investigated the cis-trans isomerization of Ala-Pro in the cationic form (i.e., NH3+-Ala-Pro-COOH modeled as NH3+-Ala-Pro-CONH2).
16592
J. Phys. Chem. B, Vol. 113, No. 52, 2009
Doshi and Hamelberg
Figure 2. Average reweighted free energy profiles as a function of the ω dihedral for (A) Ace-Pro-NHMe and (B) NMA obtained from five independent accelerated MD runs using original AMBER parameters (blue) and modified parameters (red) for peptide bond torsion dihedral constants. The free energy profiles are normalized with respect to the global trans minimum.
All of the MD simulations were carried out using the AMBER 9.0 suite of programs with the modified version (ff99SB)6 of the parm99 all-atom force field.4 The starting structures and input files were generated using the xleap program. All the model systems and test dipeptides were solvated in cubic periodic boxes with the TIP3P (transferable intermolecular potential 3P) water model.39 In another set of simulations, we also used the SPC/E water model40 for solvating the test dipeptides to evaluate the effect of the solvent model on the barrier heights and cis/trans equilibria. Each periodic box had an approximate dimension of 26 × 26 × 26 Å3 and was filled with close to 500 water molecules that were placed up to 9 Å away from the solute atoms. The starting structures were equilibrated in three steps as follows: An initial minimization of 5000 steps was performed with a steepest descent algorithm to relax the solvent while the solute atoms were kept fixed with a harmonic potential (restraining force constant of 50 kcal/mol · Å2). This was followed by a 100 ps MD simulation in which the restraining force constant was reduced to 25 kcal/mol · Å2. The last step of equilibration was also a 100 ps MD run in which both the solute and the water molecules were allowed to relax. These steps adjusted the density of the system to ∼0.986 g/mL. The production MD runs were carried out with the sander or the pmemd module (or both) of AMBER9.0 that was modified for the implementation of accelerated MD. The equilibration and production simulations were carried out in the NPT (fixed number of atoms, pressure, and temperature) ensemble at 1 bar and 300 K. During equilibration, the temperature was controlled by coupling to an external heat bath with a relaxation time of 0.2 ps, and for the production runs, temperature was regulated with a Langevin thermostat (collision frequency ) 1 ps-1). The system was maintained at the reference pressure by coupling to an external pressure bath using isotropic position scaling (coupling constant ) 1 ps). Newton’s equations of motion were integrated with a time step of 2 fs. All bonds involving hydrogen atoms were constrained using the SHAKE41 algorithm with a tolerance of 1 × 10-4. The nonbonded pair list was updated every 10 timesteps and a cutoff of 9 Å was used for the calculation of short-range nonbonded interactions. Long-range interactions were treated with the particle mesh Ewald summation.42 To study the slow process of cis-trans isomerization of the peptide bonds, we performed accelerated MD23 by adding a continuous nonnegative bias potential ∆V(r) (where ∆V(r) ) (E - V(r)2)/(R +(E - V(r))) only to the torsional potential V(r) whenever V(r) is less than a predefined boost energy, E. Simulations were performed on the original potential V(r) when V(r) > E; that is, when ∆V(r) ) 0. E and R are fixed throughout
TABLE 1: Parameters for Accelerated Molecular Dynamics of Model Systems and Test Peptides model/test peptide
av total dihedral energy (kcal/mol)
boost energy, E(kcal/mol)
R (kcal/mol)
Ace-Pro-NHMe Ala-Pro Phe-Pro NMA Ala-Phe
19.68 14.98 17.66 3.50 9.33
50 50 60 30 50
8 8 12 10 8
the simulations and dictate the extent of acceleration. E was set to ∼30-40 kcal/mol above the average total dihedral energy for each model and test peptide (see Table 1). Each configuration obtained from accelerated MD runs was reweighted by eβ∆V(r), and the probability density function of the ω dihedral was calculated using Gaussian kernel density estimation (eq 1) with a bandwidth (h) of 1.0. From the free energy profiles as a function of the ω dihedral angle (eq 2), the barriers to cis-trans isomerization and the free energy differences between the cis and the trans forms were estimated. For Ace-Pro-NHMe, accelerated MD simulations were repeated with several values of V2 and default values for the rest of the AMBER force field parameters. The value that gave the best agreement with the experimentally obtained free energy barrier for the cis-trans isomerization of the peptidyl-prolyl bond of Ace-Pro-NHMe was chosen. Using the modified V2, accelerated MD of NMA was performed to tune V1 and V2′ against experimental cis/trans equilibria and free energy barriers to rotation of the ω bond. This modified set of parameters (V1, V2, and V2′) was then used for the accelerated MD of test peptides. For each model and test peptide, the simulations were carried out for 300 ns in each of the five individual MD runs, for a total of 1.5 µs. Each time the simulations were restarted, a distinct random number seed was assigned to generate pseudorandom forces in the Langevin thermostat. The barrier heights and free energy differences between cis and trans isomers listed in Table 2 are average values from five independent MD runs. In addition, the errors reported in Table 2 are obtained from standard deviations across five independent MD runs. N
P(ω) )
∑
1 -1/2 1 e Nh i)1 √2π
(
ω - ωi
F(ω) ) -RT ln P(ω)
h
)
2
(1) (2)
Results and Discussion Reparametrization of ω Torsional Parameters Using AcePro-NHMe and NMA As Model Systems. Ace-Pro-NHMe is the simplest model for studying the isomerization of the peptide
AMBER Parameters for Peptide Bond Torsions
J. Phys. Chem. B, Vol. 113, No. 52, 2009 16593
TABLE 2: Comparison of Free Energy Barriers to ω Dihedral Rotation and Free Energy Difference between the Cis and Trans Isomers from Experiments and Accelerated MD experimental (kcal/mol) model/test peptide Ace-Pro-NHMe Ala-Pro Phe-Pro
NMA Ala-Phe
# ∆Gtransfcis b
20.40 19.80c 18.43d 20.58 ( 0.22e 18.88d 20.46 ( 0.24e 21.60h 19.10i 21.84j 21.30 ( 0.30f 21.68 ( 0.48g
from accelerated MD simulations (kcal/mol)
# ∆Gcisftrans
19.80 18.16d 20.39 ( 0.22e 19.60d 21.68 ( 0.24e 22.20h 19.65i 22.50j 18.80 ( 0.30f 18.12 ( 0.38g
∆Gcis-trans 0.57 ( 0.07
b
0.27d 0.19e -0.72d -1.22e -0.60h -0.55i -0.65j 2.50 ( 0.30f 3.35 ( 0.38g
# ∆Gtransfcis
# ∆Gcisftrans
∆Gcis-trans
20.4 ( 0.3 19.2 ( 0.3
19.3 ( 0.3 20.0 ( 0.3
1.1 ( 0.1 -0.9 ( 0.1
19.3 ( 0.7
21.2 ( 0.6
-1.9 ( 0.7
21.1 ( 0.0 20.8 ( 0.4
18.6 ( 0.1 18.5 ( 0.1
2.5 ( 0.1 2.3 ( 0.4
a Negative numbers indicate that the stability of the cis isomer is more than the trans one. For experimental numbers the references are quoted in footnotes b-h. In addition, the conditions and methods of experimental studies are mentioned. b Ref 24, 1H NMR, 298 K, D2O, pD 7.0 c Ref 42, pH jump to final pH 3.5, Activation energy obtained from Arrhenius plot of observed relaxation time vs inverse temperature. d Ref 43, 1H NMR, 298 K, D2O, pD 7.0. e Ref 27, dynamic capillary electrophoresis, 298 K, H2O, pH 9.5. f Ref 16, NMR, 333K. g Ref 45, 1H NMR, 13 C- NMR, 298 K, H2O:D2O mixture, pH 5.9. h Ref 44, D2O, pH 8.4, relaxation after ad hoc dissolution probed by capillary zone electrophoresis. i Ref 44, relaxation after in-column incubation probed by capillary zone electrophoresis, pH 8.4. j Ref 44, relaxation after ad hoc dissolution probed by 1H NMR, 283 K.
bond preceding proline in proline-containing peptides. Due to the pyrrolidine ring, the amide nitrogen is disubstituted, and the different combinations of four atoms that define the ω dihedral do not involve the H-N-C-O path (Figure 1). Therefore, the ω torsional potential is described only by the 2-fold rotational term, and the isomerization barrier is purely defined by V2. Hence, the one-dimensional free energy profiles along the ω dihedral obtained from accelerated MD simulations of Ace-Pro-NHMe help us in investigating the sheer contribution of V2 to the free energy barrier for cis/trans conversion. Figure 2A shows the free energy profile generated after accelerated MD simulations of Ace-Pro-NHMe with the original AMBER parameters for ω torsion; that is, V2 ) 20 kcal/mol. The default AMBER parameters yielded a ∆G#tfc ) 14.8 ( 0.1 kcal/mol that is underestimated from those obtained from experiments by ∼6 kcal/mol (Table 2, Figure 2A). This clearly showed that the ω dihedral constants of AMBER needed reparametrization. When V2 was modified to 28 kcal/mol, the free energy barrier height ∆G#tfc agreed remarkably with experimental estimates (Table 2). The free energy difference between the two isomers was also quite close to experimental ∆Gc-t, given the average error of 0.5 kcal/mol from all of the simulations. Modifying V2 does not affect the resulting cis/trans equilibrium that is partially governed by nonbonded and electrostatic interactions. In addition, the 5-6% of cis Xaa-Pro peptide bonds reported in surveys of three-dimensional structures of proteins13–15 corresponds to a relative stability of ∼1.6-1.8 kcal/mol that is not in accord with estimates from equilibrium experiments on Xaa-Pro dipeptides. We next investigated the cis-trans isomerization of the ω bond (Figure 1) in NMA by an accelerated MD approach using the reoptimized V2. This time, we recalibrated the dihedral constants V1 and V2′ for the H-N-C-O torsion (Figure 1). V1 and V2′ for the H-N-C-O torsion are used in AMBER force field to maintain the correct free energy difference between the cis and the trans isomers for non-proline peptide ω bonds. A comparison of the original and revised AMBER parameters is shown in Table 3. Accelerated MD simulations with default AMBER parameters V1 ) 4.0, V2′ ) 5.0, and V2 ) 20.0 kcal/ mol resulted in ∆G#tfc ) 17.1 ( 0.1 kcal/mol and a free energy difference of 4.6 ( 0.3 kcal/mol between the cis and trans isomers (∆Gc-t). Thus, neither the experimental isomerization
TABLE 3: Comparison of Original and Revised AMBER Parameters for ω Torsional Potential Vn/2 (kcal/mol)
torsion X-C-N-X H-N-C-O H-N-C-O
phase no. of bond optimized periodicity, angle, γ (degrees) paths original (this work) n 4 1 1
10 2.5 2.0
14 2.9 1.0
2 2′ 1
180 180 0
barrier nor the cis/trans equilibrium was reproduced for NMA. Our modified set of parameters V1 ) 2.0, V2′ ) 5.8, and V2 ) 28.0 kcal/mol yielded the free energy profile shown in Figure 2B, and ∆G#tfc, ∆G#cft, and ∆Gc-t obtained from these are listed in Table 2. Our results are within 0.2 kcal/mol of the respective experimental values, which is less than the average error from simulations as well as experimental error of 0.3 kcal/mol. Testing the Transferability of the Modified ω Torsional Parameters on Xaa-Pro and Xaa-non-Pro Dipeptides. Cis-trans isomerization of Ala-Pro and Phe-Pro has been probed by various methods, including proton NMR, capillary zone electrophoresis, and pH jump studies.27,43–46 Table 2 also summarizes the free energy barrier heights for cis-trans isomerization and free energy differences between the cis and trans states for these peptides estimated from experiments as well as from accelerated MD simulations using the modified dihedral constant; that is, V2 ) 28.0 kcal/mol. In the case of Phe-Pro, the cis isomer is found to be more stable than the trans state from both experiments and simulations. For Ala-Pro, the simulations calculated the stability of the cis isomer to be slightly more than that of trans, in contrast to experimental findings. However, this observation does not mean much, since the relative free energy difference between the cis and the trans basins calculated from experiment is very small (
