Comparison of Secondary Structure Formation Using 10 Different Force Fields in Microsecond Molecular Dynamics Simulations

We have compared molecular dynamics (MD) simulations of a β-hairpin forming peptide derived from the protein Nrf2 with 10 biomolecular force fields using trajectories of at least 1 μs. The total simulation time was 37.2 μs. Previous studies have shown that different force fields, water models, simulation methods, and parameters can affect simulation outcomes. The MD simulations were done in explicit solvent with a 16-mer Nrf2 β-hairpin forming peptide using Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, GROMOS96 53a6, CHARMM27, and OPLS-AA/L force fields. The effects of charge-groups, terminal capping, and phosphorylation on the peptide folding were also examined. Despite using identical starting structures and simulation parameters, we observed clear differences among the various force fields and even between replicates using the same force field. Our simulations show that the uncapped peptide folds into a native-like β-hairpin structure at 310 K when Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, or GROMOS96 53a6 were used. The CHARMM27 simulations were able to form native hairpins in some of the elevated temperature simulations, while the OPLS-AA/L simulations did not yield native hairpin structures at any temperatures tested. Simulations that used charge-groups or peptide capping groups were not largely different from their uncapped counterparts with single atom charge-groups. On the other hand, phosphorylation of the threonine residue located at the β-turn significantly affected the hairpin formation. To our knowledge, this is the first study comparing such a large set of force fields with respect to β-hairpin folding. Such a comprehensive comparison will offer useful guidance to others conducting similar types of simulations.     

Performances of Five Representative Force Fields on Gaseous Amino Acids with Di erent Termini

There is a growing interest in the study of structures and properties of biomolecules in gas phase. Applications of force fields are highly desirable for the computational efficiency of the gas phase study. To help the selection of force fields, the performances of five repre-sentative force fields for gaseous neutral, protonated, deprotonated and capped amino acids are systematically examined and compared. The tested properties include relative conforma-tional energies, energy differences between cis and trans structures, the number and strength of predicted hydrogen bonds, and the quality of the optimized structures. The results of BHandHLYP/6-311++G(d,p) are used as the references. GROMOS53A6 and ENCADS are found to perform poorly for gaseous biomolecules, while the performance of AMBER99SB, CHARMM27 and OPLSAA/L are comparable when applicable. Considering the general availability of the force field parameters, CHARMM27 is the most recommended, followed by OPLSAA/L, for the study of biomolecules in gas phase.     

Parameterization of Ca+2-protein interactions for molecular dynamics simulations

Molecular dynamics simulations of Ca+2 ions near protein were performed with three force fields: GROMOS96, OPLS-AA, and CHARMM22. The simulations reveal major, force-field dependent, inconsistencies in the interaction between the Ca+2 ions with the protein. The variations are attributed to the nonbonded parameterizations of the Ca+2-carboxylates interactions. The simulations results were compared to experimental data, using the Ca+2-HCOO- equilibrium as a model. The OPLS-AA force field grossly overestimates the binding affinity of the Ca+2 ions to the carboxylate whereas the GROMOS96 and CHARMM22 force fields underestimate the stability of the complex. Optimization of the Lennard-Jones parameters for the Ca+2-carboxylate interactions were carried out, yielding new parameters which reproduce experimental data.     

Evaluating nonpolarizable nucleic acid force fields: A systematic comparison of the nucleobases hydration free energies and chloroform‐to‐water partition coefficients

Nucleic acid force fields have been shown to reproduce structural properties of DNA and RNA very well, but comparative studies with respect to thermodynamic properties are rare. As a test for thermodynamic properties, we have computed hydration free energies and chloroform‐to‐water partition coefficients of nucleobases using the AMBER‐99, AMBER‐gaff, CHARMM‐27, GROMOS‐45a4/53a6 and OPLS‐AA force fields. A mutual force field comparison showed a very large spread in the calculated thermodynamic properties, demonstrating that some of the parameter sets require further optimization. The choice of solvent model used in the simulation does not have a significant effect on the results. Comparing the hydration free energies obtained by the various force fields to the adenine and thymine experimental values showed a very large deviation for the GROMOS and AMBER parameter sets. Validation against experimental partition coefficients showed good agreement for the CHARMM‐27 parameter set. In view of mutation studies, differences in partition coefficient between two bases were also compared, and good agreement between experiments and calculations was found for the AMBER‐99 parameter set. Overall, the CHARMM‐27 parameter set performs best with respect to the thermodynamic properties tested here. © 2012 Wiley Periodicals, Inc.     

Modelling of carbohydrate–aromatic interactions: <Emphasis Type="Italic">ab initio </Emphasis>energetics and force field performance

Summary Aromatic amino acid residues are often present in carbohydrate-binding sites of proteins. These binding sites are characterized by a placement of a carbohydrate moiety in a stacking orientation to an aromatic ring. This arrangement is an example of CH/π interactions. Ab initio interaction energies for 20 carbohydrate–aromatic complexes taken from 6 selected ultra-high resolution X-ray structures of glycosidases and carbohydrate-binding proteins were calculated. All interaction energies of a pyranose moiety with a side chain of an aromatic residue were calculated as attractive with interaction energy ranging from −2.8 to −12.3 kcal/mol as calculated at the MP2/6-311+G(d) level. Strong attractive interactions were observed for a wide range of orientations of carbohydrate and aromatic ring as present in selected X-ray structures. The most attractive interaction was associated with apparent combination of CH/π interactions and classical H-bonds. The failure of Hartree–Fock method (interaction energies from +1.0 to −6.9 kcal/mol) can be explained by a dispersion nature of a majority of the studied complexes. We also present a comparison of interaction energies calculated at the MP2 level with those calculated using molecular mechanics force fields (OPLS, GROMOS, CSFF/CHARMM, CHEAT/CHARMM, Glycam/AMBER, MM2 and MM3). For a majority of force fields there was a strong correlation with MP2 values. RMSD between MP2 and force field values were 1.0 for CSFF/CHARMM, 1.2 for Glycam/AMBER, 1.2 for GROMOS, 1.3 for MM3, 1.4 for MM2, 1.5 for OPLS and to 2.3 for CHEAT/CHARMM (in kcal/mol). These results show that molecular mechanics approximates interaction energies very well and support an application of molecular mechanics methods in the area of glycochemistry and glycobiology.     

Modern protein force fields behave comparably in molecular dynamics simulations

Several molecular dynamics simulations were performed on three proteins--bovine apo-calbindin D9K, human interleukin-4 R88Q mutant, and domain IIA of bacillus subtilis glucose permease--with each of the AMBER94, CHARMM22, and OPLS-AA force fields as implemented in CHARMM. Structural and dynamic properties such as solvent-accessible surface area, radius of gyration, deviation from their respective experimental structures, secondary structure, and backbone order parameters are obtained from each of the 2-ns simulations for the purpose of comparing the protein portions of these force fields. For one of the proteins, the interleukin-4 mutant, two independent simulations were performed using the CHARMM22 force field to gauge the sensitivity of some of these properties to the specific trajectory. In general, the force fields tested performed remarkably similarly with differences on the order of those found for the two independent trajectories of interleukin-4 with CHARMM22. When all three proteins are considered together, no force field showed any consistent trend in variations for most of the properties monitored in the study.     

Are Protein Force Fields Getting Better? A Systematic Benchmark on 524 Diverse NMR Measurements

Recent hardware and software advances have enabled simulation studies of protein systems on biophysically-relevant timescales, often revealing the need for improved force fields. Although early force field development was limited by the lack of direct comparisons between simulation and experiment, recent work from several labs has demonstrated direct calculation of NMR observables from protein simulations. Here we quantitatively evaluate recent molecular dynamics force fields against a suite of 524 chemical shift and J coupling ((3)JH(N)H(α), (3)JH(N)C(β), (3)JH(α)C', (3)JH(N)C', and (3)JH(α)N) measurements on dipeptides, tripeptides, tetra-alanine, and ubiquitin. Of the force fields examined (ff96, ff99, ff03, ff03*, ff03w, ff99sb*, ff99sb-ildn, ff99sb-ildn-phi, ff99sb-ildn-nmr, CHARMM27, OPLS-AA), two force fields (ff99sb-ildn-phi, ff99sb-ildn-nmr) combining recent side chain and backbone torsion modifications achieve high accuracy in our benchmark. For the two optimal force fields, the calculation error is comparable to the uncertainty in the experimental comparison. This observation suggests that extracting additional force field improvements from NMR data may require increased accuracy in J coupling and chemical shift prediction. To further investigate the limitations of current force fields, we also consider conformational populations of dipeptides, which were recently estimated using vibrational spectroscopy.     

Electrostatic solvation free energy of amino acid side chain analogs: implications for the validity of electrostatic linear response in water

Electrostatic free energies of solvation for 15 neutral amino acid side chain analogs are computed. We compare three methods of varying computational complexity and accuracy for three force fields: free energy simulations, Poisson-Boltzmann (PB), and linear response approximation (LRA) using AMBER, CHARMM, and OPLS-AA force fields. We find that deviations from simulation start at low charges for solutes. The approximate PB and LRA produce an overestimation of electrostatic solvation free energies for most of molecules studied here. These deviations are remarkably systematic. The variations among force fields are almost as large as the variations found among methods. Our study confirms that success of the approximate methods for electrostatic solvation free energies comes from their ability to evaluate free energy differences accurately.     

Effects of all-atom force fields on amyloid oligomerization: replica exchange molecular dynamics simulations of the Aβ(16-22) dimer and trimer

The aim of this work is to investigate the effects of molecular mechanics force fields on amyloid peptide assembly. To this end, we performed extensive replica exchange molecular dynamics (REMD) simulations on the monomer, dimer and trimer of the seven-residue fragment of the Alzheimer's amyloid-β peptide, Aβ(16-22), using the AMBER99, GROMOS96 and OPLS force fields. We compared the force fields by analysing the resulting global and local structures as well as the free energy landscapes at 300 K. We show that AMBER99 strongly favors helical structures for the monomer and does not predict any β-sheet structure for the dimer and trimer. In contrast, the dimer and trimer modeled by GROMOS96 form antiparallel β-sheet structures, while OPLS predicts diverse structures. Overall, the free energy landscapes obtained by three force fields are very different, and we also note a weak structural dependence of our results on temperature. The implications of this computational study on amyloid oligomerization, fibril growth and inhibition are also discussed.     

Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models

We present an extensive study on hydration thermodynamic properties of analogues of 13 amino acid side chains at 298 K and 1 atm. The hydration free energies DeltaG, entropies DeltaS, enthalpies DeltaH, and heat capacities Deltac(P)() were determined for 10 combinations of force fields and water models. The statistical sampling was extended such that precisions of 0.3, 0.8, 0.8 kJ/mol and 25 J/(mol K) were reached for DeltaG, TDeltaS, DeltaH, and Deltac(P)(), respectively. The three force fields used in this study are AMBER99, GROMOS 53A6, and OPLS-AA; the five water models are SPC, SPC/E, TIP3P, TIP4P, and TIP4P-Ew. We found that the choice of water model strongly influences the accuracy of the calculated hydration entropies, enthalpies, and heat capacities, while differences in accuracy between the force fields are small. On the basis of an analysis of the hydrophobic analogues of the amino acid side chains, we discuss what properties of the water models are responsible for the observed discrepancies between computed and experimental values. The SPC/E water model performs best with all three biomolecular force fields.     

Strength of hydrogen bonds in α helices

The intramolecular helix backbone >C(DOUBLE BOND)O---H—N< hydrogen bonding energies in poly(L -alanine) α helices have been estimated both in vacuum and in an aqueous environment using the parameter sets of five of the most commonly used force fields for modeling biomolecules, namely AMBER, CHARMM, ECEPP, GROMOS, and OPLS. The relative capabilities of these force fields in describing the H-bonding interactions with different dielectric continuum models have been assessed. A modified Hingerty–Lavery function is proposed for the treatment of electrostatic interactions of biomolecules in an aqueous environment. The helix backbone H-bonding energies predicted by this function (∼−1 kcal/mol) correspond closely with the experiment. © 1997 John Wiley & Sons, Inc. J Comput Chem 18: 1245–1252     

Free energies of hydration from thermodynamic integration: Comparison of molecular mechanics force fields and evaluation of calculation accuracy

Four commonly used molecular mechanics force fields, CHARMM22, OPLS, CVFF, and GROMOS87, are compared for their ability to reproduce experimental free energies of hydration (ΔG_hydr) from molecular dynamics (MD) simulations for a set of small nonpolar and polar organic molecules: propane, cyclopropane, dimethylether, and acetone. ΔG_hydr values were calculated by multiconfiguration thermodynamic integration for each of the different force fields with three different sets of partial atomic charges: full charges from an electrostatic potential fit (ESP), and ESP charges scaled by 0.8 and 0.6. All force fields, except for GROMOS87, give reasonable results for ΔG_hydr · if partial atomic charges of appropriate magnitude are assigned. For GROMOS87, the agreement with experiment for hydrocarbons (propane and ethane) was improved considerably by modifying the repulsive part of the carbon-water oxygen Lennard-Jones potential. The small molecules studied are related to the chemical moieties constituting camphor (C₁₀H_l6O). By invoking force-field transferability, we calculated the ΔG_hydr for camphor. With the modified GROMOS force field, a ΔG_hydr within 4 kJ/mol of the experimental value of −14.8 kJ/mol was obtained. Camphor is one of the largest molecules for which an absolute hydration free energy has been calculated by molecular simulation. The accuracy and reliability of the thermodynamic integration calculations were analyzed in detail and we found that, for ΔG_hydr calculations for the set of small molecules in aqueous solution, molecular dynamics simulations of 0.8–1.0 ns in length give an upper statistical error bound of 1.5 kJ/mol, whereas shorter simulations of 0.25 nm in length given an upper statistical error bound of 3.5 kJ/mol. © 1997 by John Wiley & Sons, Inc.     

Revision of AMBER Torsional Parameters for RNA Improves Free Energy Predictions for Tetramer Duplexes with GC and iGiC Base Pairs

All-atom force fields are important for predicting thermodynamic, structural, and dynamic properties of RNA. In this paper, results are reported for thermodynamic integration calculations of free energy differences of duplex formation when CG pairs in the RNA duplexes r(CCGG)(2), r(GGCC)(2), r(GCGC)(2), and r(CGCG)(2) are replaced by isocytidine-isoguanosine (iCiG) pairs. Agreement with experiment was improved when ε/ζ, α/γ, β, and χ torsional parameters in the AMBER99 force field were revised on the basis of quantum mechanical calculations. The revised force field, AMBER99TOR, brings free energy difference predictions to within 1.3, 1.4, 2.3, and 2.6 kcal/mol at 300 K, respectively, compared to experimental results for the thermodynamic cycles of CCGG → iCiCiGiG, GGCC → iGiGiCiC, GCGC → iGiCiGiC, and CGCG → iCiGiCiG. In contrast, unmodified AMBER99 predictions for GGCC → iGiGiCiC and GCGC → iGiCiGiC differ from experiment by 11.7 and 12.6 kcal/mol, respectively. In order to test the dynamic stability of the above duplexes with AMBER99TOR, four individual 50 ns molecular dynamics (MD) simulations in explicit solvent were run. All except r(CCGG)(2) retained A-form conformation for ≥82% of the time. This is consistent with NMR spectra of r(iGiGiCiC)(2), which reveal an A-form conformation. In MD simulations, r(CCGG)(2) retained A-form conformation 52% of the time, suggesting that its terminal base pairs may fray. The results indicate that revised backbone parameters improve predictions of RNA properties and that comparisons to measured sequence dependent thermodynamics provide useful benchmarks for testing force fields and computational methods.     

Assessment of biomolecular force fields for molecular dynamics simulations in a protein crystal

Different biomolecular force fields (OPLS‐AA, AMBER03, and GROMOS96) in conjunction with SPC, SPC/E and TIP3P water models are assessed for molecular dynamics simulations in a tetragonal lysozyme crystal. The root mean square deviations for the C_a atoms of lysozymes are about 0.1 to 0.2 nm from OPLS‐AA and AMBER03, smaller than 0.4 nm from GROMOS96. All force fields exhibit similar pattern in B‐factors, whereas OPLS‐AA and AMBER03 accurately reproduce experimental measurements. Despite slight variations, the primary secondary structures are well conserved using different force fields. Water diffusion in the crystal is approximately ten‐fold slower than in bulk phase. The directional and average water diffusivities from OPLS‐AA and AMBER03 along with SPC/E model match fairly well with experimental data. Compared to GROMOS96, OPLS‐AA and AMBER03 predict larger hydrophilic solvent‐accessible surface area of lysozyme, more hydrogen bonds between lysozyme and water, and higher percentage of water in hydration shell. SPC, SPC/E and TIP3P water models have similar performance in most energetic and structural properties, but SPC/E outperforms in water diffusion. While all force fields overestimate the mobility and electrical conductivity of NaCl, a combination of OPLS‐AA for lysozyme and the Kirkwood‐Buff model for ions is superior to others. As attributed to the steric restraints and surface interactions, the mobility and conductivity in the crystal are reduced by one to two orders of magnitude from aqueous solution. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Calculation of the electrostatic charges and energies for intercalation of aromatic drug molecules with DNA

In this work, we analyzed the influence of the charge model on the magnitudes of atomic charges and electrostatic energies for the binding of aromatic drug molecules with DNA. The dependence of the charge and energy on the level of theory (HF, DFT (B3LYP), MP2, semi‐empirical methods), basis set (STO‐3G, 3‐21G, 6‐31G, 6‐31G*, 6‐31G**), method of charge computation (Mulliken, Natural Population Analysis, CHelpG, Merz–Kollman), and force field charge (CHARMM27, AMBER99) has been tracked for typical aromatic drugs of different structure and charge state. Recommendations and restrictions have been formulated for the use of particular approaches in charge/electrostatic energy calculations. © 2009 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Interaction specific binding hotspots in Endonuclease colicin-immunity protein complex from MD simulations

The binding of Endonuclease colicin 9 (E9) by Immunity protein 9 (Im9) was found to involve some hotspots from helix III of Im9 on protein-protein interface that contribute the dominant binding energy to the complex.In the current work,MD simulations of the WT and three hotspot mutants (D51A,Y54A and Y55A of Im9) of the E9-Im9 complexes were carried out to investigate specific interaction mechanisms of these three hotspot residues.The changes of binding energy between the WT and mutants of the complex were computed by the MM/PBSA method using a polarized force field and were in excellent agreement with experiment values,verifying that these three residues were indeed hotspots of the binding complex.Energy decomposition analysis revealed that binding by D51 to E9 was dominated by electrostatic interaction due to the presence of the carboxyl group of Asp51 which hydrogen bonds to K89.For binding by hotspots Y54 and Y55,van der Waals interaction from the aromatic side chain of tyrosine provided the dominant interaction.For comparison,calculation by using the standard (nonpolarizable) AMBER99SB force field produced binding energy changes from these mutations in opposite direction to the experimental observation.Dynamic hydrogen bond analysis showed that conformations sampled from MD simulation in the standard AMBER force field were distorted from the native state and they disrupted the inter-protein hydrogen bond network of the protein-protein complex.The current work further demonstrated that electrostatic polarization plays a critical role in modulating protein-protein binding.     

Comparison of the AMBER,CHARMM, COMPASS,GROMOS, OPLS,TraPPE and UFF force fields for prediction of vapor–liquid coexistence curves and liquid densities

Monte Carlo simulations in the isobaric–isothermal and Gibbs ensembles are used to compute liquid densities and vapor–liquid coexistence curves for a series of small organic molecules for the AMBER-96, CHARMM22, COMPASS, GROMOS 43A1, OPLS-aa, TraPPE-UA, and UFF force fields. The simulation results are compared with experimental measurements to provide an assessment of the accuracy expected when using these force fields to study unknown molecules.     

Choice of Force Fields and Water Models for Sampling Solution Conformations of Bacteriophage T4 Lysozyme

A protein may exist as an ensemble of different conformations in solution, which cannot be represented by a single static structure. Molecular dynamics (MD) simulation has become a useful tool for sampling protein conformations in solution, but force fields and water models are important issues. This work presents a case study of the bacteriophage T4 lysozyme (T4L). We have found that MD simulations using a classic AMBER99SB force field and TIP4P water model cannot well describe hinge-bending domain motion of the wild-type T4L at the timescale of one microsecond. Other combinations, such as a residue-specific force field called RSFF2+ and a dispersion-corrected water model TIP4P-D, are able to sample reasonable solution conformations of T4L, which are in good agreement with experimental data. This primary study may provide candidates of force fields and water models for further investigating conformational transition of T4L.     

Rapid and accurate evaluation of the binding energies and the individual N-H···O=C, N-H···N, C-H···O=C, and C-H···N interaction energies for hydrogen-bonded peptide-base complexes

The binding energies of thirty-six hydrogen-bonded peptide-base complexes, including the peptide backbone-ase complexes and amino acid side chain-base complexes, are evaluated using the analytic potential energy function established in our lab recently and compared with those obtained from MP2, AMBER99, OPLSAA/L, and CHARMM27 calculations. The comparison indicates that the analytic potential energy function yields the binding energies for these complexes as reasonable as MP2 does, much better than the force fields do. The individual N H…O=C, N H…N, C H…O=C, and C H…N attractive interaction energies and C=O…O=C, N H…H N, C H…H N, and C H…H C repulsive interaction energies, which cannot be easily obtained from ab initio calculations, are calculated using the dipole-dipole interaction term of the analytic potential energy function. The individual N H…O=C, C H…O=C, C H…N attractive interactions are about 5.3±1.8, 1.2±0.4, and 0.8 kcal/mol, respectively, the individual N H … N could be as strong as about 8.1 kcal/mol or as weak as 1.0 kcal/mol, while the individual C=O…O=C, N H…H N, C H…H N, and C H…H C repulsive interactions are about 1.8±1.1, 1.7±0.6, 0.6±0.3, and 0.35±0.15 kcal/mol. These data are helpful for the rational design of new strategies for molecular recognition or supramolecular assemblies.