Assessing atomistic and coarse-grained force fields for protein-lipid interactions: the formidable challenge of an ionizable side chain in a membrane

Ionizable amino acid side chains play important roles in membrane protein structure and function, including the activation of voltage-gated ion channels, where it has been previously suggested that charged side chains may move through the hydrocarbon core of the membrane. However, all-atom molecular dynamics simulations have demonstrated large free energy barriers for such lipid-exposed motions. These simulations have also revealed that the membrane will deform due to the presence of a charged side chain, leading to a complex solvation microenvironment for which empirical force fields were not specifically parametrized. We have tested the ability of the all-atom CHARMM, Drude polarizable CHARMM, and a recent implementation of a coarse-grained force field to measure the thermodynamics of arginine-membrane interactions as a function of protonation state. We have employed model systems to attempt to match experimental bulk partitioning and quantum mechanical interactions within the membrane and found that free energy profiles from nonpolarizable and polarizable CHARMM simulations are accurate to within 1-2 kcal/mol. In contrast, the coarse-grained simulations failed to reproduce the same membrane deformations, exhibit interactions that are an order of magnitude too small, and thus, have incorrect free energy profiles. These results illustrate the need for careful parametrization of coarse-grained force fields and demonstrate the utility of atomistic molecular dynamics for providing quantitative thermodynamic and mechanistic analysis of protein-lipid interactions.     

Molecular dynamics simulation of hydrated DPPC monolayers using charge equilibration force fields

We present results of molecular dynamics simulations of a model DPPC-water monolayer using charge equilibration (CHEQ) force fields, which explicitly account for electronic polarization in a classical treatment of intermolecular interactions. The surface pressure, determined as the difference between the monolayer and pure water surface tensions at 323 K, is predicted to be 22.92 ±1.29 dyne/cm, just slightly below the broad range of experimental values reported for this system. The surface tension for the DPPC-water monolayer is predicted to be 42.35 ±1.16 dyne/cm, in close agreement with the experimentally determined value of 40.9 dyne/cm. This surface tension is also consistent with the value obtained from DPPC monolayer simulations using state-of-the-art nonpolarizable force fields. The current results of simulations predict a monolayer-water potential difference relative to the pure water-air interface of 0.64 ±0.02 Volts, an improved prediction compared to the fixed-charge CHARMM27 force field, yet still overestimating the experimental range of 0.3 to 0.45 Volts. As the charge equilibration model is a purely charge-based model for polarization, the current results suggest that explicitly modeled polarization effects can offer improvements in describing interfacial electrostatics in such systems.     

Structure, thermodynamics, and liquid-vapor equilibrium of ethanol from molecular-dynamics simulations using nonadditive interactions

We present a molecular-dynamics simulation study of the bulk and liquid-vapor interfacial properties of ethanol using a polarizable force field based on the fluctuating charge (FQ) formalism, as well as the nonpolarizable CHARMM22 force field. Both models are competitive with respect to the prediction of ambient liquid properties such as liquid density, enthalpy of vaporization, dielectric constant, and self-diffusion constants. The polarizable model predicts an average condensed-phase dipole moment of 2.2 D associated with an induced liquid-phase dipole moment of 0.6 D; though qualitatively in agreement with earlier nonadditive models as well as recent Car-Parinello calculations, the current FQ model underestimates the condensed-phase dipole moment. In terms of liquid structure, both models are in agreement with recent neutron-diffraction results of liquid ethanol structure, although the polarizable model predicts the hydroxyl-hydrogen-hydroxyl-hydrogen structure factor in closer agreement with the experimental data. In terms of interfacial properties, both models predict ambient surface tension to within 4% of the experimental value of 22.8 dyncm, while overestimating the surface excess entropy by almost a factor of 2. Both models display the characteristic preferential orientation of interfacial molecules. The polarizable model allows for a monotonic variation of the average molecular dipole moment from the bulk value to that of the vapor phase. Consequently, there is a dramatic difference in the surface potential predicted by the polarizable and nonpolarizable models. The polarizable model estimates a surface potential of -209+/-3 mV, while the nonpolarizable model yields a value of -944+/-10 mV. Finally, based on the vapor-liquid equilibrium simulation data from several temperatures, we estimate the critical properties of both models. As observed with other FQ models for associating fluids (such as water and methanol), and counter to what one would anticipate by modeling more physically the electrostatic response to local environment, the current FQ model underestimates the critical temperature and overestimates the critical density of ethanol; moreover, the FQ model is, in this respect, equivalent to the underlying fixed-charge model. These results further suggest the need to revisit polarizable models in terms of quantitative vapor-liquid equilibrium prediction.     

A theoretical study of aqueous solvation of K comparing ab initio, polarizable, and fixed-charge models

The hydration of K(+) is studied using a hierarchy of theoretical approaches, including ab initio Born-Oppenheimer molecular dynamics and Car-Parrinello molecular dynamics, a polarizable force field model based on classical Drude oscillators, and a nonpolarizable fixed-charge potential based on the TIP3P water model. While models based more directly on quantum mechanics offer the possibility to account for complex electronic effects, polarizable and fixed-charges force fields allow for simulations of large systems and the calculation of thermodynamic observables with relatively modest computational costs. A particular emphasis is placed on investigating the sensitivity of the polarizable model to reproduce key aspects of aqueous K(+), such as the coordination structure, the bulk hydration free energy, and the self diffusion of K(+). It is generally found that, while the simple functional form of the polarizable Drude model imposes some restrictions on the range of properties that can simultaneously be fitted, the resulting hydration structure for aqueous K(+) agrees well with experiment and with more sophisticated computational models. A counterintuitive result, seen in Car-Parrinello molecular dynamics and in simulations with the Drude polarizable force field, is that the average induced molecular dipole of the water molecules within the first hydration shell around K(+) is slightly smaller than the corresponding value in the bulk. In final analysis, the perspective of K(+) hydration emerging from the various computational models is broadly consistent with experimental data, though at a finer level there remain a number of issues that should be resolved to further our ability in modeling ion hydration accurately.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

Polarizable empirical force field for sulfur‐containing compounds based on the classical Drude oscillator model

Condensed‐phase computational studies of molecules using molecular mechanics approaches require the use of force fields to describe the energetics of the systems as a function of structure. The advantage of polarizable force fields over nonpolarizable (or additive) models lies in their ability to vary their electronic distribution as a function of the environment. Toward development of a polarizable force field for biological molecules, parameters for a series of sulfur‐containing molecules are presented. Parameter optimization was performed to reproduce quantum mechanical and experimental data for gas phase properties including geometries, conformational energies, vibrational spectra, and dipole moments as well as for condensed phase properties such as heats of vaporization, molecular volumes, and free energies of hydration. Compounds in the training set include methanethiol, ethanethiol, propanethiol, ethyl methyl sulfide, and dimethyl disulfide. The molecular volumes and heats of vaporization are in good accordance with experimental values, with the polarizable model performing better than the CHARMM22 nonpolarizable force field. Improvements with the polarizable model were also obtained for molecular dipole moments and in the treatment of intermolecular interactions as a function of orientation, in part due to the presence of lone pairs and anisotropic atomic polarizability on the sulfur atoms. Significant advantage of the polarizable model was reflected in calculation of the dielectric constants, a property that CHARMM22 systematically underestimates. The ability of this polarizable model to accurately describe a range of gas and condensed phase properties paves the way for more accurate simulation studies of sulfur‐containing molecules including cysteine and methionine residues in proteins. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Water Permeation Through DMPC Lipid Bilayers using Polarizable Charge Equilibration Force Fields

We investigate permeation energetics of water entering a model dimyristoylphosphatidylcholine (DMPC) bilayer via molecular dynamics simulations using polarizable Charge Equilibration (CHEQ) models. Potentials of mean force show 4.5-5.5 kcal/mol barriers for water permeation into bilayers. Barriers are highest when water coordination within the bilayer is prevented, and also when using force fields that accurately reproduce experimental alkane hydration free energies. The magnitude of the average water dipole moment decreases from 2.6 Debye (in bulk) to 1.88 Debye (in membrane interior). This variation correlates with the change in a water molecule's coordination number.     

Long-range Lennard-Jones and electrostatic interactions in interfaces: application of the isotropic periodic sum method

Molecular dynamics (MD) simulations of heptane/vapor, hexadecane/vapor, water/vapor, hexadecane/water, and dipalmitoylphosphatidylcholine (DPPC) bilayers and monolayers are analyzed to determine the accuracy of treating long-range interactions in interfaces with the isotropic periodic sum (IPS) method. The method and cutoff (rc) dependences of surface tensions, density profiles, water dipole orientation, and electrostatic potential profiles are used as metrics. The water/vapor, heptane/vapor, and hexadecane/vapor interfaces are accurately and efficiently calculated with 2D IPS (rc=10 A). It is demonstrated that 3D IPS is not practical for any of the interfacial systems studied. However, the hybrid method PME/IPS [Particle Mesh Ewald for electrostatics and 3D IPS for Lennard-Jones (LJ) interactions] provides an efficient way to include both types of long-range forces in simulations of large liquid/vacuum and all liquid/liquid interfaces, including lipid monolayers and bilayers. A previously published pressure-based long-range LJ correction yields results similar to those of PME/IPS for liquid/liquid interfaces. The contributions to surface tension of LJ terms arising from interactions beyond 10 A range from 13 dyn/cm for the hexadecane/vapor interface to approximately 3 dyn/cm for hexadecane/water and DPPC bilayers and monolayers. Surface tensions of alkane/vapor, hexadecane/water, and DPPC monolayers based on the CHARMM lipid force fields agree very well with experiment, whereas surface tensions of the TIP3P and TIP4P-Ew water models underestimate experiment by 16 and 11 dyn/cm, respectively. Dipole potential drops (DeltaPsi) are less sensitive to long-range LJ interactions than surface tensions. However, DeltaPsi for the DPPC bilayer (845+/-3 mV proceeding from water to lipid) and water (547+/-2 mV for TIP4P-Ew and 521+/-3 mV for TIP3P) overestimate experiment by factors of 3 and 5, respectively, and represent expected deficiencies in nonpolarizable force fields.     

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree–Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self‐consistent field method, Møller–Plesset perturbation theory method, and time‐dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self‐consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations‐All Atom (OPLS‐AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented. © 2013 Wiley Periodicals, Inc.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

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.     

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.     

The CHARMM–TURBOMOLE interface for efficient and accurate QM/MM molecular dynamics,free energies,and excited state properties

The quantum mechanical (QM)/molecular mechanical (MM) interface between Chemistry at HARvard Molecular Mechanics (CHARMM) and TURBOMOLE is described. CHARMM provides an extensive set of simulation algorithms, like molecular dynamics (MD) and free energy perturbation, and support for mature nonpolarizable and Drude polarizable force fields. TURBOMOLE provides fast QM calculations using density functional theory or wave function methods and excited state properties. CHARMM–TURBOMOLE is well‐suited for extended QM/MM MD simulations using first principles methods with large (triple‐ζ) basis sets. We demonstrate these capabilities with a QM/MM simulation of Mg²⁺(aq), where the MM outer sphere water molecules are represented using the SWM4‐NDP Drude polarizable force field and the ion and inner coordination sphere are represented using QM PBE, PBE0, and MP2 methods. The relative solvation free energies of Mg²⁺ and Zn²⁺ were calculated using thermodynamic integration. We also demonstrate the features for excited state properties. We calculate the time‐averaged solution absorption spectrum of indole, the emission spectrum of the indole excited state, and the electronic circular dichroism spectrum of an oxacepham. © 2014 Wiley Periodicals, Inc.     

Ion solvation thermodynamics from simulation with a polarizable force field

Thermodynamic measurements of the solvation of salts and electrolytes are relatively straightforward, but it is not possible to separate total solvation free energies into distinct cation and anion contributions without reference to an additional extrathermodynamic assumption. The present work attempts to resolve this difficulty using molecular dynamics simulations with the AMOEBA polarizable force field and perturbation techniques to directly compute absolute solvation free energies for potassium, sodium, and chloride ions in liquid water and formamide. Corresponding calculations are also performed with two widely used nonpolarizable force fields. The simulations with the polarizable force field accurately reproduce in vacuo quantum mechanical results, experimental ion-cluster solvation enthalpies, and experimental solvation free energies for whole salts, while the other force fields do not. The results indicate that calculations with a polarizable force field can capture the thermodynamics of ion solvation and that the solvation free energies of the individual ions differ by several kilocalories from commonly cited values.     

CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model

A fluctuating charge (FQ) force field is applied to molecular dynamics simulations for six small proteins in explicit polarizable solvent represented by the TIP4P-FQ potential. The proteins include 1FSV, 1ENH, 1PGB, 1VII, 1H8K, and 1CRN, representing both helical and beta-sheet secondary structural elements. Constant pressure and temperature (NPT) molecular dynamics simulations are performed on time scales of several nanoseconds, the longest simulations yet reported using explicitly polarizable all-atom empirical potentials (for both solvent and protein) in the condensed phase. In terms of structure, the FQ force field allows deviations from native structure up to 2.5 A (with a range of 1.0 to 2.5 A). This is commensurate to the performance of the CHARMM22 nonpolarizable model and other currently existing polarizable models. Importantly, secondary structural elements maintain native structure in general to within 1 A (both helix and beta-strands), again in good agreement with the nonpolarizable case. In qualitative agreement with QM/MM ab initio dynamics on crambin (Liu et al. Proteins 2001, 44, 484), there is a sequence dependence of average condensed phase atomic charge for all proteins, a dependence one would anticipate considering the differing chemical environments around individual atoms; this is a subtle quantum mechanical feature captured in the FQ model but absent in current state-of-the-art nonpolarizable models. Furthermore, there is a mutual polarization of solvent and protein in the condensed phase. Solvent dipole moment distributions within the first and second solvation shells around the protein display a shift towards higher dipole moments (increases on the order of 0.2-0.3 Debye) relative to the bulk; protein polarization is manifested via the enhanced condensed phase charges of typical polar atoms such as backbone carbonyl oxygens, amide nitrogens, and amide hydrogens. Finally, to enlarge the sample set of proteins, gas-phase minimizations and 1 ps constant temperature simulations are performed on various-sized proteins to compare to earlier work by Kaminsky et al. (J Comp Chem 2002, 23, 1515). The present work establishes the feasibility of applying a fully polarizable force field for protein simulations and demonstrates the approach employed in extending the CHARMM force field to include these effects.     

Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability: theory and applications

A current emphasis in empirical force fields is on the development of potential functions that explicitly treat electronic polarizability. In the present article, the commonly used methodologies for modeling electronic polarization are presented along with an overview of selected application studies. Models presented include induced point-dipoles, classical Drude oscillators, and fluctuating charge methods. The theoretical background of each method is followed by an introduction to extended Lagrangian integrators required for computationally tractable molecular dynamics simulations using polarizable force fields. The remainder of the review focuses on application studies using these methods. Emphasis is placed on water models, for which numerous examples exist, with a more thorough discussion presented on the recently published models associated with the Drude-based CHARMM and the AMOEBA force fields. The utility of polarizable models for the study of ion solvation is then presented followed by an overview of studies of small molecules (e.g., CCl₄, alkanes, etc.) and macromolecule (proteins, nucleic acids and lipid bilayers) application studies. The review is written with the goal of providing a general overview of the current status of the field and to facilitate future application and developments.     

Toward accurate solvation dynamics of lanthanides and actinides in water using polarizable force fields: from gas-phase energetics to hydration free energies

In this contribution, we focused on the use of polarizable force fields to model the structural, energetic, and thermodynamical properties of lanthanides and actinides in water. In a first part, we chose the particular case of the Th(IV) cation to demonstrate the capabilities of the AMOEBA polarizable force field to reproduce both reference ab initio gas-phase energetics and experimental data including coordination numbers and radial distribution functions. Using such model, we predicted the first polarizable force field estimate of Th(IV) solvation free energy, which accounts for −1,638 kcal/mol. In addition, we proposed in a second part of this work a full extension of the SIBFA (Sum of Interaction Between Fragments Ab initio computed) polarizable potential to lanthanides (La(III) and Lu(III)) and to actinides (Th(IV)) in water. We demonstrate its capabilities to reproduce all ab initio contributions as extracted from energy decomposition analysis computations, including many-body charge transfer and discussed its applicability to extended molecular dynamics and its parametrization on high-level post-Hartree–Fock data.     

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.