|Home | About | Journals | Submit | Contact Us | Français|
The development of the CHARMM additive all-atom lipid force field (FF) is traced from the early 1990’s to the most recent version (C36) published in 2010. Though simulations with early versions yielded useful results, they failed to reproduce two important quantities: a zero surface tension at the experimental bilayer surface area, and the signature splitting of the deuterium order parameters in the glycerol and upper chain carbons. Systematic optimization of parameters based on high level quantum mechanical data and free energy simulations have resolved these issues, and bilayers with a wide range of lipids can be simulated in tensionless ensembles using C36. Issues associated with other all-atom lipid FFs, success and limitations in the C36 FF and ongoing developments are also discussed.
Lipids are amphiphilic molecules with polar heads and hydrocarbon tails. This attribute allows them to form liquid crystals that are both thermotropic and lyotropic; i.e., their phase behavior is modulated by temperature and solvent. Biological lipids are also quite variable. Their head groups can be highly charged, zwitterionic, or absent charged groups; the tails are different lengths, have varying degrees of unsaturation, and some are even branched. These features of lipids give cell membranes remarkable complexity and sensitivity to their surroundings, and the ability to carry out a wide range of functions such as: act as permeability barriers between the inside and outside of cells and compartments of cells; provide specialized environments for proteins and modulate enzymatic reactions; participate in cell signaling pathways; and undergo huge deformations such as pore formation and fusion, vesiculation and tabulation.
The functions noted above are being increasingly studied by classical molecular dynamics (MD) and Monte Carlo (MC) computer simulation carried out using models of varying detail. The models range from all-atom (or heavy atom), to coarse grained (where 3 or 4 heavy atoms are combined to form a single particle), to mesoscopic (where an entire lipid is represented by one or two particles), with the treatment of solvent ranging from explicit to continuum representations. Structural and dynamical details of lipid/lipid, lipid/protein, and lipid/solvent interaction are presently the realm of all-atom simulations. Slower processes such as long-wave length undulation, pore formation, and self-assembly are examined by coarse grained models. Very long length and time scale processes, sometimes termed “membrane remodeling,” are accessed by mesoscopic models. At the core of each model is a force field (FF), which specifies how the particles interact. The aforementioned sensitivity of membranes to temperature, solvent, and the detailed structure of component lipids has made development of lipid FFs particularly challenging. Additionally, because coarse grained and mesoscopic models are often parametrized at least in part using the results of all-atom models, errors can propagate between FFs.
The usefulness of a simulation is determined not only by the quality of the force field, but also by understanding its limitations. The focus of this review is on the CHARMM all-atom FF1,2. Following an overview and a brief history, issues associated with other commonly used FFs are discussed and the successes and limitations of the most recent CHARMM FF, C36,3 are summarized. This review closes with a discussion of ongoing lipid FF development efforts, including the explicit treatment of electronic polarization.
It is a fair generalization (at least at present) to assert that all-atom simulations primarily explore the thermotropic character of lipids, as opposed to the lyotropic; i.e., a particular phase is specified in the simulation (e.g., bilayer) and transitions to another solvent induced phase (e.g. micelle) are not the subject of the simulation. Lipid bilayers in cells are primarily in the liquid-crystalline (Lα, or fluid) phase; the analogous phase for lipid monolayers is the liquid-expanded phase. The individual lipids in both phases are ordered with respect the surface normal, but are highly dynamic and otherwise do not occupy well defined positions.
To investigate the physical properties of lipids a number of experimental methods have been applied including NMR, X-ray and neutron diffraction, and fluorescence techniques. Such experiments supply detailed information about the dynamics or distribution of various chemical moieties in lipids as well as the surrounding solvent as a function of the mono- or bilayer depth. However, the inherent disorder along the plane typically leads to some model dependence when extracting atomic detail from experiment. All-atom molecular dynamics simulations generate conformations of the system in atomic detail and are a natural complement to experiments. The results of simulations have been used to test the models assumed in analysis of experiment, or to develop new models to refine the interpretation of experiment. Conversely, the experiments provide an invaluable source of data to validate the simulation methodology, including the force fields. While early simulations focused on small patches of pure bilayers of simple lipids, particularity 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC, Figure 1), the last decade has seen many studies of more complex systems including imbedded or surface peptides and proteins.
Development of a FF proceeds by: (i) specification of a potential energy function, V(), where is a generalized position vector; (ii) fitting of the parameters of V() based on experimental or quantum mechanical data, including simulations of both simple test systems and lipid mono- and bilayers; validation is subsequently performed by simulations of large systems not included in the training set and comparison with experiment. Over the last 20 years a number of lipid force fields have been developed. Methodological and computational improvements have allowed for those force fields to be more rigorously tested, and many have undergone significant improvements, as detailed below.
The CHARMM additive lipid force field has been designed to be compatible with CHARMM additive force fields1 for the other biopolymers and drug-like compounds. The form is
The bond and angles terms are harmonic, with force constants Kb and Kθ, and equilibrium values b0 and θ0. The dihedral potential is a sum of sinusoids with force constant Kϕ, j, multiplicity n, offset δ, and j can range from 1 to 6. A Lennard-Jones (LJ) potential models the van der Waals interactions, where εij is the potential energy minimum between two particles separated by rij, and Rmin,ij is the position of this minimum. Lastly, qi and qj are the partial atomic charges for the Columbic term; the dielectric constant εD equals 1 in explicit solvent simulations. These partial atomic charges are fixed such that Eq. (1) describes an additive FF. Polarizable FFs are discussed below. While the form of the energy function in Eq. (1) is similar to that of other commonly used biomolecular FFs (eg. AMBER,4 OPLS,5 GROMOS6), the FF itself is unique because of the approach to parameter optimization. Accordingly, it is highly recommended that FFs of different origins not be mixed due to inconsistencies in the treatment of the nonbond terms; e.g., the CHARMM lipid FF should only be used with the remainder of the CHARMM additive FF.
The following are some underappreciated points concerning the CHARMM and most other all- or heavy-atom FFs:
The original CHARMM lipid FF, C22, by Schlenkrich et al.15 was developed in the early 1990’s. Parameter optimization was based on small molecules for which a range of target data exists; these included butane, methylacetate (MAS), dimethylphosphate (DMP), and tetramethylammonium (TMA). This data included experimental pure solvent heats of vaporization and densities, and ab initio QM results for geometries, vibrational spectra, conformational energies, and interactions with water. The QM calculations were limited to the Hartree-Fock level by computational considerations, though the use of HF/6-31G(d) interactions with water still remains the standard for the CHARMM additive FF. The small molecule parameters were then combined to create three phospholipids, lauroylpropanediol-phosphorycholine (LPPC), 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), and 1,2-dilauroyl-sn-glycero-1-phosorylethanolamine (DLPE), and the resultant FF was validated using crystal minimizations and MD simulations. The dynamics simulations were limited to 100 ps in the NVE ensemble, so a robust test of the ability of the nonbond parameters to reproduce the unit cell parameters was not performed. A 190 ps NVE simulation of a DPPC bilayer, which was not part of the optimization process, led to an overall stable structure and successfully reproduced the “fast” NMR T1 relaxation times.16 This relatively short simulation required 6 months on an IBM3090, an advanced machine for the time (the equivalent simulation would presently take only several hours on a typical laboratory cluster). Subsequent longer NPAT simulations (600–800 ps) at a range of surface areas17 demonstrated that the deuterium order parameters (|SCD|) in the bilayer interior agreed well with experiment18 at 62.9 Å2/lipid, an area close to the value 64 Å2/lipid determined by x-ray diffraction at the same temperature.19 The model was also extended to 1,2-dioleoyl-sn-phosphatidylcholine (DOPC),20 and an NPAT trajectory yielded excellent agreement with position distributions from neutron diffraction data.
While the successes of these and other early studies established the utility of bilayer simulations,21 deficiencies in the FF were already apparent. In particular, the 800 ps DPPC bilayer simulation at 62.9 Å2/lipid yielded a surface tension of 39 dyn/cm/side. Sound statistical mechanical arguments indicated the that surface tension calculated from a simulation should be 0 for a perfectly flat bilayer at the equilibrium surface area,22–24 and only a few dyn/cm if undulations were taken into account.25 (Estimates based on the standard summation over wave vectors26 by Feller and Pastor27 implied that the suppression of undulations should yield sizeable surface tensions for simulation sized systems. These estimates were recently found to be too high because the integral must be taken over square, not circular, domains (Lerner and Pastor, unpublished).) The high surface tension of DPPC bilayers simulated with C22 implies that if the same bilayer were simulated at NPT (the equivalent of γ=0), the surface area would contract and the chains would become gel-like. The order parameters near the water interface for systems simulated with the C22 FF were also in error (Table 1). In particular, the inequivalence of the chains and precise orientation of the head group defined by the signature splitting at the acyl chain C2 and glycerol G1 positions was not reproduced. Accordingly, subsequent parametrizations aimed at improving the order parameter profile and lowering the surface tension. However, because a zero surface tension can be obtained in many ways, the approach was to refine the FF parameters against other data, and then evaluate γ. The area collapse could be avoided by either simulating at constant area, or with an applied surface tension set to yield the experimental surface area, so there was a reasonable short-term way around the problem.
C27, published in 2000,30 was the first systematic refinement of C22. It involved improvements in the alkane LJ and torsional parameters, and to the partial atomic charges and torsional parameters of the phosphate moiety. QM surfaces for butane were generated with explicit treatment of electron correlation. The FF was explicitly parameterized with long range electrostatics, and was validated in an 11 ns simulation of a DPPC bilayer. C27 was subsequently used for a range of studies many of which pushed the MD time scale into the 100 ns regimen. At these longer time scales a small but systematic overestimate of the chain order parameters became evident; this emphasizes the importance of obtaining proper sampling during an MD simulation to overcome inherent bias associated with the starting configuration of the system. Subsequently, additional optimization of the alkane torsion parameters ensued including high level QM calculations on larger model compounds, such as hexane and heptane. The resulting set, C27r,12 lead to excellent agreement with the aliphatic chain order parameters, though issues with splitting of the glycerol and C2 order parameters persisted (Table 1). Perhaps more importantly, the C27 and C27r parameter sets both yielded surface tensions of approximately 20 dyn/cm/side for DPPC bilayers at the experimental surface area; i.e., it was still inadvisable to carry out simulations at NPT.
At this stage, three problems with the force field were considered: the nonbond parameters of the head groups (those associated with DMP and TMA); the torsion parameters for the head group and esterified glycerol linker region; the nonbond parameters associated with the ester moieties. Free energy of aqueous solvation calculations indicated that the TMA and DMP parameters are sufficiently accurate. Simultaneously, extensive QM calculations of larger model compounds representative of the head group/linker region were performed, and the results used to refit selected torsion parameters. This fitting also included energetic considerations concerning conformations of the head group/linker region that would lead to splitting of the relevant order parameters. This process achieved the goal of order parameter splitting, although the positive surface tension persisted. Neutron diffraction experiments on DOPC at low hydration indicated a significant underestimate by simulation of penetration of water into the linker region, and thereby strongly implicated the ester nonbond parameters. Accordingly, reevaluation of these parameters, based on MAS, was undertaken. The final model reproduced MAS pure solvent properties while the dipole moment and the free energy of aqueous hydration were systematically overestimated (i.e., ΔGaq was too favorable). Several MAS models with these general characteristics were developed and tested in extended NPT simulation of different lipids from which a final set was selected. This parameter set, termed C36,3 yields the correct surface area/head group and order parameters (Table 1); the following section describes some other properties.
A number of other lipid FF were developed over the same time period, as reviewed in the introduction by Klauda et al.3 Of these, the GROMOS force field has seen the most extensive use, and one aspect of its development forms an interesting contrast to the CHARMM path. The simulation of sodium-decanoate/decanol/water system by Egberts and Berendsen31 published in 1988 is considered the first “modern” MD simulation of a lipid bilayer, in that it contained acyl chains, head groups, and water, and did not require additional constraints to keep the bilayer in place. These authors32 later simulated a DPPC in water using the same potential function (with additional terms for the PC headgroup). However, they simulated at NPT, in contrast to Venable et al.16 who carried out their simulations at NVE at nearly the same time. The NPT trajectory immediately revealed a problem in the potential because the bilayer contracted to a gel-like configuration (at a temperature well above the phase transition). To overcome this problem, the magnitude of all the charges were decreased by a factor of 2, leading to agreement with the fluid phase experimental data. This approach was justified by the incorrect dielectric constant of the water model, but may be considered questionable and requires additional physical justification. While the FF and the follow-up FF of Berger et al.6 enabled users of GROMOS to obtain reasonable results for DPPC with NPT bilayer simulations years before CHARMM users were able to do so, the correction was not applicable to other lipid types. In contrast, the C36 FF, in which adjustments to the FF were performed at the model compound level and then applied to full lipids, has been shown to work well with 6 different lipid types (see below), though limitations with other lipids cannot be excluded.
A notable feature of C36 is that it yields accurate head group areas at NPT for a range of lipids: DPPC, DMPC, POPC (at low and high hydration), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); the agreement is within 2% on average, with a maximum difference of 5% for POPC. Deuterium order parameters, X-ray and neutron diffraction profiles, frequency dependent 13C NMR relaxation rates, and isothermal area compressibilities are all satisfactorily reproduced. The FF also includes parameters for cholesterol and polyunsaturated hydrocarbons. This battery of lipids will allow simulations of complex membranes as well as more heterogeneous systems that include other biomolecules, such as proteins or carbohydrates.
It must be emphasized that it is difficult to determine surface areas per lipid experimentally,19 and that revisions are common. Is it therefore not advisable to unconditionally rely on the surface areas from NPT simulations using C36 or any other FF. A prudent approach is to simulate at more than one surface area either directly in the NPnAT ensemble, or indirectly in the NPnγT by changing the applied surface tension. If the phenomenon under investigation is relatively insensitive to the surface area, then the robustness of the study has been demonstrated. If it is not, extra caution is advised.
Naturally, C36 contains limitations. Given the attention paid to surface tension, the most vexing is the inconsistency of bilayers versus monolayers. The results described so far for bilayers were obtained with a 12 Å truncation of the LJ interactions. When this truncation scheme is applied to DPPC monolayers, the surface tension at 64 Å2/lipid is underestimated by 15 dyn/cm.3 A similar underestimation is obtained for the surface tension of alkane/air interfaces for this cutoff, emphasizing the role of long-range LJ interactions in interfacial systems.33 For alkanes the situation is remedied by using a method that includes long-ranged LJ terms such as IPS/DFFT (Isotropic Periodic Sum, evaluated with a discrete fast Fourier transform),33,34 and a similar result is obtained for DPPC monolayers: the calculated surface tension of 44 dyn/cm is close to the experimental value of 41 dyn/cm. However, when simulations of bilayers are carried out using IPS at 64 Å2/lipid DPPC lipid bilayer, γ=14 dyn/cm/side. This implies that the areas would contract when simulated in the NPT ensemble. Consequently, it is not possible to simulate bilayers and monolayers self-consistently with C36. Of course, it is possible to simulate monolayers with the 12 Å cutoff and examine trends, but the absolute values will be offset. At this time, the solution has been to settle on parameters that yield zero surface tension for bilayers. It is reasonable to simulate monolayers with both IPS and a 12 Å cutoff, and compare the results. Since most simulation programs do not support IPS, this approach for monolayers may not be useful to groups not using CHARMM. It should also be emphasized that isotropic long-range LJ corrections are not appropriate for the highly anisotropic systems such as lipid mono- and bilayers. Another alternative that is appropriate for anisotropic systems are corrections to the LJ forces updated periodically during an MD simulation. For example, the LJ forces can be calculated using an extended cutoff (e.g., 30 Å) every ps with the difference in forces between the extended and standard truncation distance used as a correction for the intervening simulation time.35 Such a correction could readily be implemented in most simulation packages.
Another limitation of C36 is that the calculated dipole potential for DPPC bilayers, approximately 700 mV, substantially overestimates the experimental value of 250 mV.36 Most models show similar discrepancies, and it is likely that polarizable models (next section) will be required to obtain better agreement with experiment. For example, there is a 200 mV drop between the bilayer center and the end of the acyl chain region. This arises from the overestimation of the partial atomic charges in the all-atom model, and too low dielectric constant in the aliphatic region. (United models do not show this drop because there are no partial atomic charges in the aliphatic region.)
Preliminary (unpublished) comparisons of experiment37 and simulation of POPC and 1-palmitoyl-2-oleoyl-sn-phosphatidylserine (POPS) bilayers in NaCl solution indicate that Na+ binding is too tight with C36; i.e., the zeta potential is too low. Analogous problems for other FF were recently summarized by Khandelia et al.38 Given the importance of ions in biological systems, this is a matter of concern. However, due to the additive approximation in the FF, the treatment of ions is inherently problematic. A potential correction for this problem may be in the use of atom-pair specific LJ parameters that will correct the ion-lipid interactions without impacting the remainder of the balance of the nonbond interactions.
The decision to publish a FF does not mean that testing ends. It’s a time to simulate new systems and to compare old ones with recent experiments with an eye to the next revision.
Recent neutron diffraction studies39 of 1,2-dipropyl-sn-phosphatidylcholine (C3-PC) in water have yielded detailed information on hydration of the isolated head group that provide useful targets for simulation. As another example, the lack of “supporting structure” of a bilayer makes micelles potentially good targets for testing FF when the quality of experimental data is high. Low angle X-ray scattering results of dodecylphospocholine (DPC) micelles40 appear to be such a target.
Development of the CHARMM FF was greatly aided by experimental data on surface area, isothermal compressibility, and NMR deuterium order parameters and T1 relaxation times on simple lipids. It would be very helpful to have similar experimental data on more complex lipids and mixtures for further development.
Turning to methodology, the partial atomic charges in an additive FF represent a significant compromise in that target data ranging from the gas phase to nonpolar to polar environments cannot be reproduced with a single set of fixed charges.41 This is not surprising because the charge distribution need to adjust to the environment. Hence it is likely that some of the deficiencies and internal inconsistencies present in the additive model will be remedied by adopting models that explicitly treat electronic polarization. Towards this end polarizable force fields for lipids have been presented in the context of the fluctuating charge42 and classical Drude oscillator models.41 Notably, an MD simulation using a preliminary polarizable Drude model was able to more accurately reproduce the dipole potential for a DPPC monolayer.43 Insights gained from polarizable force fields may facilitate further optimization of additive models.
An ongoing problem for most of the present FFs is the reliance on 1980’s water models such as TIP3P (for CHARMM) and SPC/E44 (for GROMOS). For example, the TIP3P model underestimates the viscosity of water by approximately a factor of 3 at room temperature, and poorly reproduces the temperature and solute dependence of viscosity.10 These errors impact the simulated dynamic properties of solutes in solution. While corrections can be made to account for this issue,10 they represent further approximations. However, correcting these and related problems is a significant challenge because the water model is deep within the fabric of a FF and cannot be replaced simply. While next-generation polarizable FFs may be able to overcome these problems (e.g., the DRUDE SWM4-NDP model45 reproduces a range of gas and condensed phase target data), it may be considered desirable to develop new water models with improved properties that are consistent with available additive FFs.
Coarse-grained models have achieved remarkable success despite their simplicity. This is partly because they capture the aspects of the hydrophobic interaction (lipid tails and water are parametrized not to mix) and steric repulsion. However, they are limited by a lack of detail. It is thus reasonable to consider multiscale modeling, in which a patch of membrane is described with all-atom models and the surrounding area is coarse grained. The potentials will need to be expanded to include such possibilities.
Finally, efforts are ongoing to generate parameters for more complex lipids including cardiolipin, phosphoinositols, and sphingolipids. Towards this end, an extensive set of parameters for carbohydrates has recently been developed46–49 and represents an important extension of the CHARMM all-atom additive FF. In combination with the available protein and lipid portions of the FF, the availability of the carbohydrate parameters will allow for simulation studies of glycoproteins and glycolipids. These capabilities now allow the CHARMM additive FF, including the CHARMM General FF50 for small molecules, to represent a comprehensive model for the study of heterogeneous biomolecular systems using simulation technologies.
We thank our many collaborators involved in this now two-decade project, especially Bernard Brooks, Scott Feller, Jeffery Klauda, Benoit Roux, Doug Tobias, and Richard Venable. Experimentalists Michael Brown, Klaus Gawrisch, Stuart McLaughlin, John Nagle, and Stephen White have also been indispensible for providing data, often before publication, and helping us understand it. Financial support from the NIH (GM070855, GM072558 and GM015101 to ADM), and the University of Maryland Computer-Aided Drug Design Center are appreciated. This research was supported in part by the Intramural Research Program of the NIH, National Heart, Lung and Blood Institute.
Richard W. Pastor received a B.A. in philosophy from Hamilton College (1973), a M.S. in chemistry from Syracuse University (1977), and Ph.D. in Biophysics from Harvard University (1984). He did research and review at the Food and Drug Administration from 1984 to 2007, and moved to the National Institutes of Health in 2007. He presently is Chief of the Membrane Biophysics Section in the Laboratory of Computational Biology, National Heart Lung Blood Institute. http://www.lobos.nih.gov/mbs/index.shtml
Alexander D. MacKerell, Jr. received a Ph.D. in Biochemistry in 1985 from Rutgers University, which was followed by postdoctoral fellowships in the Department of Medical Biophysics, Karolinska Intitutet, Stockholm, Sweden and the Department of Chemistry, Harvard University. In 1992 he assumed his faculty position in the School of Pharmacy, University of Maryland where he is currently the Grollman-Glick Professor of Pharmaceutical Sciences and the Director of University of Maryland Computer-Aided Drug Design Center. http://mackerell.umaryland.edu/MacKerell_Lab.html