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
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
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 (). 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.
Table 1 Deuterium order parameters, |SCD|, for glycerol carbons, and carbon 2 of the sn-1 and sn-2 chains of DPPC bilayers at 323 K from experiment,28,29 and four successive CHARMM lipid parameter sets. (See for nomenclature.)
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 (). 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 (); 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.