Lipid21: Complex Lipid Membrane Simulations with AMBER

1In recent years, the convergence of lipid membrane force field development and advances in computer hardware have allowed the molecular dynamics simulation of numerous membrane phenomena, such as bilayer phase transitions, vesicle dynamics, and the behavior of realistic cell membranes. Whereas coarse-grained force fields are approaching simulations of membrane patches that are of biologically relevant time scales and dimensions, such simulations are beyond the reach of atomistic force fields in the absence of specialized hardware. However, atomistic force fields come to the fore when studying processes in atomic detail, such as drug-membrane interactions, membrane protein behavior, or membrane-peptide partitioning. To gain the most relevant insights, atomistic lipid force fields must be parametrized appropriately, such that agreement to experiment is optimal.

2AMBER is both a molecular dynamics software suite and a series of interoperable force fields covering proteins, nucleic acids, carbohydrates and a selection of lipid types. Beginning with Lipid11, followed by GAFFlipid and then Lipid14 the AMBER lipid force field has gone through a number of iterations but until recently only supported simulation of phosphatidylcholine (PC), phosphatidylethanolamine (PE), and cholesterol lipids. Although Lipid14 allowed tensionless simulations of bilayers, deficiencies were reported, such as poor agreement with NMR headgroup order parameters, a mismatch of the DPPC phase transition temperature, and the observation of gel phase behavior of DMPC at long simulation times.

3In this work, we update the AMBER lipid force field to address these issues, while also extending the coverage of lipid species. Due to the modular nature of the force field, parameter derivation is only required for individual headgroup and tail units. Head group torsion parameters are simultaneously fitted to QM conformational energies. Partial charges are derived for all lipid head groups. Hydrocarbon chain parameters are updated, allowing better agreement with phase transition temperatures. We run multiple long validation simulations for each lipid species and find suitable agreement with biophysical measurements. Finally, model lipid raft simulations are performed to investigate the impact of increasing PUFA concentrations on cholesterol molecules.

4We also document instances when lipid simulations with AMBER give undesired results—simulations of anionic lipids are found to be sensitive to the cation type and force field model. We recommend best practices for lipid simulations with AMBER, which now cover a range of lipid types with the modular Lipid21 force field.

5The initial parameter set comprises lipid tail van der Waals parameters derived in the fitting of Lipid14 to reproduce densities and heats of vaporization of long chain hydrocarbons, and lipid tail partial charges, also derived in the fitting of Lipid14. The remaining parameters were initially covered by GAFF version 1.7 and the AMBER phosphate parameters.

6A summary of all parameters fitted in this work is as follows: the alkane C–C–C angle in lipid tails with atom types cD–cD–cD was found to have a large influence on bilayer fluidity and was fit to quantum mechanical energies (see Figure 1). To accommodate this fit, all hydrocarbon torsions in lipid tails were refit (atom types cD–cD–cD–cD, cD–cB–cB–cD, cB–cD–cD–cD, cB–cB–cD–cD, cB–cB–cD–cB) (see Supporting Information (SI) Figures S3–S7). All headgroup partial charges are derived at the MP2/cc-pVTZ level with the polarizable continuum model using the Lipid11 capping strategy, covering phosphatidylcholine (PC), phosphatidylglycerol (PG), phosphatidylserine (PS), phosphatidylethanolamine (PE), phosphatidic acid (PH−), and sphingomyelin (SPM) head groups. Lipid tail partial charges are adopted from Lipid14. Finally, headgroup torsion terms are fit to quantum mechanical energies (atom types cC–oS–cA–cA, oS–cA–cA–oS, cA–cA–cA–oT, cA–cA–oT–pA, cA–oT–pA–oT, oT–cA–cA–nA, cB–cA–cA–nN) (see Figure 2).

7All QM single-point energies, angle scans, and torsion scans were performed using Gaussian 09. Angle and torsion scans on model hydrocarbon molecules were performed at the MP2/cc-pVTZ level. Angles were scanned from 90 to 150° in 1° increments. Torsions were scanned over 360° in 10° increments. Single-point energy calculations on model glyceride, ceramide, and phosphatidylcholine molecules were performed at the MP2/cc-pVDZ level, with application of the polarizable continuum model to create an implicit solvent environment for phosphatidylcholine only. To obtain relevant glyceride, ceramide, and phosphatidylcholine conformations, POPC and PSM simulations were performed for 100 ns using initial Lipid21 parameters, and 1000 random lipid structures were extracted. Prior to the QM single-point energy calculation, each molecule was minimized using AMBER20 for 1000 steps, with the first 500 steps using the steepest descent and the final 500 steps using the conjugate gradient method, with restraints of 5000 kcal/mol/rad on each of the torsions being fitted. The GBn generalized Born model (igb = 7) was used during minimization of model phosphatidylcholine molecules only. We repeated this process on separate trajectories to create a test set, used for parameter validation. Parameters for the model glyceride, ceramide, and phosphatidylcholine molecules consisted of initial Lipid21 parameters and partial charges derived at the MP2/cc-pVDZ level after optimization of a single molecular conformation, allowing a two-stage RESP fit.

8All torsion parameters were fitted to the QM energies using the Pyevolve Genetic Algorithm. First, conformational energies were evaluated using AMBER20 with the torsional barrier heights set to zero for torsions of interest. When fitting to a 1D torsion scan, the model fragment was minimized for a small number of steps for each point along the scan, with a restraint applied to the torsion under study.

9The Genetic Algorithm was then used to identify torsion barrier height and phase shift parameters which minimize the least-squares fitness function between QM and MM energies:where N is the number of conformations, EQM is the reference energy from QM and EMM is the MM energy from AMBER, calculated as the “zeroed” torsion(s) plus the torsion energy contribution using fitted torsion parameters determined via the standard AMBER torsion equation:where the sum runs over periodicity n, Vn is the barrier height, θn is the torsion angle and γn is the phase shift. When fitting 1D torsion scans, the periodicity n = 1–5. When fitting single torsions to collections of glyceride and ceramide conformations, n = 1–3. When simultaneously fitting phosphatidylcholine torsions, periodicity was retained from the starting torsion parameters, with the exception of the cA–cA–oT–pA and cA–cA–cA–oT torsions, to which additional period terms were introduced. The phase shift γn was allowed to take terms 0 or 180° during fitting of hydrocarbon torsions. During fitting of headgroup torsions, phase shift γn was allowed to vary in 60° increments to accommodate headgroup stereochemistry.

10Partial charge fitting was performed for all lipid headgroup units using the Lipid11 capping strategy, at a higher level of QM theory. For each headgroup type, a 100 ns bilayer simulation was performed using Lipid21 fitted parameters, and 200 random lipids were extracted. The methyl headgroup capping procedure was applied as in Lipid11, and a single-point QM calculation was performed at MP2/cc-pVTZ with the polarizable continuum model. Partial charges were then derived for each of the head groups using a two-stage RESP fit; final charges are the average of the 200 individual fits. This procedure provides Boltzmann-averaged implicitly polarized charges over a conformational ensemble.

11Due to the different level of QM theory, methyl acetate capping group partial charges were rederived at the MP2/cc-pVTZ level with the polarizable continuum model; updated partial charges are available in the SI Figure S1. The capping procedure used for fitting of sphingomyelin charges is also available in the SI Figure S2.

12To prepare a starting structure, a DPPC bilayer was equilibrated at 300 K for 300 ns. The final coordinates were then used for repeat heating simulations, run in the NPT ensemble for 100 ns while increasing the temperature from 300 to 350 K, at a rate of 0.5 K/ns. Three independent heating runs were performed. The area per lipid was monitored for each simulation to obtain an estimate of the phase transition temperature.

13Twenty homogeneous lipid bilayer systems were simulated to validate the structural properties of the membranes, covering zwitterionic, anionic, PUFA, and sphingomyelin lipid species (see Table 1). Initial coordinates were obtained from the CHARMM-GUI. To ensure relaxation of each bilayer, a simulated annealing protocol was applied to the initial coordinates prior to equilibration and production runs. Specifically, bilayers were heated to 393 K in a 1 ns NPT simulation. Each membrane was then simulated at 393 K for 10 ns in the NPT ensemble. Systems were then cooled to the target temperature (Table 1) in a 1 ns NPT run and equilibrated for 100 ns. Following equilibration, three independent production runs of 300 ns were performed for each lipid species.

14All simulations performed in this work used a similar protocol. First, each system was minimized for 10000 steps, with the first 5000 steps using steepest descent and the final 5000 steps using the conjugate gradient method. Heating to 100 K was performed in a 5 ps NVT simulation with restraints of 5 kcal/mol/Å2 on nonsolvent molecules. The same restraints were maintained during a 100 ps NPT simulation and systems heated to the target temperature (Table 1). All restraints were removed and the PME was used to treat all electrostatic interactions with a real space cutoff of 10 Å; a long-range analytical dispersion correction was applied to the energy and pressure. A constant pressure of 1 atm was maintained using a semi-isotropic Berendsen barostat for lipid systems (x- and y- box vectors were coupled, with the z-vector allowed to change freely) and a pressure relaxation time of 1 ps. Temperature was controlled by the Langevin thermostat, with a collision frequency of γ = 1 ps–1. Water was modeled using the TIP3P water model and counterions using the Åqvist parameters. All simulations used AMBER20 PMEMD CUDA on GPU cards using the SPFP precision model., Single exploratory production runs of 300 ns were also executed with the Monte Carlo barostat for pressure coupling, keeping all other settings identical. Concerning single 1 μs simulations: settings were identical to those described, with the exception that hydrogen mass repartitioning was used to allow a 4 fs time step.

15Three systems with a raft-like composition were prepared and simulated, with increasing mole fraction of the PUFA lipid DAPC (see Table 2). These underwent heating and equilibration of 100 ns as previously described, followed by two repeat 1 μs production runs at 310 K and application of hydrogen mass repartitioning to allow a 4 fs time step.

16The area per lipid of each system was calculated directly from the simulation box size as done previously. Volume per lipid was determined in a similar manner, using the simulation box size and the average volume of a TIP3P water at the target temperature over a 50 ns NPT simulation. Bilayer thicknesses were determined from electron density profiles calculated with CPPTRAJ. NMR order parameters were also calculated with CPPTRAJ. X-ray and neutron scattering form factors were determined with the SIMtoEXP program. To study the cholesterol tilt angle distribution in the raft-like simulations, the z-position of the oxygen atom of each cholesterol was extracted using CPPTRAJ, as was the z-position position of the sterol carbon connecting the acyl chain of the cholesterol (see SI Figure S8). The angle between the bilayer normal and this vector was then calculated. The z-coordinate of cholesterol oxygen atoms was monitored for transit events, with any visit to −2 < z < 2 Å from the bilayer center being counted as an event.

17The C–C–C angle parameter was found to have an influence on bilayer fluidity (Lipid21 atom types cD–cD–cD). Therefore, we refitted this angle parameter to high level gas-phase QM scans. The updated angle parameters are seen to better capture the energy minima about 112° and the parabolic energy well at higher angle stretches (see Figure 1).

18To accommodate the updated C–C–C angle parameter, hydrocarbon torsions were also refit at the same level of QM theory—these fits are detailed in the SI.

19Glycerol, ceramide, and phosphatidylcholine headgroup torsions were refit to better reproduce QM single-point energies (see Figure 2). To obtain relevant conformations, 1000 lipids were randomly extracted from 100 ns runs of POPC and PSM bilayers. Each molecule was then minimized using AMBER with restraints on the torsion(s) of interest. The QM energy was then evaluated at the MP2/cc-pVDZ level (with the PCM model applied to phosphatidylcholine molecules only, to mimic the solvent environment lipid head groups experience). Torsions of interest were then fitted simultaneously using the Genetic Algorithm protocol detailed in the Methods section. Regarding the quality of the torsion fitting, glycerol energies could be corrected to reproduce QM energies well, as seen from both the training and test sets. The phosphatidylcholine results are only a marginal improvement, with little change in RMSE, although there is a shift in correlation from r2 = 0.03 to r2 = 0.31 on the training set and r2 = 0.02 to r2 = 0.18 on the test set, potentially due to the complex electrostatic nature of these zwitter-ionic model molecules. Although it was possible to obtain better reproduction of QM energies by introducing additional periodic terms to each torsion, initial bilayer simulations were poor, indicating the more complex torsion sets suffered from overfitting. Finally, sphingomyelin was not covered previously; however, the torsion fit results indicate that the ceramide torsion reproduces the QM conformational energy landscape well. It was found that this torsion fitting procedure not only improved the bulk phase properties of lipid bilayer simulations but also headgroup order parameter agreement with NMR measurements, similar to recent Slipids results.

20As can be seen from Table 3, simulations with Lipid21 capture well the bulk phase properties of different lipid species as determined by experiment. The performance of PC lipid types is similar to Lipid14, while the area per lipid of POPE comes into good agreement with the experimental result, which was previously under-predicted with Lipid14. Regarding the charged lipids, areas per lipid for all PG and PS lipid types agree well with experiment. However, area per lipid for both DSPG and DSPC are slightly low, suggesting that long acyl chains cause slight condensing of membranes with Lipid21.

21Simulations of PUFA lipids (SDPC, DAPC) are close to available experimental data for area per lipid and bilayer thickness of SDPC; the SDPC area per lipid is approximately 3 Å2 below experiment. The area per lipid of DAPC is similar to the Charmm36 parametrization results. Finally, the results for sphingomyelin are encouraging, although the area per lipid is 4–5 Å2 below recently determined experimental values, suggesting there is room for further improvement.

22All production simulations used the Berendsen barostat for pressure coupling. We also ran single exploratory simulations of 300 ns per lipid type using the Monte Carlo barostat, which allows a speed-up in simulation throughput, keeping all other simulation settings identical. As shown in Table S1, we find that areas per lipid are depressed for all lipid types, dropping by an average of 3.78 Å2. We therefore recommend the Berendsen barostat pressure coupling for lipid bilayer simulations in AMBER.

23In the next sections, lipid ordering is analyzed, providing an additional comparison to experiment.

24One of the motivations to reparameterize the PC lipid model was poor agreement of headgroup order parameters with available NMR data, since the headgroup order parameters provide an accurate experimental picture of bilayer structure.

25Calculated headgroup order parameters for POPC, DPPC, POPS, and POPG from Lipid21 simulations are shown in Figure 3, along with comparison to experiment and in the case of POPC and DPPC, results for previous Lipid14 simulations. Where experimental data are available, POPC and DPPC results agree very well with experiment, showing a marked improvement over Lipid14. Concerning the charged lipids POPS and POPG, headgroup order parameters trend with experiment yet leave room for improvement. In particular, order parameters of the α C–H vector in POPS are significantly far from experiment. However, the overall results indicate that bilayer structure agrees suitably with the NMR experiments.

26The ordering of the lipid chains can also be calculated allowing comparison to NMR experiments. Figure 4 details the results for POPC and DPPC lipids, finding good agreement with experiment, although the splitting of the sn-2 carbon-2 position does not reach the level determined experimentally and the POPC and DPPC carbon chains are too ordered near the headgroup region.

27A similar analysis can be performed for the PSM sphingomyelin model, for which NMR data are available for the N-linked chain.Figure 5 details the Lipid21 results, finding good agreement with experiment for the N-linked chain, with the only exception being the carbon-2 position, which does not show the level of splitting observed in NMR—a similar behavior to the Lipid21 POPC and DPPC models. However, the suitable agreement indicates the bilayer structure of PSM is close to that studied in NMR experiments.

28Finally, small-angle X-ray scattering form factors (SAXS) provide another detailed reference of lipid bilayer structure. The results for POPC, DPPC, POPS, and POPG simulated with Lipid21 and comparison to experimental data are given in Figure 6. For all models, agreement with experiment is suitable, in particular the positions of the minima. A mismatch between experimental and simulated SAXS profiles is however observed around the second maxima, being particularly pronounced for POPS.

29Also shown are the SAXS and small-angle neutron scattering (SANS) form factors for PSM collected at 318 K and 100% D2O, with comparison to Lipid21 simulations (328 K) in Figure 7. Both scattering profiles compare well to experiment, with the SAXS results showing minima at equivalent positions as experiment. As with the POPS result, the second SAXS maxima is overpredicted in comparison to experiment.

30Following updating of hydrocarbon bonded parameters, fitting of headgroup torsions to single-point QM energies and partial charge derivation, the initial Lipid21 parameter set was tested for the ability to reproduce the DPPC experimental melting point (see Figure 8). A lipid force field comparison from Pluháčková et al. indicated that the Lipid14 DPPC melting point is too high, estimating a value of 343 K, far above the experimental DPPC melting point of 314 K., Melting point scans were performed with a similar methodology using the initial Lipid21 model, in which bilayers are gradually heated from 300 to 350 K, at a rate of 0.5 K/ns. Plots of the area per lipid as a function of temperature show an inflection point at approximately 335 K for the initial Lipid21 parameter set. In order to better tune the force field to bilayer melting points, we scale down the 1–4 Lennard-Jones interaction of the acyl chain torsion parameter cD–cD–cD–cD, from the standard AMBER setting of 0.5 scaling to 0.167 (via the SCNB setting). By default, the Amber force field scales down 1–4 nonbonded terms to reduce the exaggeration of short-term repulsion caused by the Lennard-Jones potential (to model van der Waals interactions) and Coulomb potential (to model electrostatic interactions). The 1–4 scaling factor was previously proposed as a parameter by which to tune the melting point and other liquid phase properties of N-alkanes. The Lennard-Jones 1–4 scaling of 0.167 was found to bring the DPPC melting point down to 319 K, much closer to the experimental value of 314 K. As such, this scaling factor, which is included automatically when building the topology file, is incorporated into the final Lipid21 parameter set. It should be noted that this bilayer heating method is expected to overestimate the true melting point temperature by 5–10 K.

31Validation simulations for the anionic lipids PG, PS, and PH– used potassium counterions modeled with the Åqvist (ff99) force field to charge neutralize the simulation box, as required by PME to treat long-range electrostatics. The high concentration of cations in the water layer is likely physically unrealistic; further, the predominant extracellular cation is sodium, with potassium found with higher concentration inside the cell. The ff99 cation model was previously found to be defective, leading to formation of unphysical salt crystals when modeling high salt concentrations; Joung and Cheatham therefore considerably revised ion parameters for AMBER simulations. However, during validation runs of anionic lipids using either K+, Na+, or Ca2+ counterions and JC parameters, bilayers were observed to condense and adopt an area per lipid below experiment (see Figure 9). In fact, even Na+ with ff99 resulted in similar behavior. We therefore used K+ counterions and ff99 parameters in the present work.

32Unphysical lipid-ion pairing behavior is a known artifact in simulations of charged lipid bilayers and can result in condensing of the membrane and structural headgroup effects. Ion binding to lipid head groups has been an active area of research both computationally and experimentally, with the current view being that monovalent cations do not bind specifically to PC head groups at submolar concentrations (except Li+), whereas multivalent ions (such as Ca2+) do show specific binding. In the case of anionic head groups such as PS, monovalent ions show only weak binding, while Li+ and multivalent ions are able to form strong dehydrated complexes.

33Correct ion pairing behavior has been difficult to capture in molecular dynamics simulations, with a number of methods implemented to improve the description of lipid–ion interactions. Venable et al. derived pair-specific ion parameters for sodium, modifying the van der Waals interaction between Na+ and lipid headgroup oxygen atoms, allowing better agreement with NMR experiments for PG and PS simulations. Melcr et al. modified the Lipid14 POPC model to better capture cation interactions using the “electronic continuum correction” method, which involves scaling of the cation charge and lipid headgroup partial charges to account for electronic polarizability. This approach resulted in better agreement with NMR headgroup order parameter changes upon the addition of Na+ or Ca2+ ions (the “molecular electrometer” concept).

34Such reparameterization is evidentially required for Lipid21 but is beyond the scope of the present work. Rather, we recommend K+ cations with ff99 parameters for simulations of entirely anionic lipid types. For typical simulations, bilayers are likely to be predominantly composed of PC lipids with lower fractions of anionic lipid types. Anticipating this, we ran μs simulations of POPC and DMPC with 0.15 M NaCl using JC parameters, as are typically paired with AMBER protein or DNA force fields. As demonstrated in Figure 10, this does not significantly alter the bilayer structure, unlike the effect observed for anionic lipids.

35The ability to model diverse lipid species with Lipid21 enables simulations of complex membranes of different lipid compositions. Cholesterol has a known aversion to poly unsaturated fatty acid lipids such as DAPC, adopting a larger tilt angle in DAPC bilayers than DPPC. Molecular dynamics (MD) studies have also found cholesterol to display a greater flip-flop rate in bilayers doped with PUFA lipids. In the extreme, cholesterol in pure DAPC bilayers shows a preference to reside at the membrane interior, as studied with neutron diffraction and MD.

36We constructed and simulated three raft-like membranes, containing POPC, PSM, cholesterol, and increasing fractions of the PUFA lipid DAPC (see Table 2) with Lipid21 to investigate cholesterol tilt angle and transit events as a function of DAPC doping. These underwent heating and equilibration of 100 ns at 310 K followed by two repeat production runs of 1 μs. Figure 11 plots the most probable cholesterol tilt angle, taken over all 40 cholesterol copies and two repeat 1 μs simulations, as a function of increasing DAPC doping. Also shown are the total number of cholesterol transit events over a combined 2 μs of simulation as a function of increasing DAPC, with a transit event being considered as any visit of cholesterol oxygen atoms to the bilayer center.

37It can be observed from Figure 11 that the most probable cholesterol tilt angle in 0 mol % and 30 mol % DAPC raft-like bilayers is very similar, approximately 13°. At 60 mol % DAPC, this increases to approximately 14°. The raw probability distribution plots (SI Figure S9) reveal a broadening of the probability density, indicating destabilization of cholesterol molecules with increasing DAPC doping. Figure 11 also depicts the total number of cholesterol transit events as a function of DAPC doping. In 0 mol % DAPC bilayers, cholesterol molecules are very stable, displaying no transit events. At 30 mol % DAPC, there are two transit events. Analysis of the raw z-coordinate plots (SI Figure S10) show two complete flip-flop events, one from each 1 μs repeat simulation. Finally, at 60 mol % DAPC, cholesterol molecules are further destabilized, showing a total of five transit events, a mixture of complete flip-flop events and transient visits to the bilayer center. The increase in both tilt angle and transit events indicate that with increasing mol % of DAPC, Lipid21 qualitatively reproduces the observation of the perturbing effect of DAPC on cholesterol molecules.

38The current work extends and improves the AMBER lipid force field, increasing coverage to anionic lipids, PUFA lipids, and sphingomyelin. The modular nature of the force field allows easy expansion to other lipid species. The PC model is improved; in particular, headgroup NMR order parameters are better reproduced in comparison to Lipid14. Bulk bilayer structural properties find good comparison to experiment for all lipids studied, matching well areas per lipid, bilayer thickness, NMR order parameters and X-ray scattering profiles. Simulations of raft-like membranes demonstrate an expected perturbation of cholesterol molecules as a function of increasing PUFA content.

39The parametrization strategy required that the hydrocarbon parameters be updated, such that the DPPC melting point is better captured. Lipid headgroup partial charges are derived at a higher level of QM theory, including use of the polarizable continuum model allowing implicit inclusion of polarization. Lipid headgroup torsion parameters are updated to better reproduce QM energies.

40Further work does remain concerning membrane simulations in AMBER. Although results for sphingomyelin reproduce available experimental data reasonably well, the area per lipid remains 4–5 Å too low. Head group NMR order parameters are improved for PC lipids; however, PG and PS do not find such a close match to experimental data. The type of cation used to charge neutralize anionic bilayer systems and the force field used to model ions is found to have a dramatic effect, artificially condensing POPS and POPG systems unless potassium counterions with Åqvist parameters are used. Correctly capturing the interaction of cations with different lipids species is an active area of research, with a number of methods proposed. Clearly, such work is required to update Lipid21. However, for simulations of physiologically relevant membranes, typically consisting of PC lipids and smaller fractions of other species, we recommend NaCl/KCl ions using JC parameters, as are typically combined with AMBER protein or DNA force fields. Indeed, 1 μs simulations of POPC and DMPC with 0.15 M NaCl using JC parameters reproduce experimental areas per lipid. Additionally, we identify that the Monte Carlo barostat results in the depression of areas per lipid for all lipid types studied in this work (see SI Table S1) and recommend the Berendsen barostat for pressure coupling for lipid bilayer simulations with AMBER.

41Lipid21 significantly advances lipid membrane simulations with AMBER. The modular nature of the force field makes further lipid coverage achievable with minimal parametrization. We find suitable comparison to experiment for the lipids investigated. We also identify areas for improvement and recommend best practices for membrane simulations using Lipid21. The parameter set is available in the Supporting Information or can be downloaded from https://github.com/callumjd/lipid21.