AMBER Parm99 Research Articles

Optimization of an AMBER force field for the artificial nucleic acid, LNA, and benchmarking with NMR of L(CAAU).

Locked Nucleic Acids (LNAs) are RNA analogues with an O2′-C4′ methylene bridge which locks the sugar into a C3′-endo conformation. This enhances hybridization to DNA and RNA, making LNAs useful in microarrays and potential therapeutics. Here, the LNA, L(CAAU), provides a simplified benchmark for testing the ability of molecular dynamics (MD) to approximate nucleic acid properties. LNA χ torsions and partial charges were parametrized to create AMBER parm99_LNA. The revisions were tested by comparing MD predictions with AMBER parm99 and parm99_LNA against a 200 ms NOESY NMR spectrum of L(CAAU). NMR indicates an A-Form equilibrium ensemble. In 3000 ns simulations starting with an A-form structure, parm99_LNA and parm99 provide 66% and 35% agreement, respectively, with NMR NOE volumes and 3J-couplings. In simulations of L(CAAU) starting with all χ torsions in a syn conformation, only parm99_LNA is able to repair the structure. This implies methods for parametrizing force fields for nucleic acid mimics can reasonably approximate key interactions and that parm99_LNA will improve reliability of MD studies for systems with LNA. A method for approximating χ population distribution on the basis of base to sugar NOEs is also introduced.

The nuclear magnetic resonance of CCCC RNA reveals a right-handed helix, and revised parameters for AMBER force field torsions improve structural predictions from molecular dynamics.

The sequence dependence of RNA energetics is important for predicting RNA structure. Hairpins with Cn loops are consistently less stable than hairpins with other loops, which suggests the structure of Cn regions could be unusual in the “unfolded” state. For example, previous nuclear magnetic resonance (NMR) evidence suggested that polycytidylic acid forms a left-handed helix. In this study, UV melting experiments show that the hairpin formed by r(5′GGACCCCCGUCC) is less stable than r(5′GGACUUUUGUCC). NMR spectra for single-stranded C4 oligonucleotide, mimicking the unfolded hairpin loop, are consistent with a right-handed A-form-like helix. Comparisons between NMR spectra and molecular dynamics (MD) simulations suggest that recent reparametrizations, parm99χ_YIL and parm99TOR, of the AMBER parm99 force field improve the agreement between structural features for C4 determined by NMR and predicted by MD. Evidently, the force field revisions to parm99 improve the modeling of RNA energetics and therefore structure.

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica-exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free-energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica-exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local-minimum free-energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid-β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β-helix, α-helix, 3(10)-helix, β-hairpin, and β-sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free-energy landscape.

Folding simulations of three proteins having all α-helix, all β-strand and α/β-structures

We performed folding simulations of three proteins using four force fields, AMBER parm96, AMBER parm99, CHARMM 27 and OPLS-AA/L, in order to examine the features of these force fields. We studied three proteins, protein A (all α-helix), cold-shock protein (all β-strand) and protein G (α/β-structures), for the folding simulations. For the simulation, we used the simulated annealing molecular dynamics method, which was performed 50 times for each protein using the four force fields. The results showed that the secondary-structure-forming tendencies are largely different among the four force fields. AMBER parm96 favours β-bridge structures and extended β-strand structures, and AMBER parm99 favours α-helix structures and 310-helix structures. CHARMM 27 slightly favours α-helix structures, and there are also π-helix and β-bridge structures. OPLS-AA/L favours α-helix structures and 310-helix structures.

Controlling the secondary-structure-forming tendencies of proteins by a backbone torsion-energy term

We examined a new backbone torsion-energy term proposed by us in the force field for protein systems. This torsion-energy term is represented by a double Fourier series in two variables, namely the backbone dihedral angles φ and ψ. It gives a natural representation of the torsion energy in the Ramachandran space in the sense that any two-dimensional energy surface periodic in both φ and ψ can be expanded by the double Fourier series. We can then easily control secondary-structure-forming tendencies by modifying the torsion-energy surface. For instance, we can increase or decrease the α-helix-forming-tendencies by lowering or raising the torsion-energy surface in the α-helix region and likewise increase or decrease the β-sheet-forming tendencies by lowering or raising the surface in the β-sheet region in the Ramachandran space. We applied this torsion-energy modification method to six force fields, AMBER parm94, AMBER parm96, AMBER parm99, CHARMM27, OPLS-AA and OPLS-AA/L, and demonstrated that our modifications of the torsion-energy terms resulted in the expected changes of secondary-structure-forming tendencies by performing folding simulations of α-helical and β-hairpin peptides.

Determination method of the balance of the secondary-structure-forming tendencies of force fields

We propose a new method of optimisation of backbone torsion-energy parameters in the force field for molecular simulations of protein systems. This method is based on the idea of balancing the secondary-structure-forming tendencies, namely, those of α-helix and β-sheet structures. We perform a minimisation of the backbone dihedral angle-based root-mean-square deviation of the helix and β structure regions in many protein structures. As an example, we optimised the backbone torsion-energy parameters of AMBER parm96 force field using 100 protein molecules from the Protein Data Bank. We then performed folding simulations of α-helical and β-hairpin peptides, using the optimised force field. The results imply that the new force-field parameters give structures more consistent with the experimental implications than the original AMBER parm96 force field.

Quantum mechanical studies on model α‐pleated sheets

Pauling and Corey proposed a pleated-sheet configuration, now called alpha-sheet, as one of the protein secondary structures in addition to alpha-helix and beta-sheet. Recently, it has been suggested that alpha-sheet is a common feature of amyloidogenic intermediates. We have investigated the stability of antiparallel beta-sheet and two conformations of alpha-sheet in solution phase using the density functional theoretical method. The peptides are modeled as two-strand acetyl-(Ala)(2)-N-methylamine. Using stages of geometry optimization and single point energy calculation at B3LYP/cc-pVTZ//B3LYP/6-31G* level and including zero-point energies, thermal, and entropic contribution, we have found that beta-sheet is the most stable conformation, while the alpha-sheet proposed by Pauling and Corey has 13.6 kcal/mol higher free energy than the beta-sheet. The alpha-sheet that resembles the structure observed in molecular dynamics simulations of amyloidogenic proteins at low pH becomes distorted after stages of geometry optimization in solution. Whether the alpha-sheets with longer chains would be increasingly favorable in water relative to the increase in internal energy of the chain needs further investigation. Different from the quantum mechanics results, AMBER parm94 force field gives small difference in solution phase energy between alpha-sheet and beta-sheet. The predicted amide I IR spectra of alpha-sheet shows the main band at higher frequency than beta-sheet.

Polarization and Polarizability Assessed by Protein Amide Acidity

Hydroxide-catalyzed exchange rate constants were determined for those amides of FK506-binding protein (FKBP12), ubiquitin, and chymotrypsin inhibitor 2 (CI2) that are solvent-accessible in the high-resolution X-ray structures. When combined with previous hydrogen exchange results for the rubredoxin from Pyrococcus furiosus, the acidity of these amides was calculated by continuum dielectric methods as a function of the nonpolarizable electrostatic parameter set, internal dielectric, and the charge distribution of the peptide anion. The CHARMM22 parameter set with an internal dielectric value of 3 and an ab initio-derived anion charge distribution yielded an rmsd value of 7 for the 56 amide exchange rate constants ranging from 10(0.67) to 10(9.0) M(-1) s(-1). The OPLS-AA parameter set yielded comparably robust predictions, while that of PARSE, AMBER parm99, and AMBER ff03 performed more poorly. The small value for the optimal internal dielectric, combined with the brief lifetime of the peptide anion intermediate and the uniformity of the correlation between predicted and observed amide acidities, is consistent with electronic polarizability providing the dominant contribution to dielectric shielding. By construction, nonpolarizable force fields do not model electric field attenuation by electronic polarizability. Accurate prediction of the total electrostatic energy by such force fields necessitates the hyperpolarization of the atomic charge values in order to match the average electric field energy density (1/2)epsilon(tau)E(2)(tau) when epsilon(tau) is set to the in vacuo dielectric value of 1. The resulting predictions of the experimental hydrogen exchange data demonstrate the substantial systematic errors in the predicted electrostatic potential that can arise when dielectric shielding due to electronic polarizability is neglected.

Temperature and Pressure Dependence of Alanine Dipeptide Studied by Multibaric−Multithermal Molecular Dynamics Simulations

We applied the multibaric-multithermal (MUBATH) molecular dynamics (MD) algorithm to an alanine dipeptide in explicit water. The MUBATH MD simulation covered a wide range of conformational space and sampled the states of PII, C5, alphaR, alphaP, alphaL, and C7(ax). On the other hand, the conventional isobaric-isothermal simulation was trapped in local-minimum free-energy states and sampled only a few of them. We calculated the partial molar enthalpy difference DeltaH and partial molar volume difference DeltaV among these states by the MUBATH simulation using the AMBER parm99 and AMBER parm96 force fields and two sets of initial conditions. We compared these results with those from Raman spectroscopy experiments. The Raman spectroscopy data of DeltaH for the C5 state against the PII state agreed with both MUBATH data with the AMBER parm96 and parm99 force fields. The partial molar enthalpy difference DeltaH for the alphaR state and the partial molar volume difference DeltaV for the C5 state by the Raman spectroscopy agreed with those for the AMBER parm96 force field. On the other hand, DeltaV for the alphaR state by the Raman spectroscopy was consistent with our AMBER-parm99 force-field result. All the experimental results fall between those of simulations using AMBER parm96 and parm99 force fields, suggesting that the ideal force-field parameters lie between those of AMBER parm96 and parm99.

A Billion-fold Range in Acidity for the Solvent-Exposed Amides of Pyrococcus furiosus Rubredoxin

The exchange rates of the static solvent-accessible amide hydrogens of Pyrococcus furiosus rubredoxin range from near the diffusion-limited rate to a billion-fold slower for the non-hydrogen-bonded Val 38 (eubacterial numbering). Hydrogen exchange directly monitors the kinetic acidity of the peptide nitrogen. Electrostatic solvation free energies were calculated by Poisson-Boltzmann methods for the individual peptide anions that form during the hydroxide-catalyzed exchange reaction to examine how well the predicted thermodynamic acidities match the experimentally determined kinetic acidities. With the exception of the Ile 12 amide, the differential exchange rate constant for each solvent-exposed amide proton that is not hydrogen bonded to a backbone carbonyl can be predicted within a factor of 6 (10 (0.78)) root-mean-square deviation (rmsd) using the CHARMM22 electrostatic parameter set and an internal dielectric value of 3. Under equivalent conditions, the PARSE parameter set yields a larger rmsd value of 1.28 pH units, while the AMBER parm99 parameter set resulted in a considerably poorer correlation. Either increasing the internal dielectric value to 4 or reducing it to a value of 2 significantly degrades the quality of the prediction. Assigning the excess charge of the peptide anion equally between the peptide nitrogen and the carbonyl oxygen also reduces the correlation to the experimental data. These continuum electrostatic calculations were further analyzed to characterize the specific structural elements that appear to be responsible for the wide range of peptide acidities observed for these solvent-exposed amides. The striking heterogeneity in the potential at sites along the protein-solvent interface should prove germane to the ongoing challenge of quantifying the contribution that electrostatic interactions make to the catalytic acceleration achieved by enzymes.

Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/ γ Conformers

We present here the parmbsc0 force field, a refinement of the AMBER parm99 force field, where emphasis has been made on the correct representation of the α/ γ concerted rotation in nucleic acids (NAs). The modified force field corrects overpopulations of the α/ γ = (g+,t) backbone that were seen in long (more than 10 ns) simulations with previous AMBER parameter sets (parm94-99). The force field has been derived by fitting to high-level quantum mechanical data and verified by comparison with very high-level quantum mechanical calculations and by a very extensive comparison between simulations and experimental data. The set of validation simulations includes two of the longest trajectories published to date for the DNA duplex (200 ns each) and the largest variety of NA structures studied to date (15 different NA families and 97 individual structures). The total simulation time used to validate the force field includes near 1 μs of state-of-the-art molecular dynamics simulations in aqueous solution.

Molecular dynamics study of the conformational dynamics and energetics of some large‐ring cyclodextrins (CDn, n = 24, 25, 26, 27, 28, 29)

Molecular dynamics simulations in water solution were performed on six large-ring cyclodextrins (LR-CDs) with a degree of polymerization 24, 25, 26, 27, 28, and 29. The AMBER parm99 force field and explicit water molecules (TIP3P) were used in the simulations. The present research was aimed at further extending our knowledge on the structural dynamics and the energetics of this new class of compounds that may eventually provide chiral cavities suitable for formation of inclusion complexes with small molecules, and, accordingly, to serve as host structures for chiral recognition. The study focused on several representatives flanking CD26-the largest LR-CD for which X-ray data is available. Both the monitoring of the structural variations during the simulations as well as the analyses of energy balances are indicative for high flexibility of the macrorings. Slight differences of the overall preferred shapes were detected with diminishing the size of the macromolecules from CD29 to CD24. An elongated cavity (CD28) or a double parallel strand in different specific representations are the dominating motifs in the LR-CDs studied: with loops at the two ends (CD25, CD28, CD29), with a loop at one end (CD25), twisted (CD26, CD27) or twisted with an open portion in the middle (CD24), helical (CD24, CD25), or linking two loops from one of their sides (CD27). Two loops connected by an arc (CD28, CD29) and a cavity with the shape of an extended rectangular (CD24, CD28) appear preferentially during the conformational interconversions of the two larger CDs, whereas helical motifs are present in the smaller macrorings: an extended helix with ends linked by an arc (CD24), helical turn and helical portion (CD26, CD27). A triple propeller conformation or three symmetrical loops of almost equal size were also detected for CD26 and CD29, respectively. The present results further support the hypothesis for the existence of more than one cavity in large-ring cyclodextrins and suggest preferred conformations in water solution for the LR-CDs with degree of polymerization from 24 to 29.

Large‐ring cyclodextrins. Further support for the preferred conformations of CD26 in water solution: Molecular dynamics studies on CD26‐derived conformations of CDn (n = 24, 25, 27, 28, 29, 30)

AbstractMolecular dynamics (MD) simulations in water solution were performed on six large ring cyclodextrins (LR‐CDs) with a degree of polymerization (DP) 24, 25, 27, 28, 29, and 30, using as starting structures geometries derived from CD26, which was studied previously by MD. The present research was aimed at further extending our knowledge on the structural dynamics and the energetics of this new class of compounds, focusing on several representatives flanking CD26, the largest LR‐CD for which X‐ray data are available. The most recent parameterization for carbohydrates, Glycam‐04, and explicit water molecules (TIP3P) were used in the 10.0‐ns simulations. Preceding 10.0‐ns simulations of CD26, with the AMBER parm99 force field, revealed that the overall shape of the macromolecule did not change significantly during the last 5.0‐ns simulations. We now further test for the probability of prevalence of this conformation by using a different force field and by examining the effect of the size of the macro‐ring on the preferred conformations of LR‐CDs with DP within a range of 24–30. Indeed, we found CD26‐like geometries present for CD25, CD27, and CD28, which could be interpreted as an indication for an enhanced stability of the CD26 conformation in water solution. Again, we witnessed high flexibility of the macro‐cycles, and different local structural motifs were detected: twisted helical double parallel strand (CD25, CD30), single helix of three turns (CD29), a spiral (CD25, CD27, CD30), open‐shaped loops and helical turns (CD28; “CD26‐like geometry”), a rectangular cavity (CD24), a triple‐lobe fragment (CD29) or a portion shaped like a jaw (CD24, CD27), and an overall shape of the macromolecule that resembled a butterfly (CD27) or “X” (CD28). Thus, we conclude that the same conformation of CD26 could be also the preferred conformation in water solution of other LR‐CDs with DP in the close vicinity of 26. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2007

Secondary-Structure Design of Proteins by a Backbone Torsion Energy

We propose a new backbone-torsion-energy term in the force field for protein systems. This torsion-energy term is represented by a double Fourier series in two variables, the backbone dihedral angles phi and psi. It gives a natural representation of the torsion energy in the Ramachandran space in the sense that any two-dimensional energy surface periodic in both phi and psi can be expanded by the double Fourier series. We can then easily control secondary-structure-forming tendencies by modifying the torsion-energy surface. For instance, we can increase/decrease the alpha-helix-forming-tendencies by lowering/raising the torsion-energy surface in the alpha-helix region and likewise increase/decrease the beta-sheet-forming tendencies by lowering/raising the surface in the beta-sheet region in the Ramachandran space. We applied our approach to AMBER parm94 and AMBER parm96 force fields and demonstrated that our modifications of the torsion-energy terms resulted in the expected changes of secondary-structure-forming-tendencies by performing folding simulations of alpha-helical and beta-hairpin peptides.

Theoretical Study of the Unusual Protonation Properties of the Active Site Cysteines in Thioredoxin

The thiol/disulfide oxidoreductases of the thioredoxin family have, in the active site, two cysteines that can be in a reduced or an oxidized form. One of the cysteines in the reduced state is deprotonated, and it is called nucleophilic cysteine. The pK(a) of this cysteine is different from that of a normal cysteine and varies widely among the different enzymes of this family. However, the factors responsible for the different degrees of stabilization of nucleophilic cysteine thiolate are not fully understood. Here, we have studied the well-known hypothesis of proton sharing between the active site thiols by performing a linear transit scan for the transfer of the proton between the active site cysteines. We used a two-layered (DFT/MM) ONIOM formalism, with the active site region treated at the B3LYP/6-31+G(d) level and the remains of the protein treated with the Amber Parm94 force field. The solvation free energy was accounted for with a continuum solvent model, by solving the Poisson-Boltzmann equation using the program Delphi. We have obtained excellent agreement with the experimental data available in the literature. Besides refuting the proton sharing hypothesis, our results include a value of 14.0 for the pK(a) of the buried cysteine, a quantity that has not been possible to obtain experimentally but which has been proven to be higher than 11. Additionally, this study also provides detailed information on the very interesting and so far unknown fact that the contribution of the enzymatic structure (8.3 kcal/mol) prevails in relation to that of the solvent (0.60 kcal/mol) concerning the differential stabilization of the active site thiolates.

AMBER-based hybrid force field for conformational sampling of polypeptides

AMBER-based hybrid force field was developed, with mixing the mainchain torsion energies of AMBER parm94 and parm96 force fields, with introducing a weighting factor to control the mixing ratio. First, dependence of an adiabatic Φ– Ψ potential–energy map of alanine dipeptide on the weighting factor was investigated, and the map was compared with an ab initio Φ– Ψ map. Second, we did enhanced conformational sampling of two six-residue peptides known to form a helical or turn conformation in solution. The hybrid force field reproduced the ab initio energy better than either parm94 or parm96, and exhibited the secondary-structure preference depending on the amino acid sequence.

PROTEIN FORCE-FIELD PARAMETERS OPTIMIZED WITH THE PROTEIN DATA BANK I: FORCE-FIELD OPTIMIZATIONS

We optimized five existing sets of force-field parameters for protein systems by our recently proposed method. The five force fields are AMBER parm94, AMBER parm96, AMBER parm99, CHARMM version 22, and OPLS-AA. The method consists of minimizing the sum of the square of the force acting on each atom in the proteins with the structures from the Protein Data Bank (PDB). We selected the partial-charge and backbone torsion-energy parameters for this optimization, and 100 molecules from the PDB were used. We gave detailed comparisons of the optimized force fields and found that there is a tendency of convergence towards the same function for the torsion-energy term.

PROTEIN FORCE-FIELD PARAMETERS OPTIMIZED WITH THE PROTEIN DATA BANK II: COMPARISONS OF FORCE FIELDS BY FOLDING SIMULATIONS OF SHORT PEPTIDES

In Paper I of this series, the formulations of the optimization method of existing force-field parameters for protein systems have been presented. We then applied it to five sets of force-field parameters, namely, AMBER parm94, AMBER parm96, AMBER parm99, CHARMM version 22, and OPLS-AA. In order to test the validity of these force fields, the folding simulations of α-helical and β-hairpin peptides have been performed with each of the original and optimized force-field parameters. We found that all five modified force-field parameters gave both α-helical and β-hairpin structures more consistent with the experimental implications than the original force fields.

Insights into the Structure of Large-Ring Cyclodextrins through Molecular Dynamics Simulations in Solution

Molecular dynamics simulations were carried out using the AMBER parm99 force field and explicit water molecules (TIP3P) to gain insight into the structural deformations and energetics of several large-ring cyclodextrins (with a degree of polymerization 26, 30, 55, 70, 85, and 100) in solution. The structures displayed by CD26 during the MD simulation (10.0 ns) did not correspond to the conformation in the crystalline state. The two "flips" present in the macroring of CD26 in the crystal state disappeared after 1.5 and 3.0 ns simulations, respectively. The larger CDs bear a considerable degree of flexibility. They display different modes of folding and cavity-like regions of different sizes and shapes: circular and elongated loops of variable size, orientated in different fashion; portions of a double helical strand with the two single helices parallel to each other (CD30, CD70); a helix of three turns and a serpentiform portion containing six loops (CD55); a cone-shaped spiral region (CD70); a rounded dendritic fold with several arbitrarily oriented small loops on the surface of the clustering (CD85); three spiral portions and a tendency for bending into two (CD100). These results support the hypothesis for the existence of more than one cavity in large-ring cyclodextrins.

Optimization of protein force-field parameters with the Protein Data Bank

We propose a novel method to optimize existing force-field parameters for protein systems. The method consists of minimizing the summation of the square of the force acting on each atom in the proteins with the structures from the Protein Data Bank. We performed this optimization to the partial-charge and torsion-energy parameters of the AMBER parm94 force field, using 100 molecules from the Protein Data Bank. We then performed folding simulations of α-helical and β-hairpin peptides. The optimized force-field parameters gave structures more consistent with the experimental implications than the original AMBER force field.