Many-body polarization effects and the membrane dipole potential

Drude Polarizable Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Saturated and Monounsaturated Zwitterionic Lipids

Accurate empirical force fields of lipid molecules are a critical component of molecular dynamics simulation studies aimed at investigating properties of monolayers, bilayers, micelles, vesicles, and liposomes, as well as heterogeneous systems, such as protein-membrane complexes, bacterial cell walls, and more. While the majority of lipid force field-based simulations have been performed using pairwise-additive nonpolarizable models, advances have been made in the development of the polarizable force field based on the classical Drude oscillator model. In the present study, we undertake further optimization of the Drude lipid force field, termed Drude2023, including improved treatment of the phosphate and glycerol linker region of PC and PE headgroups, additional optimization of the alkene group in monounsaturated lipids, and inclusion of long-range Lennard-Jones interactions using the particle-mesh Ewald method. Initial optimization targeted quantum mechanical (QM) data on small model compounds representative of the linker region. Subsequent optimization targeted QM data on larger model compounds, experimental data, and dihedral potentials of mean force from the CHARMM36 additive lipid force field using a parameter reweighting protocol. The use of both experimental and QM target data during the reweighting protocol is shown to produce physically reasonable parameters that reproduce a collection of experimental observables. Target data for optimization included surface area/lipid for DPPC, DSPC, DMPC, and DLPC bilayers and nuclear magnetic resonance (NMR) order parameters for DPPC bilayers. Validation data include prediction of membrane thickness, scattering form factors, electrostatic potential profiles, compressibility moduli, surface area per lipid, water permeability, NMR T1 relaxation times, diffusion constants, and monolayer surface tensions for a variety of saturated and unsaturated lipid mono- and bilayers. Overall, the agreement with experimental data is quite good, though the results are less satisfactory for the NMR T1 relaxation times for carbons near the ester groups. Notable improvements compared to the additive C36 force field were obtained for membrane dipole potentials, lipid diffusion coefficients, and water permeability with the exception of monounsaturated lipid bilayers. It is anticipated that the optimized polarizable Drude2023 force field will help generate more accurate molecular simulations of pure bilayers and heterogeneous systems containing membranes, advancing our understanding of the role of electronic polarization in these systems.

The lung surfactant activity probed with molecular dynamics simulations

The surface of pulmonary alveolar subphase is covered with a mixture of lipids and proteins. This lung surfactant plays a crucial role in lung functioning. It shows a complex phase behavior which can be altered by the interaction with third molecules such as drugs or pollutants. For studying multicomponent biological systems, it is of interest to couple experimental approach with computational modelling yielding atomic-scale information. Simple two, three, or four-component model systems showed to be useful for getting more insight in the interaction between lipids, lipids and proteins or lipids and proteins with drugs and impurities. These systems were studied theoretically using molecular dynamic simulations and experimentally by means of the Langmuir technique. A better understanding of the structure and behavior of lung surfactants obtained from this research is relevant for developing new synthetic surfactants for efficient therapies, and may contribute to public health protection.

Harnessing Deep Learning for Optimization of Lennard-Jones Parameters for the Polarizable Classical Drude Oscillator Force Field

The outcomes of computational chemistry and biology research, including drug design, are significantly influenced by the underlying force field (FF) used in molecular simulations. While improved FF accuracy may be achieved via inclusion of explicit treatment of electronic polarization, such an extension must be accompanied by optimization of van der Waals (vdW) interactions, in the context of the Lennard-Jones (LJ) formalism in the present study. This is particularly challenging due to the extensive nature of chemical space combined with the correlated nature of LJ parameters. To address this challenge, a deep learning (DL)-based parametrization framework is developed, allowing for sampling of wide ranges of LJ parameters targeting experimental condensed phase thermodynamic properties. The present work utilizes this framework to develop the LJ parameters for atoms associated with four distinct groups covering 10 different atom types. Final parameter selection was facilitated by quantum mechanical data on rare-gas interactions with the training set molecules. The chosen parameters were then validated through experimental hydration free energies and condensed phase thermodynamic properties of validation set molecules to confirm transferability. The ultimate outcome of utilizing this framework is a set of LJ parameters in the context of the polarizable Drude FF, which demonstrated improvement in the reproduction of both experimental pure solvent and crystal properties and hydration free energies of the molecules compared to the additive CHARMM General FF (CGenFF) including the ability of the Drude FF to accurately reproduce both experimental pure solvent properties and hydration free energies. The study also shows how correlations between difference in the reproduction of condensed phase data between model compounds may be used to direct the selection of new atom types and training set molecules during FF development.

A simple fragment-based method for van der Waals corrections over density functional theory

Modeling intermolecular noncovalent interactions between large molecules remains a challenge for the electron structure theory community due to the high cost. Fragment-based methods usually fare well in reducing the cost...

Preparing and Analyzing Polarizable Molecular Dynamics Simulations with the Classical Drude Oscillator Model

Molecular dynamics (MD) simulations performed with force fields that include explicit electronic polarization are becoming more prevalent in the field. The increasing emergence of these simulations is a result of continual refinement against a range of theoretical and empirical target data, optimization of software algorithms for higher performance, and availability of graphical processing unit hardware to further accelerate the simulations. Polarizable MD simulations are likely to be most impactful in biomolecular systems in which heterogeneous environments or unique microenvironments exist that would lead to inaccuracies in simulations performed with fixed-charge, nonpolarizable force fields. The further adoption of polarizable MD simulations will benefit from tutorial material that specifically addresses preparing and analyzing their unique features. In this chapter, we introduce common protocols for preparing routine biomolecular systems containing proteins, including both a globular protein in aqueous solvent and a transmembrane model peptide in a phospholipid bilayer. Details and example input files are provided for preparation of the simulation system using CHARMM, performing the simulations with OpenMM, and analyzing interesting dipole moment properties in CHARMM.

An ω-3, but Not an ω-6 Polyunsaturated Fatty Acid Decreases Membrane Dipole Potential and Stimulates Endo-Lysosomal Escape of Penetratin

Although the largely positive intramembrane dipole potential (DP) may substantially influence the function of transmembrane proteins, its investigation is deeply hampered by the lack of measurement techniques suitable for high-throughput examination of living cells. Here, we describe a novel emission ratiometric flow cytometry method based on F66, a 3-hydroxiflavon derivative, and demonstrate that 6-ketocholestanol, cholesterol and 7-dehydrocholesterol, saturated stearic acid (SA) and ω-6 γ-linolenic acid (GLA) increase, while ω-3 α-linolenic acid (ALA) decreases the DP. These changes do not correlate with alterations in cell viability or membrane fluidity. Pretreatment with ALA counteracts, while SA or GLA enhances cholesterol-induced DP elevations. Furthermore, ALA (but not SA or GLA) increases endo-lysosomal escape of penetratin, a cell-penetrating peptide. In summary, we have developed a novel method to measure DP in large quantities of individual living cells and propose ALA as a physiological DP lowering agent facilitating cytoplasmic entry of penetratin.

Molecular Dynamics Simulations Based on Polarizable Models Show that Ion Permeation Interconverts between Different Mechanisms as a Function of Membrane Thickness

Different mechanisms have been proposed to explain the permeation of charged compounds through lipid membranes. Overall, it is expected that an ion-induced defect permeation mechanism, where substantial membrane deformations accompany ion movement, should be dominant in thin membranes but that a solubility-diffusion mechanism, where ions partition into the membrane core with large associated dehydration energy costs, becomes dominant in thicker membranes. However, while this physical picture is intuitively reasonable, capturing the interconversion between these two permeation mechanisms in molecular dynamics (MD) simulations based on atomic models is challenging. In particular, simulations relying on nonpolarizable force fields are artificially unfavorable to the solubility-diffusion mechanism, as induced polarization of the nonpolar hydrocarbon is ignored, causing overestimated free energy costs for charged molecules to enter into this region of the membrane. In this study, all-atom MD simulations based on nonpolarizable and polarizable force fields are used to quantitatively characterize the permeation process for the arginine side chain analog methyl-guanidinium through bilayer membranes of mono-unsaturated phosphatidylcholine lipids with and without cholesterol, resulting in thicknesses spanning from ∼24 to ∼42 Å. With simulations based on a nonpolarizable force field, ion translocation can take place solely through an ion-induced defect mechanism, with free energy barriers increasing linearly from 14 to 40 kcal/mol, depending on the thickness. However, with simulations based on a polarizable force field, ion translocation is predominantly dominated by an ion-induced defect mechanism in thin membranes, which progressively converts to a solubility-diffusion mechanism as the membranes get thicker. The transition between the two mechanisms occurs at a thickness of ∼29 Å, with lipid tails of 22 or more carbon atoms. This situation appears to represent the upper limit for ion-induced defect permeation within the current polarizable models. Beyond this thickness, it becomes energetically preferable for the ion to dehydrate and partition into the membrane core-a phenomenon that cannot be captured using the nonpolarizable models. Induced electronic polarizability therefore leads not just to a shift in permeation energetics but to an interconversion between two strikingly different physical mechanisms. The result highlights the importance of induced polarizability in modeling lipid membranes.

Statistical mechanics of polarizable force fields based on classical Drude oscillators with dynamical propagation by the dual-thermostat extended Lagrangian

Polarizable force fields based on classical Drude oscillators offer a practical and computationally efficient avenue to carry out molecular dynamics (MD) simulations of large biomolecular systems. To treat the polarizable electronic degrees of freedom, the Drude model introduces a virtual charged particle that is attached to its parent nucleus via a harmonic spring. Traditionally, the need to relax the electronic degrees of freedom for each fixed set of nuclear coordinates is achieved by performing an iterative self-consistent field (SCF) calculation to satisfy a selected tolerance. This is a computationally demanding procedure that can increase the computational cost of MD simulations by nearly one order of magnitude. To avoid the costly SCF procedure, a small mass is assigned to the Drude particles, which are then propagated as dynamic variables during the simulations via a dual-thermostat extended Lagrangian algorithm. To help clarify the significance of the dual-thermostat extended Lagrangian propagation in the context of the polarizable force field based on classical Drude oscillators, the statistical mechanics of a dual-temperature canonical ensemble is formulated. The conditions for dynamically maintaining the dual-temperature properties in the case of the classical Drude oscillator are analyzed using the generalized Langevin equation.

Interfacial Water Structure at Zwitterionic Membrane/Water Interface: The Importance of Interactions between Water and Lipid Carbonyl Groups

In this work, atomistic molecular dynamics (MD) simulations of palmitoyl-oleoyl-phosphatidylcholine (POPC) bilayer were carried out to investigate the effect of water models on membrane dipole potential, which is primarily associated with the preferential orientation of molecular dipoles at the membrane-water interface. We discovered that the overestimation of the dipole potential by the TIPS3P water model can be effectively reduced by the TIP4P water model. On the one hand, the TIP4P water model decreases the negative contribution of lipid to the dipole potential through influencing the orientation of lipid headgroups. On the other hand, the TIP4P water model reduces the positive contribution of water to the dipole potential by increasing the preference of H-down orientation (the water dipole orients toward the bilayer center). Interestingly, the TIP4P water model affects the orientation of interfacial water molecules more obviously than that of lipid headgroups, leading to the decrease in the dipole potential. Furthermore, the MD results revealed that the water close to the positively charged choline (namely, N-associated water) prefers the H-down orientation while the water around the negatively charged phosphate (namely, P-associated water) favors the H-up orientation, in support of recent experimental and MD studies. However, interfacial water molecules are more strongly influenced by the phosphate groups than by the choline groups, resulting in the net H-up orientation (the water dipole orients toward the bilayer center) in the region of lipid headgroups. In addition, it is intriguing that the preference of H-up orientation decreases when water molecules penetrate more deeply into the lipid bilayer. This is attributed to the counteracting effect of lipid carbonyl groups, and the effect varies with the lipid chains (oleoyl and palmitoyl chains), suggesting the important role of lipid carbonyl groups.

Predicting Partition Coefficients of Neutral and Charged Solutes in the Mixed SLES-Fatty Acid Micellar System

Sodium laureth sulfate (SLES) and fatty acids are common ingredients in many cosmetic products. Understanding how neutral and charged fatty acid compounds partition between micellar and water phases is crucial to achieve the optimal design of the product formulation. In this paper, we first study the formation of mixed SLES and fatty acid micelles using molecular dynamics (MD) simulations. Micelle/water partition coefficients of neutral and charged fatty acids are then calculated using COSMOmic as well as a MD approach based on the potential of mean force (PMF) calculations performed using umbrella sampling (US). The combined US/PMF approach was performed with both the additive, non-polarizable CHARMM general force field (CGenFF) and the classical Drude polarizable force field. The partition coefficients for the neutral solutes are shown to be accurately calculated with the COSMOmic and additive CGenFF US/PMF approaches, while only the US/PMF approach with the Drude polarizable force field accurately calculated the experimental partition coefficient of the charged solute. These results indicate the utility of the Drude polarizable force field as a tool for the rational development of mixed micelles.

First-Principles Many-Body Nonadditive Polarization Energies from Monomer and Dimer Calculations Only: A Case Study on Water

The many-body polarization energy is the major source of nonadditivity in strongly polar systems such as water. This nonadditivity is often considerable and must be included, if only in an average manner, to correctly describe the physical properties of the system. Models for the polarization energy are usually parametrized using experimental data, or theoretical estimates of the many-body effects. Here we show how many-body polarization models can be developed for water complexes using data for the monomer and dimer only using ideas recently developed in the field of intermolecular perturbation theory and state-of-the-art approaches for calculating distributed molecular properties based on the iterated stockholder atoms (ISA) algorithm. We show how these models can be calculated, and we validate their accuracy in describing the many-body nonadditive energies of a range of water clusters. We further investigate their sensitivity to the details of the polarization damping models used. We show how our very best polarization models yield many-body energies that agree with those computed with coupled-cluster methods, but at a fraction of the computational cost.

Bottom-up derived flexible water model with dipole and quadrupole moments for coarse-grained molecular simulations

A bottom-up CG water model is developed to capture the electrostatic multipoles, structural correlation and thermodynamics of water.

Accurate Biomolecular Simulations Account for Electronic Polarization

In this perspective, we discuss where and how accounting for electronic many-body polarization affects the accuracy of classical molecular dynamics simulations of biomolecules. While the effects of electronic polarization are highly pronounced for molecules with an opposite total charge, they are also non-negligible for interactions with overall neutral molecules. For instance, neglecting these effects in important biomolecules like amino acids and phospholipids affects the structure of proteins and membranes having a large impact on interpreting experimental data as well as building coarse grained models. With the combined advances in theory, algorithms and computational power it is currently realistic to perform simulations with explicit polarizable dipoles on systems with relevant sizes and complexity. Alternatively, the effects of electronic polarization can also be included at zero additional computational cost compared to standard fixed-charge force fields using the electronic continuum correction, as was recently demonstrated for several classes of biomolecules.

Effect of Cholesterol and 6-Ketocholestanol on Membrane Dipole Potential and Sterol Flip-Flop Motion in Bilayer Membranes

A variety of experimental and theoretical approaches have been employed to investigate the sterol flip-flop motion in lipid bilayer membranes. However, the sterol effect on the dipole potential of lipid bilayer membranes is less well studied and the influence of dipole potential on sterol flip-flop motion in lipid bilayer membranes is less well understood. In our previous works, we have demonstrated the performance of our coarse-grained (CG) model in the computation of the dipole potential. In this work, five 30 μs CG simulations of dimyristoylphosphatidylcholine (DMPC) bilayers were carried out at different sterol concentrations (in a range from 10 to 50% mole fraction). Then, a comparison was made between the effects of cholesterol (CHOL) and 6-ketocholestanol (6-KC) on the dipole potential of DMPC lipid bilayers as well as the sterol flip-flop motion. Our CG simulations show that the membrane dipole potential is impacted more significantly by 6-KC than by CHOL. This finding is consistent with recent experimental studies. Meanwhile, our work suggests that the sterol-sterol interactions (in particular, electrostatic interactions) should be critical to the formation of sterol-sterol clusters, which would hinder the sterol flip-flop motion inside the lipid bilayers. This is in support of the recent experimental study on the sterol transportation in lipid bilayer membranes.

Force Fields for Small Molecules

Molecular dynamics (MD) simulations have been widely applied to computer-aided drug design (CADD). While MD has been used in a variety of applications such as free energy perturbation and long-time simulations, the accuracy of the results from those methods depends strongly on the force field used. Force fields for small molecules are crucial, as they not only serve as building blocks for developing force fields for larger biomolecules but also act as model compounds that will be transferred to ligands used in CADD. Currently, a wide range of small molecule force fields based on additive or nonpolarizable models have been developed. While these nonpolarizable force fields can produce reasonable estimations of physical properties and have shown success in a variety of systems, there is still room for improvements due to inherent limitations in these models including the lack of an electronic polarization response. For this reason, incorporating polarization effects into the energy function underlying a force field is believed to be an important step forward, giving rise to the development of polarizable force fields. Recent simulations of biological systems have indicated that polarizable force fields are able to provide a better physical representation of intermolecular interactions and, in many cases, better agreement with experimental properties than nonpolarizable, additive force fields. Therefore, this chapter focuses on the development of small molecule force fields with emphasis on polarizable models. It begins with a brief introduction on the importance of small molecule force fields and their evolution from additive to polarizable force fields. Emphasis is placed on the additive CHARMM General Force Field and the polarizable force field based on the classical Drude oscillator. The theory for the Drude polarizable force field and results for small molecules are presented showing their improvements over the additive model. The potential importance of polarization for their application in a wide range of biological systems including CADD is then discussed.

Simulations of lipid bilayers using the CHARMM36 force field with the TIP3P-FB and TIP4P-FB water models

The CHARMM36 force field for lipids is widely used in simulations of lipid bilayers. The CHARMM family of force fields were developed for use with the mTIP3P water model. This water model has an anomalously high dielectric constant and low viscosity, which limits its accuracy in the calculation of quantities like permeability coefficients. The TIP3P-FB and TIP4P-FB water models are more accurate in terms of the dielectric constant and transport properties, which could allow more accurate simulations of systems containing water and lipids. To test whether the CHARMM36 lipid force field is compatible with the TIP3P-FB and TIP4P-FB water models, we have performed simulations of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine bilayers. The calculated headgroup area, compressibility, order parameters, and X-ray form factors are in good agreement with the experimental values, indicating that these improved water models can be used with the CHARMM36 lipid force field without modification when calculating membrane physical properties. The water permeability predicted by these models is significantly different; the mTIP3P-model diffusion in solution and at the lipid-water interface is anomalously fast due to the spuriously low viscosity of mTIP3P-model water, but the potential of mean force of permeation is higher for the TIP3P-FB and TIP4P-FB models due to their high excess chemical potentials. As a result, the rates of water permeation calculated the FB water models are slower than the experimental value by a factor of 15-17, while simulations with the mTIP3P model only underestimate the water permeability by a factor of 3.

Link between Membrane Composition and Permeability to Drugs

Prediction of membrane permeability to small molecules represents an important aspect of drug discovery. First-principles calculations of this quantity require an accurate description of both the thermodynamics and kinetics that underlie translocation of the permeant across the lipid bilayer. In this contribution, the membrane permeability to three drugs, or drug-like molecules, namely, 9-anthroic acid (ANA), 2',3'-dideoxyadenosine (DDA), and hydrocortisone (HYL), are estimated in a pure 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) and in a POPC:cholesterol (2:1) mixture. On the basis of independent 2-5-μs free-energy calculations combined with a time-fractional Smoluchowski determination of the diffusivity, the estimated membrane permeabilities to these chemically diverse permeants fall within an order of magnitude from the experimental values obtained in egg-lecithin bilayers, with the exception of HYL in pure POPC. This exception is particularly interesting because the calculated permeability of the sterol-rich bilayer to HYL, in close agreement with the experimental value, is about 600 times lower than that of the pure POPC bilayer to HYL. In contrast, the permeabilities to ANA and DDA differ by less than a factor of 10 between the pure POPC and POPC:cholesterol bilayers. The unusual behavior of HYL, a large, amphiphilic compound, may be linked with the longer range perturbation of the lipid bilayer it induces, compared to ANA and DDA, suggestive of a possibly different translocation mechanism. We find that the tendency of lower permeabilities of the POPC:cholesterol bilayer relative to those of the pure POPC one is a consequence of increased free-energy barriers. Beyond reporting accurate estimates of the membrane permeability, the present contribution also demonstrates that rigorous free-energy calculations and a fractional-diffusion model are key in revealing the molecular phenomena linking the composition of a membrane to its permeability to drugs.

In silico assessment of the conduction mechanism of the Ryanodine Receptor 1 reveals previously unknown exit pathways

The ryanodine receptor 1 is a large calcium ion channel found in mammalian skeletal muscle. The ion channel gained a lot of attention recently, after multiple independent authors published near-atomic cryo electron microscopy data. Taking advantage of the unprecedented quality of structural data, we performed molecular dynamics simulations on the entire ion channel as well as on a reduced model. We calculated potentials of mean force for Ba2+, Ca2+, Mg2+, K+, Na+ and Cl- ions using umbrella sampling to identify the key residues involved in ion permeation. We found two main binding sites for the cations, whereas the channel is strongly repulsive for chloride ions. Furthermore, the data is consistent with the model that the receptor achieves its ion selectivity by over-affinity for divalent cations in a calcium-block-like fashion. We reproduced the experimental conductance for potassium ions in permeation simulations with applied voltage. The analysis of the permeation paths shows that ions exit the pore via multiple pathways, which we suggest to be related to the experimental observation of different subconducting states.

Efficient minimization of multipole electrostatic potentials in torsion space

The development of models of macromolecular electrostatics capable of delivering improved fidelity to quantum mechanical calculations is an active field of research in computational chemistry. Most molecular force field development takes place in the context of models with full Cartesian coordinate degrees of freedom. Nevertheless, a number of macromolecular modeling programs use a reduced set of conformational variables limited to rotatable bonds. Efficient algorithms for minimizing the energies of macromolecular systems with torsional degrees of freedom have been developed with the assumption that all atom-atom interaction potentials are isotropic. We describe novel modifications to address the anisotropy of higher order multipole terms while retaining the efficiency of these approaches. In addition, we present a treatment for obtaining derivatives of atom-centered tensors with respect to torsional degrees of freedom. We apply these results to enable minimization of the Amoeba multipole electrostatics potential in a system with torsional degrees of freedom, and validate the correctness of the gradients by comparison to finite difference approximations. In the interest of enabling a complete model of electrostatics with implicit treatment of solvent-mediated effects, we also derive expressions for the derivative of solvent accessible surface area with respect to torsional degrees of freedom.

Drude Polarizable Force Field for Molecular Dynamics Simulations of Saturated and Unsaturated Zwitterionic Lipids

Additive force fields are designed to account for induced electronic polarization in a mean-field average way, using effective empirical fixed charges. The limitation of this approximation is cause for serious concerns, particularly in the case of lipid membranes, where the molecular environment undergoes dramatic variations over microscopic length scales. A polarizable force field based on the classical Drude oscillator offers a practical and computationally efficient framework for an improved representation of electrostatic interactions in molecular simulations. Building on the first-generation Drude polarizable force field for the dipalmitoylphosphatidylcholine 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) molecule, the present effort was undertaken to improve this initial model and expand the force field to a wider range of phospholipid molecules. New lipids parametrized include dimyristoylphosphatidylcholine (DMPC), dilauroylphosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), dipalmitoylphosphatidylethanolamine (DPPE), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE), and 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE). The iterative optimization protocol employed in this effort led to lipid models that achieve a good balance between reproducing quantum mechanical data on model compound representative of phospholipids and reproducing a range of experimental condensed phase properties of bilayers. A parametrization strategy based on a restrained ensemble-maximum entropy methodology was used to help accurately match the experimental NMR order parameters in the polar headgroup region. All the parameters were developed to be compatible with the remainder of the Drude polarizable force field, which includes water, ions, proteins, DNA, and selected carbohydrates.

The dipole potential correlates with lipid raft markers in the plasma membrane of living cells

The dipole potential generating an electric field much stronger than any other type of membrane potential influences a wide array of phenomena, ranging from passive permeation to voltage-dependent conformational changes of membrane proteins. It is generated by the ordered orientation of lipid carbonyl and membrane-attached water dipole moments. Theoretical considerations and indirect experimental evidence obtained in model membranes suggest that the dipole potential is larger in liquid-ordered domains believed to correspond to lipid rafts in cell membranes. Using three different dipole potential-sensitive fluorophores and four different labeling approaches of raft and nonraft domains, we showed that the dipole potential is indeed stronger in lipid rafts than in the rest of the membrane. The magnitude of this difference is similar to that observed between the dipole potential in control and sphingolipid-enriched cells characteristic of Gaucher's disease. The results established that the heterogeneity of the dipole potential in living cell membranes is correlated with lipid rafts and imply that alterations in the lipid composition of the cell membrane in human diseases can lead to substantial changes in the dipole potential.

An averaged polarizable potential for multiscale modeling in phospholipid membranes

A set of average atom-centered charges and polarizabilities has been developed for three types of phospholipids for use in polarizable embedding calculations. The lipids investigated are 1,2-dimyristoyl-sn-glycero-3-phosphocholine, 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine, and 1-palmitoyl-2-oleoyl-sn-glycerol-3-phospho-L-serine given their common use both in experimental and computational studies. The charges, and to a lesser extent the polarizabilities, are found to depend strongly on the molecular conformation of the lipids. Furthermore, the importance of explicit polarization is underlined for the description of larger assemblies of lipids, that is, membranes. In conclusion, we find that specially developed polarizable parameters are needed for embedding calculations in membranes, while common non-polarizable point-charge force fields usually perform well enough for structural and dynamical studies. © 2017 Wiley Periodicals, Inc.

Molecular-level insight into the binding of arginine to a zwitterionic Langmuir monolayer

Arginine molecules bind to a DPPC monolayer, altering the interfacial electrostatic potential and the lateral mobility of the lipids, while having little effect on the compression isotherm of the monolayer.

An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications

Molecular mechanics force fields that explicitly account for induced polarization represent the next generation of physical models for molecular dynamics simulations. Several methods exist for modeling induced polarization, and here we review the classical Drude oscillator model, in which electronic degrees of freedom are modeled by charged particles attached to the nuclei of their core atoms by harmonic springs. We describe the latest developments in Drude force field parametrization and application, primarily in the last 15 years. Emphasis is placed on the Drude-2013 polarizable force field for proteins, DNA, lipids, and carbohydrates. We discuss its parametrization protocol, development history, and recent simulations of biologically interesting systems, highlighting specific studies in which induced polarization plays a critical role in reproducing experimental observables and understanding physical behavior. As the Drude oscillator model is computationally tractable and available in a wide range of simulation packages, it is anticipated that use of these more complex physical models will lead to new and important discoveries of the physical forces driving a range of chemical and biological phenomena.

Free Energy Calculations for the Peripheral Binding of Proteins/Peptides to an Anionic Membrane. 1. Implicit Membrane Models

The binding of peptides and proteins to the surface of complex lipid membranes is important in many biological processes such as cell signaling and membrane remodeling. Computational studies can aid experiments by identifying physical interactions and structural motifs that determine the binding affinity and specificity. However, previous studies focused on either qualitative behaviors of protein/membrane interactions or the binding affinity of small peptides. Motivated by this observation, we set out to develop computational protocols for bimolecular binding to charged membranes that are applicable to both peptides and large proteins. In this work, we explore a method based on an implicit membrane/solvent model (generalized Born with a simple switching in combination with the Gouy-Chapman-Stern model for a charged interface), which we expect to lead to useful results when the binding does not implicate significant membrane deformation and local demixing of lipids. We show that the binding free energy can be efficiently computed following a thermodynamic cycle similar to protein-ligand binding calculations, especially when a Bennett acceptance ratio based protocol is used to consider both the membrane bound and solution conformational ensembles. Test calculations on a series of peptides show that our computational approach leads to binding affinities in encouraging agreement with experimental data, including for the challenging example of the bringing of flexible MARCKS-ED peptides to membranes. The calculations highlight that for a membrane with a significant fraction of anionic lipids, it is essential to include the effect of ion adsorption using the Stern model, which significantly modifies the effective surface charge. This implicit membrane model based computational protocol helps lay the groundwork for more systematic analysis of protein/peptide binding to membranes of complex shape and composition.

Molecular simulation of nonfacilitated membrane permeation

This is a review. Non-electrolytic compounds typically cross cell membranes by passive diffusion. The rate of permeation is dependent on the chemical properties of the solute and the composition of the lipid bilayer membrane. Predicting the permeability coefficient of a solute is important in pharmaceutical chemistry and toxicology. Molecular simulation has proven to be a valuable tool for modeling permeation of solutes through a lipid bilayer. In particular, the solubility-diffusion model has allowed for the quantitative calculation of permeability coefficients. The underlying theory and computational methods used to calculate membrane permeability are reviewed. We also discuss applications of these methods to examine the permeability of solutes and the effect of membrane composition on permeability. The application of coarse grain and polarizable models is discussed. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

CHARMM additive and polarizable force fields for biophysics and computer-aided drug design

BACKGROUND

Molecular Mechanics (MM) is the method of choice for computational studies of biomolecular systems owing to its modest computational cost, which makes it possible to routinely perform molecular dynamics (MD) simulations on chemical systems of biophysical and biomedical relevance.

SCOPE OF REVIEW

As one of the main factors limiting the accuracy of MD results is the empirical force field used, the present paper offers a review of recent developments in the CHARMM additive force field, one of the most popular biomolecular force fields. Additionally, we present a detailed discussion of the CHARMM Drude polarizable force field, anticipating a growth in the importance and utilization of polarizable force fields in the near future. Throughout the discussion emphasis is placed on the force fields' parametrization philosophy and methodology.

MAJOR CONCLUSIONS

Recent improvements in the CHARMM additive force field are mostly related to newly found weaknesses in the previous generation of additive force fields. Beyond the additive approximation is the newly available CHARMM Drude polarizable force field, which allows for MD simulations of up to 1μs on proteins, DNA, lipids and carbohydrates.

GENERAL SIGNIFICANCE

Addressing the limitations ensures the reliability of the new CHARMM36 additive force field for the types of calculations that are presently coming into routine computational reach while the availability of the Drude polarizable force fields offers an inherently more accurate model of the underlying physical forces driving macromolecular structures and dynamics. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".

Representation of Ion-Protein Interactions Using the Drude Polarizable Force-Field

Small metal ions play critical roles in numerous biological processes. Of particular interest is how metalloenzymes are allosterically regulated by the binding of specific ions. Understanding how ion binding affects these biological processes requires atomic models that accurately treat the microscopic interactions with the protein ligands. Theoretical approaches at different levels of sophistication can contribute to a deeper understanding of these systems, although computational models must strike a balance between accuracy and efficiency in order to enable long molecular dynamics simulations. In this study, we present a systematic effort to optimize the parameters of a polarizable force field based on classical Drude oscillators to accurately represent the interactions between ions (K⁺, Na⁺, Ca²⁺, and Cl^–) and coordinating amino-acid residues for a set of 30 biologically important proteins. By combining ab initio calculations and experimental thermodynamic data, we derive a polarizable force field that is consistent with a wide range of properties, including the geometries and interaction energies of gas-phase ion/protein-like model compound clusters, and the experimental solvation free-energies of the cations in liquids. The resulting models display significant improvements relative to the fixed-atomic-charge additive CHARMM C36 force field, particularly in their ability to reproduce the many-body electrostatic nonadditivity effects estimated from ab initio calculations. The analysis clarifies the fundamental limitations of the pairwise additivity assumption inherent in classical fixed-charge force fields, and shows its dramatic failures in the case of Ca²⁺ binding sites. These optimized polarizable models, amenable to computationally efficient large-scale MD simulations, set a firm foundation and offer a powerful avenue to study the roles of the ions in soluble and membrane transport proteins.

Current status of protein force fields for molecular dynamics simulations

The current status of classical force fields for proteins is reviewed. These include additive force fields as well as the latest developments in the Drude and AMOEBA polarizable force fields. Parametrization strategies developed specifically for the Drude force field are described and compared with the additive CHARMM36 force field. Results from molecular simulations of proteins and small peptides are summarized to illustrate the performance of the Drude and AMOEBA force fields.

A systematic molecular dynamics simulation study of temperature dependent bilayer structural properties

Although lipid force fields (FFs) used in molecular dynamics (MD) simulations have proved to be accurate, there has not been a systematic study on their accuracy over a range of temperatures. Motivated by the X-ray and neutron scattering measurements of common phosphatidylcholine (PC) bilayers (Kučerka et al. BBA. 1808: 2761, 2011), the CHARMM36 (C36) FF accuracy is tested in this work with MD simulations of six common PC lipid bilayers over a wide range of temperatures. The calculated scattering form factors and deuterium order parameters from the C36 MD simulations agree well with the X-ray, neutron, and NMR experimental data. There is excellent agreement between MD simulations and experimental estimates for the surface area per lipid, bilayer thickness (DB), hydrophobic thickness (DC), and lipid volume (VL). The only minor discrepancy between simulation and experiment is a measure of (DB-DHH)/2 where DHH is the distance between the maxima in the electron density profile along the bilayer normal. Additional MD simulations with pure water and heptane over a range of temperatures provide explanations of possible reasons causing the minor deviation. Overall, the C36 FF is accurate for use with liquid crystalline PC bilayers of varying chain types and over biologically relevant temperatures.

Force Field for Peptides and Proteins based on the Classical Drude Oscillator

Presented is a polarizable force field based on a classical Drude oscillator framework, currently implemented in the programs CHARMM and NAMD, for modeling and molecular dynamics (MD) simulation studies of peptides and proteins. Building upon parameters for model compounds representative of the functional groups in proteins, the development of the force field focused on the optimization of the parameters for the polypeptide backbone and the connectivity between the backbone and side chains. Optimization of the backbone electrostatic parameters targeted quantum mechanical conformational energies, interactions with water, molecular dipole moments and polarizabilities and experimental condensed phase data for short polypeptides such as (Ala)₅. Additional optimization of the backbone φ, ψ conformational preferences included adjustments of the tabulated two-dimensional spline function through the CMAP term. Validation of the model included simulations of a collection of peptides and proteins. This 1^st generation polarizable model is shown to maintain the folded state of the studied systems on the 100 ns timescale in explicit solvent MD simulations. The Drude model typically yields larger RMS differences as compared to the additive CHARMM36 force field (C36) and shows additional flexibility as compared to the additive model. Comparison with NMR chemical shift data shows a small degradation of the polarizable model with respect to the additive, though the level of agreement may be considered satisfactory, while for residues shown to have significantly underestimated S² order parameters in the additive model, improvements are calculated with the polarizable model. Analysis of dipole moments associated with the peptide backbone and tryptophan side chains show the Drude model to have significantly larger values than those present in C36, with the dipole moments of the peptide backbone enhanced to a greater extent in sheets versus helices and the dipoles of individual moieties observed to undergo significant variations during the MD simulations. Although there are still some limitations, the presented model, termed Drude-2013, is anticipated to yield a molecular picture of peptide and protein structure and function that will be of increased physical validity and internal consistency in a computationally accessible fashion.

The Polarizable Atomic Multipole-based AMOEBA Force Field for Proteins

Development of the AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Simulation) force field for proteins is presented. The current version (AMOEBA-2013) utilizes permanent electrostatic multipole moments through the quadrupole at each atom, and explicitly treats polarization effects in various chemical and physical environments. The atomic multipole electrostatic parameters for each amino acid residue type are derived from high-level gas phase quantum mechanical calculations via a consistent and extensible protocol. Molecular polarizability is modeled via a Thole-style damped interactive induction model based upon distributed atomic polarizabilities. Inter- and intramolecular polarization is treated in a consistent fashion via the Thole model. The intramolecular polarization model ensures transferability of electrostatic parameters among different conformations, as demonstrated by the agreement between QM and AMOEBA electrostatic potentials, and dipole moments of dipeptides. The backbone and side chain torsional parameters were determined by comparing to gas-phase QM (RI-TRIM MP2/CBS) conformational energies of dipeptides and to statistical distributions from the Protein Data Bank. Molecular dynamics simulations are reported for short peptides in explicit water to examine their conformational properties in solution. Overall the calculated conformational free energies and J-coupling constants are consistent with PDB statistics and experimental NMR results, respectively. In addition, the experimental crystal structures of a number of proteins are well maintained during molecular dynamics (MD) simulation. While further calculations are necessary to fully validate the force field, initial results suggest the AMOEBA polarizable multipole force field is able to describe the structure and energetics of peptides and proteins, in both gas-phase and solution environments.

A polarizable force field of dipalmitoylphosphatidylcholine based on the classical Drude model for molecular dynamics simulations of lipids

A polarizable force field of saturated phosphatidylcholine-containing lipids based on the classical Drude oscillator model is optimized and used in molecular dynamics simulations of bilayer and monolayer membranes. The hierarchical parametrization strategy involves the optimization of parameters for small molecules representative of lipid functional groups, followed by their application in larger model compounds and full lipids. The polar headgroup is based on molecular ions tetramethyl ammonium and dimethyl phosphate, the esterified glycerol backbone is based on methyl acetate, and the aliphatic lipid hydrocarbon tails are based on linear alkanes. Parameters, optimized to best represent a collection of gas and liquid properties for these compounds, are assembled into a complete model of dipalmitoylphosphatidylcholine (DPPC) lipids that is tested against the experimental properties of bilayer and monolayer membranes. The polarizable model yields average structural properties that are in broad accord with experimental data. The area per lipid of the model is 60 Å(2), slightly smaller than the experimental value of 63 Å(2). The order parameters from nuclear magnetic resonance deuterium quadrupolar splitting measures, the electron density profile, and the monolayer dipole potential are in reasonable agreement with experimental data, and with the nonpolarizable CHARMM C36 lipid force field.

Role of the Membrane Dipole Potential for Proton Transport in Gramicidin A Embedded in a DMPC Bilayer

The membrane potential at the water/phospholipid interfaces may play a key role for proton conduction of gramicidin A (gA). Here we address this issue by Density Functional Theory-based molecular dynamics and metadynamics simulations. The calculations, performed on gA embedded in a solvated 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) model membrane environment (about 2,000 atoms), indicate that (i) the membrane dipole potential rises at the channel mouth by ∼0.4 V. A similar value has been measured for gA embedded in a DMPC monolayer; (ii) the calculated free energy barrier is located at the channel entrance, consistent with experiments comparing gA proton conduction in different bilayers. The electronic structures of the proton ligands (water molecules and peptide units) are similar to those in the bulk solvent. Based on these results, we suggest an important role of the membrane dipole potential for the free energy barrier of proton permeation of gA. This may provide a rationale for the large increase in the rate of proton conduction under application of a transmembrane voltage, as observed experimentally. Our calculations might suggest also a role for proton desolvation for the permeation process. This role has already emerged from EVB calculations on gA embedded in a model membrane.