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.     

pmx: Automated protein structure and topology generation for alchemical perturbations

Computational protein design requires methods to accurately estimate free energy changes in protein stability or binding upon an amino acid mutation. From the different approaches available, molecular dynamics‐based alchemical free energy calculations are unique in their accuracy and solid theoretical basis. The challenge in using these methods lies in the need to generate hybrid structures and topologies representing two physical states of a system. A custom made hybrid topology may prove useful for a particular mutation of interest, however, a high throughput mutation analysis calls for a more general approach. In this work, we present an automated procedure to generate hybrid structures and topologies for the amino acid mutations in all commonly used force fields. The described software is compatible with the Gromacs simulation package. The mutation libraries are readily supported for five force fields, namely Amber99SB, Amber99SB*‐ILDN, OPLS‐AA/L, Charmm22*, and Charmm36. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Origin of intrinsic 3(10)-helix versus strand stability in homopolypeptides and its implications for the accuracy of the Amber force field

Current all-atom force fields often fail to recognize the native structure of a protein as the lowest free energy minimum. One possible cause could be the mathematical form of the potential based on the assumption that the conformation of a residue is independent of its neighbors. Here, using quantum mechanical (QM) methods (MP2/6-31g**//HF/6-31g** and MP2/cc-pVDZ//cc-pVDZ//HF/cc-pVDZ), the intrinsic correctness of the gas phase terms (without solvation) of the Amber ff03 and ff99 potentials are examined by testing their ability to reproduce the relative 3(10)-helix versus extended structure stabilities in the gas phase for 1-7-residue alanine, valine, leucine, and isoleucine homopolypeptides. The 3(10)-helix versus extended state stability strongly depends on chain length and less on the amino acid identity. The helical conformation becomes lower in energy than the extended conformation for all tested peptides longer than two residues, and its stability increases with the increase of chain length. The ff03 potential better describes the 3(10)-helix versus extended state energy than ff99 and also reproduces the curvature of the relative helix-extended state energies. Therefore, the mathematical form of the Amber potential is sufficient to describe the local effect of 3(10)-helix versus extended structure stabilization in the gas phase. However, the energy curves are shifted and the backbone geometries differ compared with the QM results. This may cause significant geometric discrepancies between native and predicted structures. Therefore, extant molecular mechanics force fields, such as Amber, need refinement of their parameters to correctly describe helix-extended state energetics and geometry of major conformations.     

Symbasis of molar volumes and the cohesion energies of organic polymers and clusters of their repeat units

A method for evaluating the molar volume and cohesion energy of a polymer when its structure is a set of clusters of repeat units with hydrogen atoms terminating the broken bonds is proposed. The properties of 40 polymers were shown via Monte Carlo methods and the molecular-mechanics method with the force fields Amber99 and MM+ with corrections for the formation of hydrogen and CN dipole-pair bonds. The results agree with the results of other calculation methods and with the literature data.     

Development and testing of a general amber force field

We describe here a general Amber force field (GAFF) for organic molecules. GAFF is designed to be compatible with existing Amber force fields for proteins and nucleic acids, and has parameters for most organic and pharmaceutical molecules that are composed of H, C, N, O, S, P, and halogens. It uses a simple functional form and a limited number of atom types, but incorporates both empirical and heuristic models to estimate force constants and partial atomic charges. The performance of GAFF in test cases is encouraging. In test I, 74 crystallographic structures were compared to GAFF minimized structures, with a root-mean-square displacement of 0.26 A, which is comparable to that of the Tripos 5.2 force field (0.25 A) and better than those of MMFF 94 and CHARMm (0.47 and 0.44 A, respectively). In test II, gas phase minimizations were performed on 22 nucleic acid base pairs, and the minimized structures and intermolecular energies were compared to MP2/6-31G* results. The RMS of displacements and relative energies were 0.25 A and 1.2 kcal/mol, respectively. These data are comparable to results from Parm99/RESP (0.16 A and 1.18 kcal/mol, respectively), which were parameterized to these base pairs. Test III looked at the relative energies of 71 conformational pairs that were used in development of the Parm99 force field. The RMS error in relative energies (compared to experiment) is about 0.5 kcal/mol. GAFF can be applied to wide range of molecules in an automatic fashion, making it suitable for rational drug design and database searching.     

A molecular dynamics study of the thermal response of crystalline cellulose Iβ

Molecular dynamics simulations were performed to better understand the atomic details of thermal induced transitions in cellulose Iβ. The latest version of the GLYCAM force field series (GLYCAM06) was used for the simulations. The unit cell parameters, density, torsion angles and hydrogen-bonding network of the crystalline polymer were carefully analyzed. The simulated data were validated against the experimental results obtained by X-ray diffraction for the crystal structure of cellulose Iβ at room and high temperatures, as well as against the temperature-dependent IR measurements describing the variation of hydrogen bonding patterns. Distinct low and high temperature structures were identified, with a phase transition temperature of 475–500 K. In the high-temperature structure, all the origin chains rotated around the helix axis by about 30° and the conformation of all hydroxymethyl groups changed from tg to either gt on origin chains or gg on center chains. The hydrogen-bonding network was reorganized along with the phase transition. Compared to the previously employed GROMOS 45a4 force field, GLYCAM06 yields data in much better agreement with experimental observations, which reflects that a cautious parameterization of the nonbonded interaction terms in a force field is critical for the correct prediction of the thermal response in cellulose crystals.     

Implementation of a protein reduced point charge model toward molecular dynamics applications

A reduced point charge model was developed in a previous work from the study of extrema in smoothed charge density distribution functions generated from the Amber99 molecular electrostatic potential. In the present work, such a point charge distribution is coupled with the Amber99 force field and implemented in the program TINKER to allow molecular dynamics (MD) simulations of proteins. First applications to two polypeptides that involve α-helix and β-sheet motifs are analyzed and compared to all-atom MD simulations. Two types of coarse-grained (CG)-based trajectories are generated using, on one hand, harmonic bond stretching terms and, on the other hand, distance restraints. Results show that the use of the unrestrained CG conditions are sufficient to preserve most of the secondary structure characteristics but restraints lead to a better agreement between CG and all-atom simulation results such as rmsd, dipole moment, and time-dependent mean square deviation functions.     

Carboxymethylated cyclodextrin derivatives as chiral NMR discriminating agents

Procedures to prepare cyclodextrins with carboxymethyl groups incorporated selectively at the primary (6-position) or secondary (2-position) are described. Complexation properties of the primary and secondary carboxymethylated derivatives of α-, β-, and γ-cyclodextrins are compared to native cyclodextrins and indiscriminately substituted carboxymethylated cyclodextrins, using pheniramine, chlorpheniramine, and brompheniramine as substrates. The stoichiometry of association of these substrates with the α-cyclodextrins is 1:1, whereas with the γ-cyclodextrins, a 2:1 substrate:cyclodextrin complex forms. Data for the β-cyclodextrins suggest that there is a mix of 1:1 and 2:1 substrate–cyclodextrin complexes. The position of the carboxymethyl groups on the cyclodextrin does not appear to alter the geometry of substrate–cyclodextrin association. The effectiveness of the carboxymethylated cyclodextrins as chiral NMR discriminating agents is compared with the native cyclodextrins. In all cases, the indiscriminately substituted α-, β-, and γ-cyclodextrins are more effective at enantiodistinction with the cationic substrates than native cyclodextrins or the derivatives with carboxymethyl groups at the primary or secondary positions. Among α-, β-, and γ-indiscriminately substituted cyclodextrins, there was no clearly optimal candidate for chiral NMR discrimination studies. The indiscriminately substituted carboxymethyl cyclodextrins are effective water-soluble chiral NMR discrimination reagents for cationic substrates.     

The influence of solvation on π *←n and π *←π transition energies in molecules containing Aza or carbonyl groups

The influence of solvation on ^*n and ^* transition energies of formaldehyde, acetaldehyde, acetone, pyridine and 1,2-, 1,3-, 1,4-diazabenzenes has been investigated through CNDO calculations.A static solvation model which distinguishes a) molecules directly involved in hydrogen bonding with solute, b) the layer of molecules in contact with the solute molecule and c) the main of molecules farther from solute, is presented.Blue and red shifts due to solvent effects are correctly predicted by calculations for each model.     

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.     

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.     

Solvation properties of N-acetyl-β-glucosamine: molecular dynamics study incorporating electrostatic polarization

N-Acetyl-β-glucosamine (NAG) is an important moiety of glycoproteins and is involved in many biological functions. However, conformational and dynamical properties of NAG molecules in aqueous solution, the most common biological environment, remain ambiguous due to limitations of experimental methods. Increasing efforts are made to probe structural properties of NAG and NAG-containing macromolecules, like peptidoglycans and polymeric chitin, at the atomic level using molecular dynamics simulations. In this work, we develop a polarizable carbohydrate force field for NAG and contrast simulation results of various properties using this novel force field and an analogous nonpolarizable (fixed charge) model. Aqueous solutions of NAG and its oligomers are investigated; we explore conformational properties (rotatable bond geometry), electrostatic properties (dipole moment distribution), dynamical properties (self-diffusion coefficient), hydrogen bonding (water bridge structure and dynamics), and free energy of hydration. The fixed-charge carbohydrate force field exhibits deviations from the gas phase relative rotation energy of exocyclic hydroxymethyl side chain and of chair/boat ring distortion. The polarizable force field predicts conformational properties in agreement with corresponding first-principles results. NAG-water hydrogen bonding pattern is studied through radial distribution functions (RDFs) and correlation functions. Intermolecular hydrogen bonding between solute and solvent is found to stabilize NAG solution structures while intramolecular hydrogen bonds define glycosidic linkage geometry of NAG oligomers. The electrostatic component of hydration free energy is highly dependent on force field atomic partial charges, influencing a more favorable free energy of hydration in the fixed-charge model compared to the polarizable model.     

Force field‐dependent solution properties of glycine oligomers

Molecular simulations can be used to study disordered polypeptide systems and to generate hypotheses on the underlying structural and thermodynamic mechanisms that govern their function. As the number of disordered protein systems investigated with simulations increase, it is important to understand how particular force fields affect the structural properties of disordered polypeptides in solution. To this end, we performed a comparative structural analysis of Gly₃ and Gly₁₀ in aqueous solution from all atom, microsecond molecular dynamics (MD) simulations using the CHARMM 27 (C27), CHARMM 36 (C36), and Amber ff12SB force fields. For each force field, Gly₃ and Gly₁₀ were simulated for at least 300 ns and 1 μs, respectively. Simulating oligoglycines of two different lengths allows us to evaluate how force field effects depend on polypeptide length. Using a variety of structural metrics (e.g., end‐to‐end distance, radius of gyration, dihedral angle distributions), we characterize the distribution of oligoglycine conformers for each force field and show that each sample conformation space differently, yielding considerably different structural tendencies of the same oligoglycine model in solution. Notably, we find that C36 samples more extended oligoglycine structures than both C27 and ff12SB. © 2015 Wiley Periodicals, Inc.     

Investigating the quartz (1010)/water interface using classical and ab initio molecular dynamics

Two different terminations of the (1010) surface of quartz (α and β) interacting with water are simulated by classical (CMD) (using two different force fields) and ab initio molecular dynamics (AIMD) and compared with previously published X-ray reflectivity (XR) experiments. Radial distribution functions between hydroxyl and water show good agreement between AIMD and CMD using the ClayFF force field for both terminations. The Lopes et al. (Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D. J. Phys. Chem. B2006, 110, 2782-2792) force field (LFF), however, underestimates the extent of hydroxyl-water hydrogen bonding. The β termination is found to contain hydroxyl-hydroxyl hydrogen bonds; the quartz surface hydroxyl hydrogens and oxygens that hydrogen bond with each other exhibit greatly reduced hydrogen bonding to water. Conversely, the hydroxyl hydrogen and oxygens that are not hydrogen bonded to other surface hydroxyls but are connected to those that are show a considerable amount of hydrogen bonding to water. The electron density distribution of an annealed surface of quartz (1010) obtained by XR is in qualitative agreement with electron densities calculated by CMD and AIMD. In all simulation methods, the interfacial water peak appears farther from the surface than observed by XR. Agreement among AIMD, LFF, and XR is observed for the relaxation of the near-surface atoms; however, ClayFF shows a larger discrepancy. Overall, results show that for both terminations of (1010), LFF treats the near-surface structure more accurately whereas ClayFF treats the interfacial water structure more accurately. It is shown that the number of hydroxyl and water hydrogen bonds to the bridging Si-O-Si oxygens connecting the surface silica groups to the rest of the crystal is much greater for the α than the β termination. It is suggested that this may play a role in the greater resistance to dissolution of the β termination than that of the α termination.     

Less is more when simulating unsulfated glycosaminoglycan 3D‐structure: Comparison of GLYCAM06/TIP3P,PM3‐CARB1/TIP3P,and SCC‐DFTB‐D/TIP3P predictions with experiment

The 3D‐structure of extracellular matrix glycosaminoglycans is central to function, but is currently poorly understood. Resolving this will provide insight into their heterogeneous biological roles and help to realize their significant therapeutic potential. Glycosaminoglycan chemical isoforms are too numerous to study experimentally and simulation provides a tractable alternative. However, best practice for accurate calculation of glycosaminoglycan 3D‐structure within biologically relevant nanosecond timescales is uncertain. Here, we evaluate the ability of three potentials to reproduce experimentally observed glycosaminoglycan monosaccharide puckering, disaccharide 3D‐conformation, and characteristic solvent interactions. Temporal dynamics of unsulfated chondroitin, chondroitin‐4‐sulfate, and hyaluronan β(1→3) disaccharides were simulated within TIP3P explicit solvent unrestrained for 20 ns using the GLYCAM06 force‐field and two semi‐empirical quantum mechanics methods, PM3‐CARB1 and SCC‐DFTB‐D (both within a hybrid QM/MM formalism). Comparison of calculated and experimental properties (vicinal couplings, nuclear Overhauser enhancements, and glycosidic linkage geometries) showed that the carbohydrate‐specific parameterization of PM3‐CARB1 imparted quantifiable benefits on monosaccharide puckering and that the SCC‐DFTB‐D method (including an empirical correction for dispersion) best modeled the effects of hexosamine 4‐sulfation. However, paradoxically, the most approximate approach (GLYCAM06/TIP3P) was the best at predicting monosaccharide puckering, 3D‐conformation, and solvent interactions. Our data contribute to the debate and emerging consensus on the relative performance of these levels of theory for biological molecules. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Force field‐dependant structural divergence revealed during long time simulations of Calbindin d9k

The structural and the dynamic features of the Calbindin (CaB) protein in its holo and apo states are compared using molecular dynamics simulations under nine different force fields (FFs) (G43a1, G53a6, Opls‐AA, Amber94, Amber99, Amber99p, AmberGS, AmberGSs, and Amber99sb). The results show that most FFs reproduce reasonably well the majority of the experimentally derived features of the CaB protein. However, in several cases, there are significant differences in secondary structure properties, root mean square deviations (RMSDs), root mean square fluctuations (RMSFs), and S² order parameters among the various FFs. What is more, in certain cases, these parameters differed from the experimentally derived values. Some of these deviations became noticeable only after 50 ns. A comparison with experimental data indicates that, for CaB, the Amber94 shows overall best agreement with the measured values, whereas several others seem to deviate from both crystal and nuclear magnetic resonance data. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Direct evidence for symmetry control in cyclodextrin-water interactions

Abstract

The cyclodextrins, cyclic oligosaccharides possessing generally 6, 7, or 8 α, 1 –4 linked glucopyranose units, are widely used for the solubilisation and transport of organic molecules in aqueous media.³ Their solubilities 145 gL^?1, α-Cd; 20 gL^?1, β-CD; 220 gL^?1, γ-CD, are considerably lower than those of simple saccharides and in particular that of β-CD, the most widely available compound, is anomalously low. We have previously shown that the solubilities arise from symmetry determined interactions of the CDs with the dynamic hexagonal structure of water. In the case of the unmodified cyclodextrins all three compounds exist in solution as large aggregates, thus the solubility of these systems is controlled by the interactions between these hydrogen bonded aggregates and the hydrogen bonding networks present in water. In the case of α- or γ-CD there are favorable overlaps with the hexagonal water structure, however for the seven-fold symmetry of β-CD no such favorable interactions can occur between the odd symmetry element and the even symmetry elements of water. Hence the observed solubilities of α- and γ-CD are higher than that of β-CD.