|Home | About | Journals | Submit | Contact Us | Français|
A significant modification to the additive all-atom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gel-like structures of bilayers well above the gel transition temperature, selected torsional, Lennard-Jones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from high-level ab initio calculations on small molecules to semi-empirical QM studies on a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: DPPC, 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); simulations of a low hydration DOPC bilayer were also performed. Agreement with experimental surface area is on average within 2%, and the density profiles agree well with neutron and x-ray diffraction experiments. NMR deuterium order parameters (SCD) are well predicted with the new FF, including proper splitting of the SCD for the aliphatic carbon adjacent to the carbonyl for DPPC, POPE, and POPC bilayers. The area compressibility modulus and frequency dependence of 13C NMR relaxation rates of DPPC, and the water distribution of low hydration DOPC bilayers also agree well with experiment. Accordingly, the presented lipid FF, referred to as C36, allows for molecular dynamics simulations to be run in the tensionless ensemble (NPT), and is anticipated to be of utility for simulations of pure lipids systems as well as heterogeneous systems including membrane proteins.
This paper presents the CHARMM36 (C36) potential energy parameter set (or force field, FF) for lipids. C36 is an additive, all-atom model that replaces CHARMM27 (C27)1,2 and CHARMM27r (C27r)3,4 (a revision to C27 updating the acyl chain torsions, but leaving other parameters unchanged); it is consistent with CHARMM (Chemistry at HARvard Molecular Mechanics) parameter sets for proteins,5,6 nucleic acids,7 ethers,8 carbohydrates,9–12 and small molecules.13,14 Accordingly, the C36 FF is applicable for simulation studies of heterogenous biomolecular systems.
Development of C36 was primarily motivated by two flaws in both C27 and C27r.3,15–18 First is the large positive surface tension (30–40 dyn/cm) for fluid phase bilayers at the experimentally determined surface area per lipid. While there is not a direct measure of the surface tension, γ, of microscopic bilayers, theoretical considerations for self-assembled systems indicate that γ is identically zero19–21 when the bilayer is flat, or only a few dyn/cm22 when undulations are taken into account. Experimental measurements on macroscopic black lipid bilayers also indicate a surface tension in the neighborhood of a dyn/cm.23 The consequence of this positive surface tension with C27 and C27r is that bilayers shrink to near gel-phase areas when simulated at zero surface tension (as in the NPT, or constant number, isotropic pressure and temperature, ensemble).17,18,24 The ad hoc solution of simulating with a large applied surface tension is unsatisfactory, given the good theoretical arguments for a much lower or vanishing surface tension.
The second flaw is that simulations carried out with these sets do not reproduce the experimental deuterium order parameters, SCD, in the glycerol and upper chain regions.3 In particular, the splittings observed for a wide range of glycerophospholipids25–29 for Carbon 2 (C2) of the aliphatic chains and Carbon 1 (CG1) of the glycerol are virtually absent. The inability of these FFs to reproduce the fine structure at the lipid water interface and the inequivalence of the acyl chains raises concerns that interactions of lipids with surface active agents would likewise be incorrectly described. Other flaws in the these FFs include underestimation of the area compressibility modulus, KA,30 underhydration of the head group region,31 underestimation of the electron density in the bilayer midplane,16 and overestimation of the frequency dependence of the 13C NMR T1 of the acyl chains near the head group.27
The strategy for the reparametrization was based on first reevaluating the partial atomic charges. Recent studies by different groups17,24,32 indicated that changes in the partial atomic charges in the head group of C27 lead to substantial reductions in surface tension for 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC), 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC) bilayers. The present reevaluation began with a determination of the validity of applying charges based on small model compounds (Fig. 1) in vacuum to DPPC in a bilayer using semi-empirical (AM1 level33) quantum mechanical (QM) calculations. Based on the relative independence of charges with system size, application of the small molecule based optimization of non-bond parameters common to the CHARMM additive force fields was applied.34 In addition, high-level ab initio calculations were carried out on model compounds representative of the lipid head group and glycerol linker region to act as the basis for optimization of selected torsion angle parameters. Optimized parameters were then tested on simulations single component lipid bilayers containing DMPC, DPPC, and POPC, 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE). Simulations of DLPC, DMPC, DPPC, DOPC, POPC, and POPE bilayers were carried out at high hydration and DOPC was also subjected simulations at low hydration. Following several iterations of testing small molecule based parameters on DPPC bilayer simulations, the final parameter set was validated by comparisons with experimentally determined surface areas, compressibilities, deuterium order parameters, 13C NMR relaxation rates, and neutron/diffraction profiles on various bilayers.
Although the focus of this work is improving the all-atom CHARMM FF, it is important to review briefly the current pair-wise additive lipid force fields that are used in molecular simulations. These atomic- or near atomic-level FFs can be classified into either all-atom or united-atom models. In addition to the abovementioned CHARMM FFs,1–4 the all-atom general AMBER force field (GAFF)35 has the capability to simulate lipid membranes. GAFF has been successfully tested on pure fully-hydrated DMPC,36 DOPC,36,37 and POPC38 bilayers with a decent description of structural and dynamical properties for a force field not fit to lipid experimental data. Notably, GAFF has been shown to be the only FF to capture the carbon-2 deuterium order parameter splitting,37 but does result in lower surface areas per lipid and higher deuterium order parameters compared to experiment.36–38 Recently, Hénin et al.39 developed a united-atom acyl chain FF for CHARMM phospholipids, while maintaining an all-atom description of the lipid head group. This matches the C27r FF, while reducing the number of atoms, thus making it attractive for larger heterogeneous membranes. One of the main united-atom FFs for lipids are the GROMOS-based FFs, which have been optimized to reproduce condensed phase properties of alkanes.40 One commonly used GROMOS-based FF was developed by Berger et al.,41 which included optimization of the LJ parameters to reproduce condensed phase properties of pentadecane. More recent modifications to the united-atom GROMOS FFs are 43A1-S342 and 53A6.43 Adjustments to the charges and LJ parameters of choline methyls and phosphate oxygens resulted in good agreement for surface area per lipid, density profiles, and deuterium order parameters for saturated and unsaturated acyl chain PC lipids with GROMOS-53A6.43,44 However, the isothermal area compressibility modulus was found to be nearly twice that experimentally measured. Besides the GROMOS-based force fields, several research groups have developed their own united-atom lipid FFs, such as, Smondyrev and Berkowitz,45 OPLS-UA,46 and the adoption of the NERD FF to fatty acids by Sum et al.47 Clearly there is a wide variety of united-atom force fields for lipids available, while only two all-atom models are available, GAFF and CHARMM, with CHARMM being the only FF parametrized directly for lipid molecules. Moreover, the different methodologies used to develop these individual FFs preclude the ability to combine two FFs from varying research communities.
By way of an outline, the Methods section initially describes the nomenclature used for the lipids and small molecules. Two additional subsections of the Methods describe the QM methods used in FF development and methods used in the molecular dynamics simulations. Section 3 provides a detailed description on how the FF parameters were fit to experimental and QM data, including adjusting partial charges to model the complex environment of a lipid bilayer (3.1), and conformational energy studies on small molecules representative of various lipid moieties (3.2). Section 4 describes the results from simulations with C36 for the model lipid DPPC (4.1), and other lipids (4.2). Section 5 examines two limitations of the new FF: the dipole potential drop (5.1); and inconsistencies of surface tensions of bilayers and monolayer (5.2). Section 6 presents the discussion and conclusions.
The Methods section is divided into three subsections. In the first, lipid molecule and atom names are briefly reviewed to avoid confusion among several commonly used conventions. Section 2.2 describes the QM methods used in justifying modifications to the partial charges for certain atoms in the head group of the lipid molecules along with the optimization procedure used to produce optimized non-bond parameters for the ester moieties in the glycerol linker region. This subsection also describes the QM methods used to obtain highly accurate torsional profiles for certain head group and alkene dihedrals. The Section 2.3 describes the methods used in the molecular dynamics (MD) simulations. The C36 additive all-atom lipid force field may be downloaded from the MacKerell group web page at http://mackerell.umaryland.edu/CHARMM_ff_params.html.
Lipid molecule names are in IPUAC sn nomenclature48,49 (see http://www.chem.qmul.ac.uk/iupac/lipid/) throughout this document. For technical reasons, each atom of a molecule in a molecular mechanics FF must have a unique and simply expressed name, generally using the letters A–Z and numerals; standard sn nomenclature references do not address atom names in this level of detail. The CHARMM lipid FF uses atom and torsion names based on the Sundaralingam nomenclature:50 a significant difference is the numbering of the glycerol carbon atoms, which is inverted with respect to the sn nomenclature. For glycerol-phosphate itself, 1-phospho-D-glycerol is the more correct name, while sn lipid molecule names are based on 3-phospho-L-glycerol, which is identical in terms of absolute configuration. To help clarify the equivalent names for a given atom, both sn and FF atom names are included in Tables 2 and and5.5. The sn atom nomenclature is used in the figures and text in the remainder of the document. The Supporting Information uses FF names in the tables, and sn names in the figures. For torsions, the Sundaralingam names will be used throughout; they are based on the chain names, which are the first three letters of the Greek alphabet, corresponding to the Sundaralingam glycerol carbon numbering. The polar chain is designated the α chain, the sn-2 chain corresponds to the β chain, and the sn-1 chain corresponds to the γ chain. Finally, the polar α chain atoms of both DOPC and POPC were renumbered in the C36 lipid FF release for consistency with other PC-based lipid descriptions.
To address the transferability of partial atomic charges based on small molecules to the larger head group in lipids, a series of semiempirical QM calculations on different compounds in different environments were undertaken. These calculations used the program DivCon,51 which applies a divide and conquer approach52,53 to allow for full AM1 level calculations to be performed on systems with 10,000 or more atoms. Both the Mulliken and CM2 charges54 were obtained and compared for the different compounds and environments to determine the impact of environment on partial atomic charges. DPPC was selected as the target lipid and corresponding small molecules were chosen, i.e., methylacetate (MAS), dimethylphosphate (DMP) and ethyltrimethylammonium (ETMA) (Fig. 1). These molecules were selected as they represent the compounds used for optimization of the non-bond parameters in the C27 FF, with DMP and MAS also being used in the present work. In addition, EGLY and PGLY (Fig. 1) were analyzed to determine the potential utility of larger model compounds. All model compound results were based on gas phase calculations and charges were averaged over multiple conformations as described in the legend of Table S1. For the full lipid, calculations were performed on the DPPC molecules individually in the gas phase and, importantly, on the fully solvated lipid bilayer. The charges for individual lipids used for averaging were based on lipids that were located near the central region of the bilayer to include realistic membrane interactions. The solvated bilayer was obtained from an NPAT MD simulation of 72 DPPC molecules at high hydration and the experimental surface area in periodic boundary conditions using the C27r FF.3 From this simulation, 13 time frames were extracted and individual lipid molecules in the central region of each leaf of the bilayer were identified, yielding a total of 64 lipid conformations. Each of these conformations was subjected to an AM1 calculation isolated in the gas phase and average charges were calculated. To evaluate the impact of the condensed phase environment on the partial atomic charges, AM1 calculations using DivCon on the full bilayer-solvent systems were calculated to obtain charges in a membrane environment. Though these calculations were performed in the absence of periodic boundary conditions used in the MD simulations, the conformations of the lipids in the gas phase and the bilayer/water phase calculations were identical; only the environment of the individual lipids changed.
The dihedral portion of the CHARMM FF is mostly from direct fits to vacuum QM-based conformational energies. This approach was not satisfactory for the glycerol torsions β1θ2, and θ4, where solvent effects are expected to be complex such that gas phase potential energy surfaces may be inappropriate.55,56 As described below, these torsions were adjusted slightly to match the experimental deuterium order parameters glycerol region for DPPC following new QM estimates. (Values of SCD reported here are always positive; i.e., the absolute value sign is assumed.) The QM calculations follow the approach used for C27r where the acyl torsions were refined; see references 3,57 for details. Figure 1 shows the small molecules used to represent regions of the phospholipid to allow for high-level QM calculations. Second order Möller-Plesset perturbation theory (MP2)58 is used to obtain global minimum energy geometries of these small molecules and a single dihedral constraint is used to obtain torsional energy scans. To obtain conformational energy at the coupled cluster level, CCSD(T),59 with a large basis set, HM-IE (the Hybrid Method for Interaction Energies)60 is used to estimate CCSD(T)/cc-pVQZ,
which assumes that the effect of larger basis sets is additive. These conformational energies using HM-IE will be simply referred to as QM in this paper.
The following objective function was used initially to obtain all new torsional terms for C36,
where Ui is the energy for conformation i. The alkane torsions3,4 were also refit to this function to maintain a consistent approach for fitting dihedral parameters (a higher weight was used for the trans and gauche minima in C27r). The alkene torsions adjacent to C=C doubles bonds were fit to conformations of 2-hexene at the RIMP2/cc-pVQZ//MP2/cc-pVDZ level using a Monte Carlo Simulated Annealing approach.61
Small molecules used in the parametrization of the C36 were simulated at 298.15 K and 1 atm pressure using the new velocity Verlet integrator62 implemented in CHARMM. Standard procedures were used to include the influence of long-range electrostatics63 and Lennard-Jones (LJ) forces.64,65 Heats of vaporization and molecular volumes of MAS were calculated using a box of 128 molecules. Free energies of aqueous solvation of MAS and DMP were calculated via free energy perturbations (FEP)66 using the staged protocol developed by Deng and Roux,67 as previously described.68 More details on these methods to obtain the heat of vaporization and free energies of solvation are provided in the supplementary material.
The CHARMM69,70 and NAMD71 programs were used to simulate lipid bilayers and monolayers. The lipids studied, total simulation times, temperatures, and other conditions are listed in Table 1. A LJ switching function over 8 to 12 Å was used in MD simulations with CHARMM. For the NAMD simulations, a shorter LJ switching function was used (11 to 12 Å). Particle mesh Ewald (PME)63 was used for the long-range electrostatic contribution for both CHARMM and NAMD. The LJ and PME contributions to the energy and forces were updated every step for the CHARMM simulations and every two and four steps, respectively, for the NAMD simulations (impulse-based Verlet-I/r-RESPA method).72,73 The temperature was held constant with the Hoover thermostat74 in CHARMM and by Langevin dynamics with a weak coupling constant of 1 ps−1 in NAMD. Similarly the pressure was maintained anisotropically (NPT ensemble simulations) or in the direction perpendicular to the membrane (NPAT ensemble simulations) by the Nose-Hoover piston75,76 (CHARMM) or a Nosé-Hoover-Langevin piston77,78 (NAMD) at 1 bar. These differences between MD simulations with CHARMM and NAMD encompass typical approaches used by various research groups. Consistent for all of these simulations was the use of a 1 fs time step and constraining of the hydrogen atoms using the SHAKE algorithm.79
Starting structures for DMPC,16 DPPC,3 POPE80 and low hydration DOPC18 simulations were obtained from previous simulations with the C27 or C27r force field. For DPPC isotherm configurations other than 64 Å2 area/lipid, two NPγT ensemble MD simulations starting from 64 Å2 with applied surface tensions (γ) of −10 and 30 dyn/cm were used to contract and expand the system, respectively. The coordinate sets with unit cell edge lengths that most closely matched that for each target area were extracted from these simulations, and used as simulation starting points. The DLPC, POPC, and high hydration DOPC bilayers were built de novo from a library of single lipid conformations, as described previously.81 In this case, the library of conformations for each lipid was generated from eight 2-ns Langevin dynamics simulations in an orienting potential, with a different random seed for each replicate simulation. The P atom of the phosphate was tightly restrained to a position along the +z axis commensurate with the experimental bilayer thickness for each lipid. To prevent free rotation and limit whole molecule motions to wobbling in a cone, three restraint planes were used: the choline N atom was restrained to be 1 Å above the P atom, the ester O atoms were restrained to be 2 Å below the P atom, and the chain terminal C atoms were restrained to be within 4 Å of the z=0 plane. The restraint planes functioned as soft walls at the indicated z positions.
This section describes the procedure of adjusting certain parameters to the CHARMM FF,
where the terms have their usual meanings. The bond and angle terms are left unchanged from previous modifications. As opposed to recent charge modifications to the CHARMM lipid FF by other workers,17,24,32 the parametrization discussed in this section is consistent with the methods used in the CHARMM additive FF development strategy.20 The first subsection justifies the need for changing the non-bonded parameters of the lipid head group (last two terms in Eq. 3). The choice of the final set of new lipid head group charges and Lennard-Jones (LJ) terms are also described in this subsection. The second subsection describes the dihedral changes to the CHARMM FF in C36.
Optimization of lipid atomic charges depends on the size of the compound used to represent the full lipid molecule and its environment. A systematic test of this size dependence is illustrated in Table 2, which lists the average partial atomic charges based on the AM1/CM2 model of selected atoms in the model compounds and the differences in the partial atomic charges for the larger model compounds and the full lipid molecules. The average partial atomic charges and their standard deviations of all the compounds are presented in Table S1. The Mulliken and CM2 based charge schemes result in similar trends (Table S2). Overall, the charges for the model compound are similar to those for the larger compounds, i.e., differences are typically 0.05e or less. However, there are some significant differences for the ester charges (sn-2 chain C1, O1, and O2). For example, in the larger compounds the carbonyl oxygen (sn-2 O2) becomes more negative indicating local polarization, though the carbonyl carbon (sn-2 C1) and the ether oxygen (sn-2 O1) become more negative and more positive, respectively, which minimize this effect. However, upon going from the gas to the bilayer/water phase for the full lipid, the negative charge on the sn-2 C1 becomes more positive while that of the sn-2 O2 becomes more negative for both ester moieties. Thus, the present QM results (Table 2, S1, and S2) indicate that local charge induction occurs in the complex environment of the lipid bilayer, and that subtle differences observed as a function of model compound and environment should be considered when selecting the final non-bond models.
MD simulations of DOPC at low hydration indicated that the head group was not adequately hydrated.31 Such limitations may be probed by considering the free energy of aqueous solvation of model compounds representative of those regions of the lipids. Studies have shown that the free energy of aqueous solvation of tetraethylammonium is in excellent agreement with experiment,82 indicating that the quaternary amine parameters are satisfactory. This implies that the non-bond parameters in the linker region of the FF, including the phosphate and esterified glycerol moieties, may be problematic. Calculations were undertaken on the model compounds DMP and MAS, with respect to both QM and condensed phase calculations. Particular emphasis was put on MAS because that model had not been assessed since the original lipid force field was published in 1996,83 whereas the phosphate parameters had been investigated in detail in 2000 in conjunction with nucleic acid FF development.7
A detailed analysis on the ability of the DMP FF model by Foloppe and MacKerell,7 showed that the model satisfactorily treated the interactions with water as a function of orientation as judged by QM calculations. Accordingly, emphasis in this study was placed on evaluation of the DMP model in the condensed phase, targeting the free energy of aqueous solvation, ΔGS. A direct experimental value for ΔGS of DMP is not available; as described in the Supporting Information it was estimated to be −76±4 kcal/mol using a thermodynamic cycle.84,85 Calculation of ΔGS of DMP for C27 and C36 was performed using a free energy perturbation methodology based on the Weeks, Chandler, Anderson decoupling approach (see supplementary material).67,86 Error estimates were obtained by repeating the calculation three times with the individual simulations seeded with different initial velocities. From these calculations, −76.7±0.05 kcal/mole was obtained for C27, which is in excellent agreement in the estimated experimental value (Table 3). Based on this high level of agreement, the published DMP non-bonded parameters were maintained with the exception of LJ parameters for ester oxygens, which were modified during MAS optimization as described below. This change has only a small effect on DMP solvation energetics (cf. C27 and C36 values in Table 3), and the new value of −77.7 kcal/mol is within the statistical uncertainty of the experimental estimate.
The original optimization of the MAS non-bond parameters focused on the interactions with water, targeting QM data, and on the reproduction of its experimental pure solvent properties. 83 In the present study, these target data were reanalyzed and the optimization was extended to include the dipole moment and the experimental ΔGS. For the interactions with water, the QM MAS-water monohydrate interaction pairs were optimized at the HF/6-31G* level, in which the minimum interaction distances are offset by 0.2 Å and the interaction energies scaled by 1.16 as required to mimic local polarization effects in the condensed phase as well as limitations in the QM level of theory. Importantly, this level maintains consistency with the remainder of the CHARMM additive all-atom biological force field.5 These energies and dipole moments were used as target data to reoptimize the MAS partial atomic charges while reproduction of the pure solvent heats of vaporization and molecular volumes were used to adjust the LJ parameters. The new LJ parameters included the recently optimized oxygen LJ parameters for the ether oxygen in the ester linkage.8 The methyl LJ parameters were not adjusted. Once satisfactory models were obtained, free energy of solvation calculations were performed and used to select the final parameter set.
Table 4 lists the MAS dipole moments and interaction energies with water based on the QM calculations and C27 and C36. The C27 energies agree well with the QM target data. In contrast, the dipole moment, which was not previously included as target data, shows significant disagreement with the QM target values with respect to both the magnitude and the orientation. Accordingly, further optimization of the charges targeted the dipole moment along with the interactions with water. LJ parameters were simultaneously optimized to overcome the significant overestimation of the MAS pure solvent heat of vaporization with C27 (Table 3). This led to the development of several non-bond models for MAS. These models focused on accurate reproduction of the heat of vaporization and dipole moment with the level of agreement with the water interaction energies varied with respect to the magnitude of the energies, as compared to the QM data.
From this effort a model was selected in which the experimental ΔGS was overestimated, motivated by the QM charge analysis presented above and the need to obtain better agreement with the solvation of the head group and linker regions of DOPC at low hydration. The results indicated that a more favorable ΔGS is necessary to properly hydrate the head group region (Table 3). Interestingly, application of the charges from a previous reparametrization of the DPPC parameters17 to MAS as part of the present study yielded ΔGS =4.8 kcal/mol suggesting that the success of the previous model was also based on increased solvation of the glycerol linker region. This increased solvation was noted in the previous work with change attributed to differences in the charges of one of the methylene groups.17 At this stage several sets of ester non-bond parameters were used in MD simulations of DPPC bilayers from which the final non-bond model of the ester moiety was selected based on the surface tension of the simulation approaching zero (not shown). Table 5 compares the final lipid charges for C36 with those of C27r for the linker region. The final model yields good agreement with respect to the dipole moment, including its orientation, and the experimental pure solvent heat of vaporization and molecular volume (Table 3). This non-bond model of the ester moiety was used as the starting point for final optimization of the dihedral parameters and subsequent extensive testing in lipid bilayer simulations.
The conformational QM energies for the alkane model compounds, IPB (β4), and PB (γ4) (Figure 1) were previously calculated3,57 and potential parameters were fit to these QM energies (Table S3). The alkane torsions were refit using an objective function consistent with that used for the other torsional fits performed in the present study (Eq. 3), but this had a minimal effect on the torsional profiles. 2-hexene was also used to refit the torsion parameters of the two C-C single bonds adjacent to C=C double bonds and C36 fits agree well with the QM data (Figure S1, Supporting Information). Conformational energies for γ1 and γ3 were calculated using the MEGLY model compound and PMP was use to calculate theα4 torsional profile (Figure 2). There are significant improvements over C27r for the torsional barriers and gauche-like conformational energies, similar to previous alkane torsions. Specifically, C27r incorrectly predicts the trans state as the most stable conformation for the choline α4 torsion, but the gauche+ conformation is the most stable. The C36 fit to the QM energies resolves the inaccurate relative conformational energies that exist with the C27r FF. This change substantially influenced the calculated deuterium order parameters for lipids, as discussed in the following section.
For β1, θ2, and θ4, additional adjustments were based on MD simulations of DPPC bilayers to evaluate torsional populations that best matched deuterium order parameters for the glycerol and acyl chain C2 positions. A 27 state model, based on assuming 3 minima for each of the 3 dihedrals, was used to assist in the evaluation. Lipids were grouped into one of the 27 states based on the values of the 3 dihedrals, and the average SCD from each group was computed. The experimental SCD for a given C-H vector then could be calculated from a population weighted sum of the SCD values for each of the states. While attempts to predict optimal populations by fitting were not successful, this exercise indicated that a major reduction in the population of the trans conformation for θ2 was required, and that shifts in the conformer populations ofθ4 should be considered. Four trial parameter sets were created and evaluated by 25-ns MD simulations of DPPC in the NPAT ensemble. The set with the best match to the experimental SCD values was chosen. Comparison of QM and C36 potential energy surfaces for the β1 and θ4 surfaces based on the MEGLY model compound are shown in Figure S2 (left) of the Supporting Information. While final optimization of the associated dihedral parameters targeted the experimental SCD values, the level of agreement between the QM and final empirical model is satisfactory. This indicates that the use of the experimental target data is not leading to the unphysical sampling of high-energy regions of the potential energy surfaces. In addition, the large discrepancy between QM methods (Fig. S2, left) suggests that the uncertainty for these two torsions is high.
For the first stage of the validation of the force field parametrization, simulations of DPPC bilayers were performed at a constant surface area of 64 Å2/lipid (NPAT ensemble), which is close to the experimental value of 63±1 Å2/lipid.87 As discussed in the introduction, the deuterium order parameters, SCD, from the C27r simulations lack proper splitting of the first carbon on the aliphatic chain (sn-2 C2) and demonstrate large disagreement with the glycerol (CG1–CG3) carbons (Fig. 3, top). The results based on the C36 FF are a dramatic improvement for these carbons (Fig. 3, bottom). The SCD for C2 of the sn-1 chain should be higher than those of the sn-2 chain and simulations with C36 agree with these measurments.25,26 There is some splitting of the two C2 hydrogens of the sn-2 chain, albeit not to the extent of experiment. Order parameters for the CG1–CG3 carbons with C36 are substantially better than those from C27r. Since the three torsions that primarily control these order parameters (β1, θ2, and θ4) were empirically optimized, this agreement is expected. However, to our knowledge this is the first FF to accurately capture the chain order splitting for all carbons. The contour plots in Fig. 4 highlight the complexity of the glycerol torsion surfaces, and the qualitative differences between C27r and C36. Moreover, the θ4/β1 torsional surface of C36 compares well with that of MP2/6-31g(d) scans (Fig. S2, right).
The SCD of the α carbon of the choline is reduced with C36 and is the result of ab initio-based torsional fits to the α4 dihedral. The carbonyl C=O and ester C-O order parameters (0.26 and 0.15, respectively) for the sn-2 chain have been measured using infrared spectroscopy.88 MD simulations with C36 agree well with these experimental values (0.26/0.14) and are an improvement compared to C27r (0.18/0.13). These torsional changes for the head group have a minimal effect on the aliphatic chain and these SCD are nearly identical between DPPC bilayer simulations with C27r and C36.
The main issue with previous versions of the CHARMM lipid FF was that NPT bilayer simulations with fully saturated lipids, such as DPPC, condensed to surface areas near that of the gel state above the transition temperature.17,24 As shown in Fig. 5 and Table 6 for DPPC and the other lipids under study (discussed in detail below), this flaw is almost entirely rectified. MD NPT simulations for DPPC using two cutoff schemes (see Methods) also result in stable surface areas that are close to the experimental surface area for simulations up to 130 ns (Table 6 and Fig. S3). Since the tensionless bilayer has an average surface area close to experiment, the SCD calculated from NPT C36 simulations are similar to the NP(A=64)T simulations and agree with experiment (Fig. S4).
The order parameters are an excellent measure of molecular order but do not describe the overall structure of the bilayer. For this, the bilayer simulations are compared to x-ray diffraction data, i.e., form factors and electron density profiles as obtained by the SDP structural model 87 (Figs. 6 and and7).7). The form factors offer a direct comparison with experiment that prevents errors that might be associated with structural models or inverse Fourier transforms.16,18,89 There is excellent agreement between experimental form factors, F(q), and bilayer simulations with C36 (Fig. 6). The root mean squared deviation (RMSD) from experimental F(q) for simulations with the C27r FF is less than C36 with a surface area of 64 Å2/lipid. The minimum RMSD for C36 was for the 66 Å2/lipid, which is likely the result of a slight reduction of the electron-rich head group density (Fig. 7). There is a significant improvement in the electron density for the hydrocarbon region with C36 (Fig. 7). The density at the center of the bilayer is too high with C27r and is reduced with C36 with the NP(A=64)T and NPT simulations. The carbonyl-glycerol (CG) distribution is broader and slightly shifted with C36 compared to C27r, but the phosphate-choline (PC) distribution is improved.
The excellent agreement between simulated and experimental surface areas for DPPC is borne out by the area isotherm at 323.15K where the zero surface tension lies between 62 and 64Å2/lipid (Table 6). This simulated data also allows estimation of the area compressibility modulus, KA, from the relation30
where A0 is the equilibrium area, here set to the experimental value of 63 Å2/lipid; the factor of 2 is required when surface tensions for bilayers are expressed in dyn/cm/side Using all 5 data points (60 to 68 Å2/lipid) and their standard deviations yields KA = 207± 14 dyn/cm; KA= 236 ± 31 dyn/cm when only the three middle points (62 to 66 Å2/lipid) are used. These are in substantially better agreement with the experimental value of 234 dyn/cm90 than the 138 ± 26 dyn/cm obtained for C27r.30
Fig. 8 compares the spin lattice relaxation rates from 13C NMR27,91 experiments with those calculated with the C27,92 C27r27 and C36 FFs. The C36 calculated NMR T1 for the head group carbons (Fig. S5) are in excellent agreement with experiment and are similar to C27r (NPAT). As shown previously,27 C27r is a significant improvement over C27 due to changes in the torsional barrier for alkanes. C36 is in similar agreement for C14 to C16. Although C27r and C36 both agree with the high frequency measurements for all aliphatic carbons, C36 significantly improves the frequency dependence at low Larmor frequencies for the center of the chain and those near the glycerol (C2 and C3).
Although the C36 FF was successful at reproducing many experimental observables for DPPC bilayers, additional lipid bilayers were simulated to test the transferability of these parameters to various systems. An experimental observable for fully saturated chains of lengths 12–16 carbons in PC-containing lipids is the surface area per lipid, which has a minimum when the aliphatic chain contains about 14 carbons (DMPC).16,87,93 MD simulations with C36 and the NPT ensemble accurately capture this dependence of surface area per lipid on chain length (Table 6). Moreover, the agreement between experiment and simulation surface areas is excellent for the NPT simulations of DPPC, DMPC, and DLPC. The density profiles and experimental form factors are in similar agreement for the DMPC bilayers (Fig. 6, bottom) using both FFs. However, MD simulations with C27r require constant area simulations, whereas C36 does not.
Saturated lipid chains are common in many cellular membranes, but lipids with unsaturated chains are also prevalent. Therefore, a successful FF should also reproduce experiments on these lipids. POPC contains one chain that is fully saturated and another chain with a single double bond. The surface area per lipid for POPC using C36 and the NPT ensemble (64.7±0.2Å2/lipid) is slightly lower than the experimental estimate 68.3±1.5 Å2/lipid.94 The NMR deuterium order parameters with C36 NPT simulations are in excellent agreement with experiment28,29 (Fig. 9, top). The agreement is nearly perfect for the fully saturated chain (sn-1), including the dip carbon position 3. C36 also reproduces the splitting of the order parameters at carbon position 2 for the sn-2 chain, and the characteristic drop large around the double bond.
For a lipid bilayer with both chains unsaturated, i.e., DOPC, the C36 NPT simulations are in good agreement in the lateral and profile densities. The average surface area per lipid is slightly higher than experiment (Table 6), though within 2 standard errors of the experimental uncertainty. The overall density profiles of C36 agree favorably with the experiment (Fig. 10).87 The density is slightly higher to the right of the peak and likely the result of a broader carbonyl-glycerol distribution and increase in water hydration. The overall density in the hydrocarbon region of the bilayer is in excellent agreement with experiment. However, the individual component densities disagree. The location of the carbons associated with the double bond (CH) is closer to the bilayer center than the SDP model. The SDP model is based on a combination of x-ray and neutron scattering measurements. Surprisingly, the simulation CH density is in better agreement with the neutron results,87 where models fit to this data place the double bond closer to the bilayer center. This discrepancy between simulation and experiment will require a more detailed study to address whether this is a slight issue with simulation or the experimentally-based SDP model.
All of these abovementioned simulations were based on fully hydrated lipid bilayers. MD simulations of DOPC at 66% relative humidity (RH) or low hydration (lo) were used to further test NPT simulations with C36. Previous simulations using the C27 FF resulted in condensation and low surface areas per lipid for DOPC using the NPT ensemble (56.5±0.3 Å2/lipid),18 while using the C36 FF there is excellent agreement between simulation and experiment (Table 6).
The spatial distribution of selected moieties along the lipid normal can be determined using neutron scattering experiments. This is due to the ability of neutrons, in contrast to X-ray based methods, to detect hydrogen and to differentiate between hydrogen and deuterium. The extent of water penetration into the lipid bilayer normal can be probed using neutron scattering. The Gaussian water distributions from the experimental work of Weiner et al.95 along with those from MD simulations of DOPC using C27 and C36 are compared in Figure 11. It is evident that the distribution is broader and is shifted to distances further away from the center of the lipid bilayer for C27r as compared to experiment. These differences highlight the limitations in the C27 FF associated with the solvation of the head group region, and led to the reevaluation of the ester non-bond parameters presented above. In the resulting C36 FF, the simulated distribution is narrower and shifts back towards the center of the lipid bilayer. The better agreement with the distribution obtained from the neutron scattering experiments with respect to both mean and width indicates that the optimized parameters lead to improved solvation of the phosphatidyl and acyl ester regions of the lipid, which, in turn, helps maintain the area/lipid near the correct values.
The last test of the C36 FF was on a lipid bilayer with a different head group, i.e., phosphatidylethanolamine (PE). The average lateral surface area per lipid for MD simulations with C36 and NPT was 59.2±0.3 Å2/lipid, which agrees with x-ray measurements (Table 6).96 Since the experimental areas are at different temperatures than simulation and the resolution is lower than other lipids reported above, another way to test the validity of the lateral surface area is to compare with NMR deuterium order parameters, which are sensitive to the average area per lipid.97 The agreement between SCD values from experiment28,98 and C36 NPT simulations is excellent for the sn-1 and sn-2 chains (Fig. 9, bottom). A dip in the SCD is not observed for the carbon 3 position of POPE sn-1 chain,98 but a slight lowering of the SCD was calculated with the C36 NPT simulations. Since this non-monotonic decrease in chain order is measured for the POPC bilayers, the POPE measurement might be slightly in error for carbon 3. Aside from this slight discrepancy, C36 NPT simulations agree favorably with the SCD values and do not require NPAT simulations that at 60 Å2/lipid would result in a bilayer surface tension of 27.8±1.2 dyn/cm with the C27r FF.
The precise experimental values of the dipole potential drop (Δψ) for bilayers and monolayers remain somewhat controversial,99 and the simulated values are very sensitive to the definition of the charge distribution.100 Nevertheless, evaluation of Δψ provides a useful comparison for parameter evaluation. Figure 12 compares the dipole orientation (top) and profile of the potential (bottom) for C36 and C27r. The location of the maximum of the water dipole in the choline region is identical for the C27r and C36 FFs. However, the angle of the dipole with respect to the interfacial plane is slightly larger at this location, i.e., 15.1° (C36) vs. 12.2° (C27r). Both sets show a similar approximate 200 mV Δψ from the bilayer center (z=0) to the head group region; this is expected because the acyl chains are unmodified (and nonpolarizable). The head group to water Δψ for C36 is approximately 500 mV, which is 200 mV smaller in magnitude than with C27r. These values are, nevertheless, substantially larger than the consensus experimental range of 225–250 mV.99
Simulation results presented so far have been carried out with PME and with a force switch on the Lennard-Jones terms. Hence, while long-range electrostatic interactions have been included, the long-range LJ interactions have been neglected. This is fairly standard methodology, and most simulation programs do not even contain provisions for including long-range LJ interactions. Hence, it is arguable that a parameter set that requires these terms is less generally useful that a set that does not.
However, alkane/air interfaces and lipid monolayers (which contain an alkane-like/air interface) are sensitive to long-range LJ forces.15,101 Consequently, if a lipid FF is to be applied to monolayers, the impact of long-range LJ forces and consistency with bilayers must be considered. Such an analysis was recently carried out on DMPC for C27r using 3D-IPS/DFFT, a method developed by Wu and Brooks102 that can be applied to treat long-range electrostatic or LJ terms. While surface tensions of DMPC lipid monolayers agreed well with experiment when simulated with PME(rc=12 Å), they increased by 8–10 dyn/cm when simulations were carried out 3D-IPS/DFFT. In contrast, the bilayer surface tensions were similar when simulated using the two methods; i.e., the DMPC bilayers did not appear to be sensitive to long-range LJ interactions for C27r. To determine if this trend holds for C36, the bilayers discussed earlier and DPPC monolayers were simulated with PME and 3D-IPS/DFFT.
Table 7 compares the calculated monolayer surface tensions with experiment. In contrast to previous results just described, results with 3D-IPS/DFFT used to treat both the electrostatic and LJ long-range contributions agree with experiment, and those with PME(r=12) are substantially lower than experiment. Some of the long-range LJ interactions that lowered γ for DPPC bilayers have also lowered it for the monolayers. This result implies that bilayer surface tensions for C36 will be too high when long-range LJ are included. This is indeed the case. Simulations of DPPC at areas of 60, 62, 64, 66, and 68 Å2/lipid with IPS-3D/DFFT yield surface tensions of 6.3±0.9, 11.0±1.0, 10.7±0.7, 14.1±0.8, and 16.6±0.7 dyn/cm/side, respectively. Similarly, the surface areas of DPPC and other bilayers contract when they are simulated at NPT and 3D-IPS/DFFT (data not shown).
A CHARMM-FF consistent approach was used to modify the atomic charges and LJ parameters of the glycerol moiety of a phospholipid (Table 5). The modifications to the CHARMM FF were suggested by semi-empirical QM calculations of the full DPPC lipid in vacuum and in a realistic environment, and the incomplete hydration of the head groups observed in bilayer simulations31 This also agrees with a recent test on the CHARMM FF that suggested head group parameter modifications are needed.103 The dipole moment of methylacetate was correctly described by the C36 charge modification (Table 4) and this resulted in an increase in the dipole moment of the carbonyl of both lipid chains. Consequently, there is an increase in hydration of lipid bilayers with C36 leading to good agreement with neutron scattering experiments on DOPC (Fig. 11). The location of the water distribution is shifted slightly inward compared to experiment and may be associated with the overestimation of the free energy of solvation of methylacetate as presented above. However, this overestimation was motivated by the QM charge calculations and facilitated the improved estimate in the DPPC surface area per lipid with NPT simulations.
A total of six lipids were simulated with the NPT ensemble and the surface areas are in very good agreement with experiment (Table 6 and Fig. 5). Deviations from experimental areas average 2%, with a maximum overestimate of 2.4% (DOPC) and a maximum underestimate of 5.3% (POPC). In contrast, C27r simulations with fully saturated chains condensed to values close to gel-like areas.18,24,104 Although this is not the first effort to obtain a CHARMM-compatible lipid FF that can use a tensionless ensemble,17,24,32 C36 is the only FF that obtains these charges in a CHARMM-FF consistent manner.20 Moreover, the stability of the surface area from these simulations has been verified to 130 ns for DPPC (Fig. S3). The ability of C36 to reproduce the correct area for saturated and unsaturated chains implies that simulations with mixtures of lipids, cholesterol and proteins can be performed with the CHARMM family of FF without imposing a surface tension.
The deuterium order parameters for carbons in the glycerol and head group region of the lipid are well described using the C36 FF. Although the SCD splitting is not perfect for carbons near the glycerol moiety (C2/CG1) (Figs. 3, ,9,9, and S4), C36 is a substantial improvement over C27r, and demonstrates that C36 is fundamentally more sound in describing intramolecular conformations near the glycerol moiety. This splitting is not limited to the target DPPC lipid used to modify the C36 FF, but is also observed for other lipids, such as POPC and POPE.
Although the three abovementioned properties were initial targets for the lipid FF modifications, agreement with experimental measurements were also observed for several other properties. First, the lipid bilayer form factors and their real-space counterparts (electron densities) are in excellent agreement with experiment (Figs. 6, ,7,7, and and10).10). The hydrocarbon density at the center of the bilayer is in near perfect agreement with experiment for both saturated and unsaturated bilayers, which was not the case for simulations for C27r. Moreover, the total density profiles of C36 match well with the experimental models. There is an increase in the width of the carbonyl-glycerol distribution compared to experiment and C27r and a slight increase in water hydration and penetration into the fully hydrated bilayers. However, these discrepancies are small and suggest a more detailed investigation of the source of these differences is needed, which is beyond the scope of this lipid parametrization.
Simulations with the C36 FF result in a significant increase in the rigidity of the membrane to respond to area changes. The C36 calculated area compressibility modulus for DPPC was increased by nearly 100 dyn/cm when compared to C27r. This resulted in a value that was within statistical uncertainty of experiment.90
The calculated 13C NMR relaxation rates for the aliphatic chain and head group are improved with C36. The dependence on the Larmor frequency is better with C36 for aliphatic carbons near the head group (C2 and C3) and those at the center of the chain (Fig. 8). The slow relaxation times are decreased with C36 and thus lead to a decrease in the dependence of the relaxation rates on the Larmor frequency. Although the frequently dependent relaxation rates for C2 remain hard to obtain with simulation, values for the other aliphatic carbons agree well with the experiment.27,91
A limitation of the C36 FF is the approximate treatment of long-range LJ forces. While the inclusion of long-range LJ with 3D-IPS/DFFT in simulations with C36 improves results for monolayers, the increased surface tension renders it unsuitable for bilayers. Hence, simulations of bilayers with C36 should be carried out with PME with rc=10 or 12 Å and no long-range correction for the LJ term. These conditions, however, will lead to underestimates of the surface tensions of lipid monolayers. It is possible that a set of adjustments could be developed so that bilayers and monolayers can be treated consistently with all long-range forces included, though this is not necessarily practical given that most simulation programs do not support evaluation of long-range LJ terms.
In conclusion, C36 is a significant update to the C27/C27r lipid force field. This FF was developed in a CHARMM-FF consistent manner, assuring the utility of the model for the treatment of heterogeneous biomolecular systems. Simulations with C36 improve the hydration of lipids near the carbonyl-glycerol section of phospholipids and allow for the use of a tensionless ensemble (NPT) in MD simulations. Although there are some limitations with respect to the treatment of long-range LJ interactions, this is likely not a major concern for most applications. Future work is currently underway to add additional lipids to CHARMM, such as, sphingolipids, phosphatidylinositols, and those with branched aliphatic chains. Since C36 NPT simulations result in stable liquid crystalline bilayers, future NPT membrane simulations with peptides, membrane proteins, and lipid mixtures should not result in unrealistic gel-like structures above the gel-to-fluid transition temperature.
We thank Senthil Kandasamy from D. E. Shaw Research for sharing valuable results from tests of the C36 parameter set. Financial support from the NIH (GM 72558 and 15101 to ADM and GM 86685 to DJT), NSF (CHE-0750175 to DJT) and the University of Maryland are appreciated. This research was supported in part by the Intramural Research Program of the NIH, National Heart, Lung and Blood Institute under the advisement of Bernard Brooks (JBK). This study utilized the high-performance computational capabilities at the National Institutes of Health, Bethesda, MD (CIT Biowulf and NHLBI LoBoS clusters) and the HPCC at UMD/OIT.
Supporting Information Available. Details on the methods from small molecule simulations; tables of partial atomic changes of model compounds and DPPC from quantum mechanical calculations; table of C36 dihedral parameters that differ from C27r; figures of adiabatic maps for selected torsions; figures of surface area vs. time, and deuterium order parameters for a 130-ns simulation of DPPC; figure comparing experimental multilayer and vesicle 13C NMR T1 for a DPPC bilayer from simulations with C36 and C27r. This information is available free of charge via the Internet at http://pubs.acs.org.