Popular Force Fields Research Articles

Comparing atomistic molecular mechanics force fields for a difficult target: a case study on the Alzheimer’s amyloid β-peptide

Macromolecular function arises from structure, and many diseases are associated with misfolding of proteins. Molecular simulation methods can augment experimental techniques to understand misfolding and aggregation pathways with atomistic resolution, but the reliability of these predictions is a function of the parameters used for the simulation. There are many biomolecular force fields available, but most are validated using stably folded structures. Here, we present the results of molecular dynamics simulations on the intrinsically disordered amyloid β-peptide (Aβ), whose misfolding and aggregation give rise to the symptoms of Alzheimer’s disease. Because of the link between secondary structure changes and pathology, being able to accurately model the structure of Aβ would greatly improve our understanding of this disease, and it may facilitate application of modeling approaches to other protein misfolding disorders. To this end, we compared five popular atomistic force fields (AMBER03, CHARMM22 + CMAP, GROMOS96 53A6, GROMOS96 54A7, and OPLS-AA) to determine which could best model the structure of Aβ. By comparing secondary structure content, NMR shifts, and radius of gyration to available experimental data, we conclude that AMBER03 and CHARMM22 + CMAP over-stabilize helical structure within Aβ, with CHARMM22 + CMAP also producing elongated Aβ structures, in conflict with experimental findings. OPLS-AA, GROMOS96 53A6, and GROMOS96 54A7 produce very similar results in terms of helical and β-strand content, calculated NMR shifts, and radii of gyration that agree well with experimental data.

Toward Optimized Potential Functions for Protein–Protein Interactions in Aqueous Solutions: Osmotic Second Virial Coefficient Calculations Using the MARTINI Coarse-Grained Force Field

Coarse-grained (CG) simulation methods are now widely used to model the structure and dynamics of large biomolecular systems. One important issue for using such methods - especially with regard to using them to model, for example, intracellular environments - is to demonstrate that they can reproduce experimental data on the thermodynamics of protein-protein interactions in aqueous solutions. To examine this issue, we describe here simulations performed using the popular coarse-grained MARTINI force field, aimed at computing the thermodynamics of lysozyme and chymotrypsinogen self-interactions in aqueous solution. Using molecular dynamics simulations to compute potentials of mean force between a pair of protein molecules, we show that the original parameterization of the MARTINI force field is likely to significantly overestimate the strength of protein-protein interactions to the extent that the computed osmotic second virial coefficients are orders of magnitude more negative than experimental estimates. We then show that a simple down-scaling of the van der Waals parameters that describe the interactions between protein pseudo-atoms can bring the simulated thermodynamics into much closer agreement with experiment. Overall, the work shows that it is feasible to test explicit-solvent CG force fields directly against thermodynamic data for proteins in aqueous solutions, and highlights the potential usefulness of osmotic second virial coefficient measurements for fully parameterizing such force fields.

Thermodynamic considerations for solubility and conformational transitions of poly-N-isopropyl-acrylamide

Thermodynamic considerations based on the free energy of hydration, free energy of solvation and partition coefficient predictions of the monomer N-isopropyl-acrylamide (NIPAM) determined using various intermolecular potentials are used to elucidate the origin of hydrophobicity/hydrophilicity across the lower critical solution temperature (LCST). Thermodynamic properties are predicted for NIPAM using adaptive bias force-molecular dynamics and various popular force-fields (AMBER, OPLS-AA, CHARMM and GROMOS) at four different temperatures: below the LCST (275 K and 300 K) and above the LCST (310 K and 330 K). The effect of changes in the thermodynamic properties of the monomer NIPAM at various temperatures below and above LCST on the kinetics of conformational transition of thermo-sensitive polymers is discussed. Our findings provide encouraging prospects for understanding the LCST transition in the hydrogel poly-N-isopropylacrylamide which shows temperature-dependent conformational changes in part due to the complex interplay of hydrophilic/hydrophobic interactions.

The Force Field Toolkit: Software for the Parameterization of Small Molecules from First Principles

The maturation of molecular dynamics (MD) through much technological advancement to extend the system size and simulation timescale remains the primary driving force in enabling the method to probe biological systems. An often-overlooked, but serious and long-persisting limitation, however, is the complexity of developing parameters for novel chemical matter, such as modified amino acids or small molecule ligands. To directly address this unmet need, we have developed the Force Field Toolkit (ffTK) for the development of CHARMM-compatible parameters from first principles with a specific focus on small molecules.ffTK is designed around the general principles of the CHARMM General Force Field (CGenFF), thus yielding parameters that are directly compatible with the popular CHARMM force field for biopolymers. Written entirely in Tcl/Tk and distributed free of charge as a VMD plugin, users interact with the toolkit via graphical user interface that formalizes each step of the parameterization, and organizes them into a clear, modular workflow. These steps include the optimization of partial atomic charges and bond, angle, and dihedral parameters against QM-calculated target data. ffTK is designed to be used in conjunction with the Gaussian software package for calculating all necessary QM target data, and provides appropriate tools to create and visually inspect all required Gaussian input and output files. Optimizations are performed using downhill or simulated annealing algorithms, and allows access to an array of settings that guide the optimization. A unique aspect of ffTK is the inclusion of embedded analysis utilities to quantitatively assess the progress of the optimization, and the performance of the resulting parameters. We have applied these tools to benchmark several classes of molecules against parameters available in the current CGenFF release, and will discuss the results herein.

Critical Assessment of Current Force Fields. Short Peptide Test Case.

The applicability of molecular dynamics simulations for studies of protein folding or intrinsically disordered proteins critically depends on quality of energetic functions-force fields. The four popular force fields for biomolecular simulations, CHARMM22/CMAP, AMBER FF03, AMBER FF99SB, and OPLS-AA/L, were compared in prediction of conformational propensities of all common proteinogenic amino acids. The minimalistic model of terminally block amino acids (dipeptides) was chosen for assessment of side chain effects on backbone propensities. The precise metadynamics simulations revealed striking inconsistency of trends in conformational preferences as manifested by investigated force fields for both backbone and side chains. To trace this disapproval between force fields, the two related AMBER force fields were studied more closely. In the cases of FF99SB and FF03, we uncovered that the distinct tends were driven by different charge models. Additionally, the effects of recent correction for side chain torsion (FF99SB-ILDN) were examined on affected amino acids and exposed significant coupling between free energy profiles and propensities of backbone and side chain conformers. These findings have important consequences for further force field development.

Solvation of complex surfaces via molecular density functional theory

We show that classical molecular density functional theory, here in the homogeneous reference fluid approximation in which the functional is inferred from the properties of the bulk solvent, is a powerful new tool to study, at a fully molecular level, the solvation of complex surfaces and interfaces by polar solvents. This implicit solvent method allows for the determination of structural, orientational, and energetic solvation properties that are on a par with all-atom molecular simulations performed for the same system, while reducing the computer time by two orders of magnitude. This is illustrated by the study of an atomistically-resolved clay surface composed of over a thousand atoms wetted by a molecular dipolar solvent. The high numerical efficiency of the method is exploited to carry a systematic analysis of the electrostatic and non-electrostatic components of the surface-solvent interaction within the popular Clay Force Field (CLAYFF). Solvent energetics and structure are found to depend weakly upon the atomic charges distribution of the clay surface, even for a rather polar solvent. We conclude on the consequences of such findings for force-field development.

Ion disturbance and clustering in the NaCl water solutions

Ion clustering and the solvation properties in the NaCl solutions are explored by molecular dynamics simulations with several popular force fields. The existence of ions has a negligible disturbance to the hydrogen bond structures and rotational mobility of water beyond the first ion solvation shells, which is suggested by the local hydrogen bond structures and the rotation times of water. The potential of mean force (PMF) of ion pair in the dilute solution presents a consistent view with the populations of ion clusters in the electrolyte solutions. The aggregation level of ions is sensitive to the force field used in the simulations. The ion-ion interaction potential plays an important role in the forming of the contact ion pair. The entropy of water increases as the ion pair approaches each other and the association of ion pair is driven by the increment of water entropy according to the results from the selected force fields. The kinetic transition from the single solvent separated state to the contact ion pair is controlled by the enthalpy loss of solution.

Structural Variations of Single and Tandem Mismatches in RNA Duplexes: A Joint MD Simulation and Crystal Structure Database Analysis

Internal loops within RNA duplex regions are formed by single or tandem basepairing mismatches with flanking canonical Watson-Crick basepairs on both sides. They are the most common motif observed in RNA secondary structures and play integral functional and structural roles. In this report, we have studied the structural features of 1 × 1, 2 × 2, and 3 × 3 internal loops using all-atom molecular dynamics (MD) simulation technique with explicit solvent model. As MD simulation is intricately dependent on the choice of force-field and these are often rather approximate, we have used both the most popular force-fields for nucleic acids-CHARMM27 and AMBER94-for a comparative analysis. We find that tandem noncanonical basepairs forming 2 × 2 and 3 × 3 internal loops are considerably more stable than the single mismatches forming 1 × 1 internal loops, irrespective of the force field. We have also analyzed crystal structure database to study the conservation of these helical fragments in the corresponding sets of RNA structures. We observe that the nature of stability in MD simulations mimic their fluctuating natures in crystal data sets also, probably indicating reliable natures of both the force fields to reproduce experimental results. We also notice significant structural changes in the wobble G:U basepairs present in these double helical stretches, leading to a biphasic stability for these wobble pairs to release the deformational strains introduced by internal loops within duplex regions.

Improving the Prediction of Absolute Solvation Free Energies Using the Next Generation OPLS Force Field.

Explicit solvent molecular dynamics free energy perturbation simulations were performed to predict absolute solvation free energies of 239 diverse small molecules. We use OPLS2.0, the next generation OPLS force field, and compare the results with popular small molecule force fields-OPLS_2005, GAFF, and CHARMm-MSI. OPLS2.0 produces the best correlation with experimental data (R(2) = 0.95, slope = 0.96) and the lowest average unsigned errors (0.7 kcal/mol). Important classes of compounds that performed suboptimally with OPLS_2005 show significant improvements.

Optimization of the OPLS-AA Force Field for Long Hydrocarbons.

The all-atom optimized potentials for liquid simulations (OPLS-AA) force field is a popular force field for simulating biomolecules. However, the current OPLS parameters for hydrocarbons developed using short alkanes cannot reproduce the liquid properties of long alkanes in molecular dynamics simulations. Therefore, the extension of OPLS-AA to (phospho)lipid molecules required for the study of biological membranes was hampered in the past. Here, we optimized the OPLS-AA force field for both short and long hydrocarbons. Following the framework of the OPLS-AA parametrization, we refined the torsional parameters for hydrocarbons by fitting to the gas-phase ab initio energy profiles calculated at the accurate MP2/aug-cc-pVTZ theory level. Additionally, the depth of the Lennard-Jones potential for methylene hydrogen atoms was adjusted to reproduce the densities and the heats of vaporization of alkanes and alkenes of different lengths. Optimization of partial charges finally allowed to reproduce the gel-to-liquid-phase transition temperature for pentadecane and solvation free energies. It is shown that the optimized parameter set (L-OPLS) yields improved hydrocarbon diffusion coefficients, viscosities, and gauche-trans ratios. Moreover, its applicability for lipid bilayer simulations is shown for a GMO bilayer in its liquid-crystalline phase.

Amino acid analogues bind to carbon nanotube via π-π interactions: Comparison of molecular mechanical and quantum mechanical calculations

Understanding the interaction between carbon nanotubes (CNTs) and biomolecules is essential to the CNT-based nanotechnology and biotechnology. Some recent experiments have suggested that the π-π stacking interactions between protein's aromatic residues and CNTs might play a key role in their binding, which raises interest in large scale modeling of protein-CNT complexes and associated π-π interactions at atomic detail. However, there is concern on the accuracy of classical fixed-charge molecular force fields due to their classical treatments and lack of polarizability. Here, we study the binding of three aromatic residue analogues (mimicking phenylalanine, tyrosine, and tryptophan) and benzene to a single-walled CNT, and compare the molecular mechanical (MM) calculations using three popular fixed-charge force fields (OPLSAA, AMBER, and CHARMM), with quantum mechanical (QM) calculations using the density-functional tight-binding method with the inclusion of dispersion correction (DFTB-D). Two typical configurations commonly found in π-π interactions are used, one with the aromatic rings parallel to the CNT surface (flat), and the other perpendicular (edge). Our calculations reveal that compared to the QM results the MM approaches can appropriately reproduce the strength of π-π interactions for both configurations, and more importantly, the energy difference between them, indicating that the various contributions to π-π interactions have been implicitly included in the van der Waals parameters of the standard MM force fields. Meanwhile, these MM models are less accurate in predicting the exact structural binding patterns (matching surface), meaning there are still rooms to be improved. In addition, we have provided a comprehensive and reliable QM picture for the π-π interactions of aromatic molecules with CNTs in gas phase, which might be used as a benchmark for future force field developments.

Comparison of selected polarizable and nonpolarizable water models in molecular dynamics simulations of ice Ih

We present a molecular dynamics simulation study in which we determined the melting point of ice I(h) for the polarizable SWM4-NDP water model (Lamoureux et al., Chem. Phys. Lett., 2006, 418, 245-249) and compared the performance of several popular water force fields, both polarizable and nonpolarizable, in terms of melting temperature, stability and orientational structuring of ice. The simulations yield the melting temperature of SWM4-NDP ice as low as T(m) = 185 ± 10 K, despite the quadrupole moment of a SWM4-NDP water molecule being close to the experimental gas phase value. The results thus show that the dependence of T(m) on the molecular quadrupole, observed for the three- and four-site water models, is generally lost if polarization is explicitly included. The study also shows that adding polarizability to a planar three-charge water model increases orientational disorder in hexagonal ice. In addition, analysis of the tetrahedral order in bulk ice reveals a correlation between the pre-existing degree of orientational disorder in ice simulated using different polarizable and nonpolarizable models and the melting temperature of the models. Our findings thus suggest some new considerations regarding the role of polarization forces in a crystalline solid that may guide future development of reliable polarizable water models for ice.

Assessment of Weak Intermolecular Interactions Across QM/MM Noncovalent Boundaries

An assessment of a number of quantum mechanical/molecular mechanical (QM/MM) combinations was performed for weak intermolecular interactions across noncovalent QM/MM 'boundaries'. The popular S22 data set, comprising of a number of weak hydrogen-bonded, dispersion-bound and complexes with mixed interactions was used for the assessment. A range of QM methods was combined with a number of popular MM force fields. The single-point interaction energies, at reference geometries, are presented as deviations to accurate CCSD(T)/CBS reference values. This investigation employed both additive and subtractive QM/MM schemes. The density functional has only a negligible effect; the choice of basis set was also negligible in terms of accuracy. The importance of selecting the most appropriate MM force field for accurately describing interactions at the noncovalent 'boundary' region has a dramatic effect on the accuracy.

Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant.

The chemical composition of small organic molecules is often very similar to amino acid side chains or the bases in nucleic acids, and hence there is no a priori reason why a molecular mechanics force field could not describe both organic liquids and biomolecules with a single parameter set. Here, we devise a benchmark for force fields in order to test the ability of existing force fields to reproduce some key properties of organic liquids, namely, the density, enthalpy of vaporization, the surface tension, the heat capacity at constant volume and pressure, the isothermal compressibility, the volumetric expansion coefficient, and the static dielectric constant. Well over 1200 experimental measurements were used for comparison to the simulations of 146 organic liquids. Novel polynomial interpolations of the dielectric constant (32 molecules), heat capacity at constant pressure (three molecules), and the isothermal compressibility (53 molecules) as a function of the temperature have been made, based on experimental data, in order to be able to compare simulation results to them. To compute the heat capacities, we applied the two phase thermodynamics method (Lin et al. J. Chem. Phys.2003, 119, 11792), which allows one to compute thermodynamic properties on the basis of the density of states as derived from the velocity autocorrelation function. The method is implemented in a new utility within the GROMACS molecular simulation package, named g_dos, and a detailed exposé of the underlying equations is presented. The purpose of this work is to establish the state of the art of two popular force fields, OPLS/AA (all-atom optimized potential for liquid simulation) and GAFF (generalized Amber force field), to find common bottlenecks, i.e., particularly difficult molecules, and to serve as a reference point for future force field development. To make for a fair playing field, all molecules were evaluated with the same parameter settings, such as thermostats and barostats, treatment of electrostatic interactions, and system size (1000 molecules). The densities and enthalpy of vaporization from an independent data set based on simulations using the CHARMM General Force Field (CGenFF) presented by Vanommeslaeghe et al. (J. Comput. Chem.2010, 31, 671) are included for comparison. We find that, overall, the OPLS/AA force field performs somewhat better than GAFF, but there are significant issues with reproduction of the surface tension and dielectric constants for both force fields.

Effect of Pore Topology and Accessibility on Gas Adsorption Capacity in Zeolitic−Imidazolate Frameworks: Bringing Molecular Simulation Close to Experiment

When all cages are assumed to be accessible, popular force fields such as universal force field (UFF) and DREIDING dramatically overpredicted gas adsorption capacity in two widely studied zeolitic−imidazolate frameworks (ZIFs), ZIF-68 and -69. Instead of adjusting the force-field parameters to match the experiments, herein we show that when the pore topology and accessibility are correctly taken into account, simulations with the standard force fields agree very well with the experiments. Careful inspection shows that ZIF-68 and -69 have two one-dimensional channels, which are not interaccessible to gases. The small channel consists of alternating small (HPR) and medium (GME) cages, while the large channel comprises the large (KNO) cages. Our analysis indicates that the small channel is not accessible to gases such as CO2. So when the cages in the small channel are intentionally blocked in our simulation, the predicted adsorption capacities of CO2, CH4 and N2 at room temperature from standard force-field parameters for the framework show excellent agreement with the experimental results. In the case of H2, all cages are accessible, so simulation results without cage-blocking show excellent agreement with experiment. Due to the promising potential of ZIFs in gas storage and separation, our work here shows that pore topology and accessibility should be carefully examined to understand how gases adsorb in ZIFs.

Physical Properties at the Base for the Development of an All-Atom Force Field for Ethylene Glycol

Ethylene glycol, the simplest of the diols, is a popular solvent, an antifreeze agent, a coolant, and a precursor in polymer production. In molecular modeling it is a model compound used to develop potentials for complex systems, like sugars. Despite the fact that many force fields for ethylene glycol exist in the literature, only few of them have been designed to reproduce the macroscopic properties of glycol and its mixtures, and rather more attention has been paid to the microscopic structure of the liquid. Those potentials that reproduce the properties accurately, apply also nonstandard fudge factors, therefore are not fully compatible with any popular force field. In this paper, we present a new potential for ethylene glycol, based on the OPLS all-atom force field and fully compatible with it, as well as with popular models for water. This potential is carefully validated against a broad range of physical properties measured experimentally and published in the literature. These properties include the density, expansion coefficient, compressibility, enthalpy of vaporization, surface tension, self-diffusion coefficient, and viscosity. Therefore, the potential presented here may be used in simulations of not only pure glycol but also mixtures with water, organic solvents, ionic liquids, phase interfaces, etc.

High-Resolution, Solvent-Free Coarse-Grained Model for Protein-Lipid Interactions

Many biophysical processes involving the interaction of proteins with membranes operate at time- and length-scales that are currently unattainable by all-atom computer simulations. To cope with this difficulty, increasingly more accurate and sophisticated coarse-grained models---both for proteins and lipids---are currently being developed.In this work, we combine two high-resolution, solvent-free coarse-grained models for proteins and lipids. Proteins are modeled by four beads per amino acid, providing enough backbone resolution to avoid explicit secondary structure bias towards the native state [Bereau and Deserno, J. Chem. Phys. 130, 235106 (2009)], while the lipid model was systematically tuned to reproduce the structural and mechanical properties of phosphocholine (PC) bilayers [Wang and Deserno, J. Phys. Chem. B 114, 11207 (2010); New J. Phys. 12, 095004 (2010)]. The transferrability of the two models across amino acid sequences and lipid species permits the investigation of a wide variety of scenarios, while the absence of explicit solvent allows for studies of large-scale phenomena.The two models were cross-parametrized to reproduce atomistic potential of mean force curves for the insertion of amino acids across the bilayer. We will illustrate different features of the model by simulating a small peptide which exhibits different stable folds in and out of the bilayer. Similarities and differences with the popular MARTINI force field will be discussed.

Force fields and molecular dynamics simulations

The objective of this review is to serve as an introductory guide for the non-expert to the exciting field of Molecular Dynamics (MD). MD simulations generate a phase space trajectory by integrating the classical equations of motion for a system of N particles. Here I review the basic concepts needed to understand the technique, what are the key elements to perform a simulation and which is the information that can be extracted from it. I will start defining what is a force field, which are the terms composing a classical force field, how the parameters of the potential are optimized, and which are the more popular force fields currently employed and the lines of research to improve them. Then the Molecular Dynamics technique will be introduced, including a general overview of the main algorithms employed to integrate the equations of motion, compute the long-range forces, work on different thermodynamic ensembles, or reduce the computational time. Finally the main properties that can be computed from a MD trajectory are briefly introduced.

Predicting hydration Gibbs energies of alkyl-aromatics using molecular simulation: a comparison of current force fields and the development of a new parameter set for accurate solvation data

The Gibbs energy of hydration is an important quantity to understand the molecular behavior in aqueous systems at constant temperature and pressure. In this work we review the performance of some popular force fields, namely TraPPE, OPLS-AA and Gromos, in reproducing the experimental Gibbs energies of hydration of several alkyl-aromatic compounds--benzene, mono-, di- and tri-substituted alkylbenzenes--using molecular simulation techniques. In the second part of the paper, we report a new model that is able to improve such hydration energy predictions, based on Lennard Jones parameters from the recent TraPPE-EH force field and atomic partial charges obtained from natural population analysis of density functional theory calculations. We apply a scaling factor determined by fitting the experimental hydration energy of only two solutes, and then present a simple rule to generate atomic partial charges for different substituted alkyl-aromatics. This rule has the added advantages of eliminating the unnecessary assumption of fixed charge on every substituted carbon atom and providing a simple guideline for extrapolating the charge assignment to any multi-substituted alkyl-aromatic molecule. The point charges derived here yield excellent predictions of experimental Gibbs energies of hydration, with an overall absolute average deviation of less than 0.6 kJ mol(-1). This new parameter set can also give good predictive performance for other thermodynamic properties and liquid structural information.

Helix Formation in a Pentapeptide: Experiment and Force-field Dependent Dynamics

We used a combined approach of experiment and simulation to determine the helical population and folding pathway of a small helix forming blocked pentapeptide, Ac-(Ala)(5)-NH(2). Experimental structural characterization of this blocked peptide was carried out with far UV circular dichroism spectroscopy, FTIR, and NMR measurements. These measurements confirm the presence of the α-helical state in a buffer solution. Direct molecular dynamics and replica-exchange simulations of the pentapeptide were performed using several popular force fields with explicit solvent. The simulations yielded statistically reliable estimates of helix populations, melting curves, folding, and nucleation times. The distributions of conformer populations are used to measure folding cooperativity. Finally, a statistical analysis of the sample of helix-coil transition paths was performed. The details of the calculated helix populations, folding kinetics and pathways vary with the employed force field. Interestingly, the helix populations, folding, and unfolding times obtained from most of the studied force fields are in qualitative agreement with each other and with available experimental data, with the deviations corresponding to several kcal/mol in energy at 300 K. Most of the force fields also predict qualitatively similar transition paths, with unfolding initiated at the C-terminus. Accuracy of potential energy parameters, rather than conformational sampling may be the limiting factor in current molecular simulations.