Lipid14: The Amber Lipid Force Field.

The AMBER lipid force field has been updated to create Lipid14, allowing tensionless simulation of a number of lipid types with the AMBER MD package. The modular nature of this force field allows numerous combinations of head and tail groups to create different lipid types, enabling the easy insertion of new lipid species. The Lennard-Jones and torsion parameters of both the head and tail groups have been revised and updated partial charges calculated. The force field has been validated by simulating bilayers of six different lipid types for a total of 0.5 μs each without applying a surface tension; with favorable comparison to experiment for properties such as area per lipid, volume per lipid, bilayer thickness, NMR order parameters, scattering data, and lipid lateral diffusion. As the derivation of this force field is consistent with the AMBER development philosophy, Lipid14 is compatible with the AMBER protein, nucleic acid, carbohydrate, and small molecule force fields.

Introduction

Membranes are integral components of the cell, separating intracellular compartments from the cytosol. Such membranes consist of a back-to-back arrangement of n class="Chemical">lipid molecules, driven into a bilayer structure by the hydrophobic effect, leaving the polar class="Chemical">pan class="Chemical">lipid head groups exposed to water, and bringing the nonpolar lipid tail groups together. The composition of cell membranes is complex, with constituent species including, but not limited to, saturated and unsaturated PC and PE lipids, sphingomyelin and cholesterol, which serve as a matrix in which membrane proteins may reside.[1] Cell membranes possess functions such as regulating transport in to and out of the cell and modulating the activity of membrane embedded ion channels and proteins.[2,3]

In order to probe the many roles of membranes in the cell, membrane structures are studied experimentally using techniques such as X-ray and neutron scattering, IR/Raman, and NMR spectroscopy.[4,5] To gain atomic-level resolution, however, membranes may also be simulated computationally using molecular dynamics (MD). The validity of results obtained using MD methods depends, to a large extent, on the potential energy function, or force field, that is used.

Membranes can be simulated using all-atom, united-atom, or coarse-grained models,[6−13] with the increasing simplicity of each representation allowing access to larger models and longer time scales, at the expense of atomic detail. All-atom models may be preferred for bilayer simulations due to their ability to reproduce NMR order parameters and easy combination with all-atom protein, nucleic acid, n class="Chemical">carbohydrate, and small molecule force fields.[14−16] One such MD software class="Chemical">package that includes all-atom simulations is AMBER,[17,18] which contains extensively validated protein, nucleic, class="Chemical">pan class="Chemical">carbohydrate, and small molecule force fields. The AMBER simulation code has also been ported to NVIDIA GPU cards, making simulation speeds in excess of 100 ns per day possible for systems of 25,000 to 50,000 atoms, with systems of up to 4 million atoms possible on the latest generation hardware.[19−21] Although the AMBER force field suite includes numerous different species of biological interest, parameters for the simulation of lipids have traditionally been lacking. Recently the Lipid11 framework was introduced, which is a modular lipid AMBER force field, allowing the simulation of a number of lipids via the combination of different head and tail groups.[22] Lipid11 used force field parameters predominantly taken directly from the General Amber Force Field (GAFF).[16] Although previous studies found some success in simulating lipid bilayers using GAFF parameters,[23−25] longer time scale simulations required a surface tension term in order to keep a bilayer in the correct phase at a given temperature.[22]

In this paper we draw inspiration from previous work[26] to update n class="Chemical">Lipid11 headgroup and tail group charges and class="Chemical">parameters to enable proper tensionless simulation of class="Chemical">pan class="Chemical">lipid bilayers, thereby creating a modular AMBER lipid force field. Lennard-Jones and torsion parameters are revised, while charges are derived according to the standard AMBER convention, as was implemented in Lipid11 with minor improvements in sampling for Lipid14. As such, the resulting parameters are expected to be compatible with other force fields in the AMBER package.

We validate these parameters for a number of PC and PE n class="Chemical">lipids – six different class="Chemical">pan class="Chemical">lipid types are simulated for a total of 500 ns each, with resulting area per lipid, volume per lipid, isothermal compressibility, NMR order parameters, scattering form factors, and lipid lateral diffusion finding favorable comparison to experiment, particularly for PC lipids. We also assess the conformation of the lipid tails by examining the number of rotamers and rotamer sequences. To fully test the reproducibility of the results, a number of additional GPU and CPU runs were performed. We believe that these parameters will be transferable thus allowing the easy insertion of new lipid types into Lipid14 owing to the modular nature of the force field. The Lipid14 force field will be released with the upcoming AmberTools14. It is our intention to incorporate support for additional lipid types and other residues commonly found in bilayers, as part of the upcoming release. The derivation of these parameters will be described elsewhere.

Parameterization Strategy

Generation of Lipid14 Parameters

n class="Chemical">Lipid14 aims to extend the class="Chemical">pan class="Chemical">Lipid11[22] framework to allow accurate, tensionless simulation of numerous lipid types via a modular lipid force field. Bond, angle, torsion, and Lennard-Jones (LJ) parameters in Lipid11 were taken directly from the General Amber Force Field (GAFF).[16] Although bond and angle parameters should not require updating, it was expected that torsion and LJ parameters would need modification to realize this aim.

Previous work on n class="Chemical">lipid simulation with AMBER indicated that the Lennard-Jones and torsion class="Chemical">parameters of the class="Chemical">pan class="Chemical">lipid aliphatic tails are intricately linked to the ability of the force field to reproduce experimentally observed bilayer properties.[26] Indeed, the majority of all-atom lipid force field parameterization has involved the modification of the LJ and torsion parameters used to model the aliphatic tail regions of the lipids.[7,8,10,27]

Simulating a box of 144 n class="Chemical">pentadecane chains using the standard GAFF LJ and torsion class="Chemical">parameters at constant pressure (1 atm) and temperature (298.15 K) causes the class="Chemical">pan class="Chemical">hydrocarbon chains to ‘freeze’ into a crystalline state in under 2 ns (see Figure 1). Such a scenario has previously been observed by Klauda et al. using a box of tetradecane molecules and the AMBER99 force field.[28] Consequently, the calculated density and heat of vaporization are much higher and the diffusion much lower than experiment.

Figure 1

A box of 144 pentadecane molecules simulated in the NPT ensemble at 298.15 K using the General Amber Force Field[16] to model the carbon chains.

A box of 144 pan class="Chemical">pentadecane molecules simulated in the NPT ensemble at 298.15 K using the General Amber Force Field[16] to model the class="Chemical">pan class="Chemical">carbon chains.

In light of this, LJ and torsion parameters were modified to reproduce the experimental density (ρ) and heat of vaporization (ΔHvap) of n class="Chemical">alkanes covering a range of chain lengths. Given that both the torsion and LJ class="Chemical">parameters affected the simulated ρ and ΔHvap, these class="Chemical">parameters were altered simultaneously, with the CH2–CH2–CH2–CH2 torsion being fitted to ab initio data using class="Chemical">pan class="Chemical">Paramfit.[17,29] Satisfactory agreement was found by modifying only the LJ parameters of the hydrogen atom on the alkane chains.

The modified LJ and torsion parameters were tested by examining the density ρ, heat of vaporization ΔHvap, the diffusion D, the pan class="Chemical">13C NMR T1 relaxation times, and the trans/gauche conformer populations of a selection of class="Chemical">pan class="Chemical">hydrocarbon chains.

The parameters of the n class="Chemical">ester linkage region connecting headgroup and tail residues in the class="Chemical">pan class="Chemical">lipids were then examined using methyl acetate as a model compound. The density and heat of vaporization of methyl acetate were calculated and ΔHvap found to be in poor agreement with experiment. Hence the Lennard-Jones parameters of a number of atoms in this region were also adjusted until better agreement with experiment was achieved.

Once optimal n class="Chemical">hydrocarbon and class="Chemical">pan class="Chemical">glycerol parameters were identified, attention turned to the lipid partial charges. The Lipid11 force field used a capping procedure, separating the lipid head and tail groups into ‘residues’, thus creating a modular lipid force field.[22] The standard AMBER RESP protocol[30] was used to generate partial charges from quantum mechanical (QM) optimized structures, using six different orientations of a single conformation.

This methodology was kept for pan class="Chemical">Lipid14. However, in line with other all-atom class="Chemical">pan class="Chemical">lipid force field parameterizations,[10,31] a greater number of conformations were used per residue (twenty-five headgroup conformations, fifty tail conformations), with the partial charges calculated as an average over all conformations. The head and tail group starting structures were extracted from previous in-house bilayer simulations. This allows one to obtain Boltzmann weighted charges, introducing a temperature dependence.[10,31] The electrostatic potential (ESP) was calculated directly from the conformations extracted from a bilayer simulation, with no QM optimization performed on those structures. Charges were derived using the standard AMBER RESP protocol (at the HF/6-31G* level in gas phase).[30]

Finally, on identification of appropriate LJ parameters and calculation of n class="Chemical">lipid class="Chemical">partial charges, torsions involving the class="Chemical">pan class="Chemical">ester linkage in the glycerol region were fitted to QM scans performed on a capped lauroyl (LA) tail moiety.

Hydrocarbon Tail Parameters

The n class="Chemical">alkane CH2–CH2–CH2–CH2 torsion potential was refitted using torsion scans performed on class="Chemical">pan class="Chemical">hexane and octane molecules using an estimation of the CCSD(T)/cc-pVQZ level of theory via the HM-IE relation[32]where the small basis set (SBS) was cc-pVDZ and the large basis set (LBS) was cc-pVQZ. Consequently molecules were optimized at the MP2/cc-pVDZ level before single-point energy calculations were performed at the CCSD(T)/cc-pVDZ and MP2/cc-pVQZ levels. In order to obtain a representative torsion fit, it was ensured that torsion scans, conducted at increments of 15°, included local minima of hexane and octane.[10,28,33] The Paramfit program of AmberTools13[17,29] was used to perform a multiple molecule weighted torsion fit, with the tgt local minima of hexane and tgttt local minima of octane given a weighting of 10, all other local minima given a weighting of 4 and cis conformers given a weighting of 0.1; all other structures were assigned a weighting of 1. These weighting values have previously given good results for alkane torsion fitting.[33] Torsions were fitted using a genetic algorithm.

Lennard-Jones parameters were modified while simultaneously refitting the torsion parameters (see Table 1) until good agreement with experiment was achieved for heats of vaporization and densities for a range of n class="Chemical">alkane chains. These properties were monitored by performing liquid phase simulations of class="Chemical">pan class="Chemical">pentane, hexane, heptane, octane, decane, tridecane, and pentadecane. The alkane trajectories also enabled the calculation of a number of other properties for comparison to experiment.

Table 1

Modification of LJ and Torsion Parameters of Alkane Chains

	LJ parameters			CH2–CH2–CH2–CH2 torsion
	atom type	radius R (Å)	well-depth ε (kcal mol–1)	force constant PK (kcal mol–1)	periodicity PN	phase (deg)
Lipid11	cA	1.9080	0.1094	0.20	1	180
	hA	1.4870	0.0157	0.25	2	180
				0.18	3	0
Lipid14	cD	1.9080	0.1094	0.3112	1	180
	hL	1.4600	0.0100	–0.1233	2	180
				0.1149	3	0
				–0.2199	4	0
				0.2170	5	0

Initial charges for each n class="Chemical">hydrocarbon chain were generated using the standard AMBER RESP protocol (optimization and calculation of the ESP at HF/6-31G* level in gas phase). A box of 288 (class="Chemical">pan class="Chemical">pentane, hexane and heptane) or 144 (octane, decane, tridecane, and pentadecane) molecules were then simulated in the liquid phase at 298.15 K for 10 ns using updated LJ and torsion parameters. From each liquid phase simulation, fifty different chains were extracted and used for the charge calculation, with the ESP calculated directly from each structure at the HF/6-31G* level, and partial charges derived using the RESP fitting procedure. Charges were taken as an average over all fifty RESP fits.

The heat of vaporization was calculated according to[34]where Epot(g) is the average molecular potential energy in the gas phase, Epot(l) is the average molecular potential energy in the liquid phase, R is the gas constant, and T is temperature. In order to calculate the heat of vaporization, two simulations were performed in a similar manner to previous work.[26] A gas phase simulation consisting of a single pan class="Chemical">alkane molecule was run for 10 ns equilibration and 50 ns production under the NPT ensemble, at a temperature of 298.15 K to obtain Epot(). A liquid phase simulation consisting of a box of either 144 or 288 chains, run under the NPT ensemble with periodic boundary conditions using class="Chemical">particle mesh Ewald to treat long-range electrostatics[35] and a real sclass="Chemical">pace cutoff of 10 Å, was also performed for 60 ns with the first 10 ns removed for equilibration to obtain Epot(l). The temperature was maintained at 298.15 K using Langevin dynamics and a collision frequency of 5 ps–1. Bonds involving class="Chemical">pan class="Chemical">hydrogen were constrained using the SHAKE algorithm.[36] Pressure regulation was achieved with isotropic position scaling, a Berendsen barostat,[37] and a pressure relaxation time of 1 ps. The system was heated from 0 to 298.15 K over 20 ps, with a force constant of 20 kcal/mol/Å2 restraining the chains. This restraint was gradually decreased to 10, 5 and finally 1 kcal/mol/Å2 and the system simulated for 20 ps at each value of the force constant. The heat of vaporization was then calculated using eq 1, with results reported as block averages ± standard deviation using five equal blocks of 10 ns.

By reducing the n class="Disease">van der Waals radius (R) and well-depth (ε) of the class="Chemical">pan class="Chemical">alkane hydrogen atom type while simultaneously correcting the torsion fit, satisfactory agreement with experiment was achieved for ρ and ΔHvap for the alkane chains under study.

Three torsion scans about the C=C double bond were also performed on a n class="Chemical">cis-5-decene molecule using the MP2:CC extrapolation method – see Figure 2. When fitting the C=C double bond torsions, the LJ class="Chemical">parameters of the class="Chemical">pan class="Chemical">alkene hydrogen atom (Lipid14 atom type hB) were scaled until satisfactory agreement with experiment was found for ρ and ΔHvap of cis-2-hexene, cis-5-decene, and cis-7-pentadecene. Reducing both the van der Waals radius (R) from 1.459 to 1.25 and the well-depth (ε) from 0.015 to 0.007 resulted in far better agreement for the density values (ρ); however, it proved difficult to correct the heat of vaporization (ΔHvap) of cis-5-decene to that of the experimental value by only modifying the LJ parameters. This is due to the high charge on the double bond atoms. A similar problem has previously been encountered by Chiu et al.[13] and Jämbeck et al.[10]

Figure 2

The energy profile for rotating about selected torsions of a cis-5-decene molecule. Energy evaluated using QM and the HM-IE method (filled triangle ▲), AMBER with standard GAFF parameters (dotted line), and AMBER with Lipid14 parameters (black line). Torsion fits from the top are as follows: CH2–CH–CH–CH2, CH–CH–CH2–CH2, and CH–CH2–CH2–CH2.

The energy profile for rotating about selected torsions of a n class="Chemical">cis-5-decene molecule. Energy evaluated using QM and the HM-IE method (filled triangle ▲), AMBER with standard GAFF class="Chemical">parameters (dotted line), and AMBER with class="Chemical">pan class="Chemical">Lipid14 parameters (black line). Torsion fits from the top are as follows: CH2–CH–CH–CH2, CH–CH–CH2–CH2, and CH–CH2–CH2–CH2.

The uncorrected diffusion n class="Chemical">DPBC was calculated for each class="Chemical">pan class="Chemical">alkane using the slope of a mean-square displacement (MSD) plot versus time for the centers of mass, averaged over the trajectories of each molecule via the Einstein relationwhere nf is the number of dimensions (in this instance nf = 3), and Δr(t)2 is the distance that the molecule travels in time t. It is then possible to correct for the system size dependence of a diffusion coefficient calculated under periodic boundaries (DPBC) to yield the corrected diffusion coefficient Dcorr by adding the correction term derived by Yeh and Hummer[38]where kB is Boltzmann’s constant, T is the temperature, ε = 2.837297, η is the viscosity, and L is the length of the simulation box.

The diffusion was calculated from 100 ns NVE ensemble simulations (extended from the 50 ns NPT runs). PME was used, along with a 10 Å cutoff, at a temperature of 298.15 K. In order to avoid energy and temperature drift, it was necessary to remove the center of mass motion every 500 steps (nscm = 500), make both the shake tolerance and Ewald direct sum tolerance more stringent, and reduce the time step to 1 fs. Diffusion values were then calculated by taking the slope of the linear 2–5 ns region of the pan class="Disease">MSD versus time curve and the correction calculated using experimental viscosity values. Diffusion results are reported with standard deviations calculated by block averaging, with each 100 ns run divided into five 20 ns blocks.

The trans, gauche, pan class="Disease">end gauche, double gauche, and kinked gtg′+gtg conformer populations were evaluated from the 50 ns NPT runs by classifying torsion angles as either gauche plus (g+) 0–120°, trans (t) 120–240°, or gauche minus (g–) 240–360°.

pan class="Chemical">13C NMR T1 relaxation times were calculated from the NPT class="Chemical">pan class="Chemical">alkane simulations using the following formula for dipolar relaxation from the reorientation correlation functions of the CH vectors:[39]

This assumes motional narrowing and an effective C–H bond length of 1.117 Å,[40] with N specifying the number of protons bonded to the n class="Chemical">carbon and μ̂ the CH vector. T1 relaxation times were calculated from simulations of class="Chemical">pan class="Chemical">alkanes of four different chain lengths. These simulations were repeated using the same NPT alkane protocol at the experimentally relevant temperature of 312 K and production runs extended to 200 ns to improve sampling.

Head Group Parameters

It was found that the heat of vaporization of n class="Chemical">methyl acetate, a model compound representing the class="Chemical">pan class="Chemical">ester linkage region connecting the lipid head and tail groups, was not sufficiently close to experiment using GAFF/Lipid11 parameters. In order to correct for this discrepancy with experiment, the Lennard-Jones well-depths of the carbonyl oxygen (oC), carbonyl carbon (cC), and ester oxygen (oS) atoms were scaled until better agreement with experiment was obtained (see Table 2). The ester oxygen parameters were also applied to oxygen atoms in the phosphate region.

Table 2

Thermodynamic Properties of Methyl Acetate Simulated with GAFF/Lipid11 and Lipid14 and Comparison to Experimenta

	LJ parameters
	atom type	radius R (Å)	well-depth ε (kcal mol–1)	ΔHvap (kJ mol-1)	ρ (kg m-3)
Lipid11	oC	1.6612	0.210	39.11 ± 0.04	928.38 ± 0.09
	oS	1.6837	0.170
	cC	1.9080	0.086
Lipid14	oC	1.6500	0.140	33.0 ± 0.07	925.8 ± 0.05
	oS	1.6500	0.120
	cC	1.9080	0.070
Expt				32.29[41]	934.2[41]

All values at 298.15 K.

The density and heat of vaporization of n class="Chemical">methyl acetate were calculated using an identical procedure as for the class="Chemical">pan class="Chemical">hydrocarbon chains, with a box of 288 methyl acetate molecules applied for the liquid phase calculation. Simulations were run for 60 ns and the final 50 ns used for sampling. Results are reported as block averages (five blocks of 10 ns).

Partial Charges

RESP fitting was performed in exactly the same manner to pan class="Chemical">Lipid11,[22] using the capping groups shown in Figure 3. However a greater number of conformations were used to calculate the average charges for each unit.

Figure 3

Structure and charges of Lipid11/Lipid14 headgroup and tail group caps.[22]

Twenty-five n class="Chemical">phosphatidylcholine (PC) and twenty-five class="Chemical">pan class="Chemical">phosphatidylethanolamine (PE) head groups were extracted from previous bilayer simulations of DOPC and POPE; while fifty lauroyl (LA), myristoyl (MY), palmitoyl (PA), and oleoyl (OL) tails were extracted from bilayer simulations of DLPC, DMPC, DPPC, and DOPC, respectively. Each headgroup was then capped with a methyl group and each tail capped with an acetate moiety (Figure 3), according to the Lipid11 charge derivation methodology.[22] For each conformation, the ESP was calculated directly from each structure at the HF/6-31G* level using Gaussian 09.[42] Charges were taken as an average over all conformations for each residue. Resulting charges for the PC and PE headgroup residues and the LA, MY, PA, and OL tail group residues are detailed in the Supporting Information.

Structure and charges of pan class="Chemical">Lipid11/class="Chemical">pan class="Chemical">Lipid14 headgroup and tail group caps.[22]

Head Group Torsion Fits

Two torsions involving the n class="Chemical">ester linkage region (see Figure 4) were fitted to QM data. The scans were performed on a capped lauroyl (LA) tail residue, at 15° increments using the HM-IE method with Gaussian 09.[42] These were then fitted using class="Chemical">pan class="Chemical">Paramfit[17,29] for periodicity n = 1 to n = 5 using the genetic algorithm implemented in Paramfit.

Figure 4

A capped lauroyl tail group residue was used to fit the oS-cC-cD-cD and oC-cC-cD-cD torsions.

A capped lauroyl tail group residue was used to fit the pan class="Chemical">oS-cC-cD-cD and class="Chemical">pan class="Chemical">oC-cC-cD-cD torsions.

The energy profiles for rotating about selected torsions of a capped lauroyl tail group residue. Energy evaluated using QM and the HM-IE method (filled triangle ▲), AMBER with standard GAFF/n class="Chemical">Lipid11 class="Chemical">parameters (dotted line), and AMBER with class="Chemical">pan class="Chemical">Lipid14 parameters (black line). Torsion fits from the top are oC-cC-cD-cD and oS-cC-cD-cD.

As can be seen from Figure 5 n class="Chemical">Paramfit brings the class="Chemical">pan class="Chemical">oS-cC-cD-cD and oC-cC-cD-cD torsions into substantially better agreement with the QM data.

Figure 5

The energy profiles for rotating about selected torsions of a capped lauroyl tail group residue. Energy evaluated using QM and the HM-IE method (filled triangle ▲), AMBER with standard GAFF/Lipid11 parameters (dotted line), and AMBER with Lipid14 parameters (black line). Torsion fits from the top are oC-cC-cD-cD and oS-cC-cD-cD.

Parameterization

Hydrocarbon Parameters

The results for the n class="Chemical">alkane properties calculated using the updated torsion and LJ class="Chemical">parameters are shown in Table 3. The heat of vaporization values increase with class="Chemical">pan class="Chemical">alkane chain length, following the experimental trend and converging with experiment as the number of carbon atoms in the chain increases. The simulation values match experiment with an RMS error of 7.67%. The simulated densities are reproduced somewhat better than Δvap, with an RMS error of 2.60% when compared to experiment. The diffusion values lie close to experiment and decrease with increasing chain length, following the experimental trend; however, the RMS error between simulation and experiment remains significant at 20.92%. As is also the case with the heat of vaporization results, the main source of this discrepancy is the result for the shorter alkane chains. The better agreement with experiment of these parameters at modeling longer hydrocarbon chains may arise from the use of a high number of octane structures during torsion fitting. Furthermore, this parameter set is intended for the simulation of membrane lipids, which typically contain aliphatic tails that are ten carbon atoms or greater in length.

Table 3

Thermodynamic and Dynamic Properties of a Selection of Alkane Chains Simulated Using the Updated Lipid14 Parameters and Comparison to Experimenta

	ΔHvap (kJ mol-1)	ρ (kg m-3)	DPBC (10-5 cm2 s-1)	Dcorr (10-5 cm2 s-1)
Pentane
Lipid14	23.03 ± 0.16	592.45 ± 0.16	6.45 ± 0.56	7.1 ± 0.56
Expt	26.43[41]	626.2[41]		5.45[43]
Hexane
Lipid14	28.54 ± 0.1	636.3 ± 0.09	4.55 ± 0.29	5.02 ± 0.29
Expt	31.56[41]	656,[44] 660.6[41]		4.21[43]
Heptane
Lipid14	33.37 ± 0.11	667.31 ± 0.14	3.47 ± 0.23	3.85 ± 0.23
Expt	36.57[41]	679.5[41]		3.12[43]
Octane
Lipid14	38.67 ± 0.31	690.96 ± 0.10	2.11 ± 0.15	2.46 ± 0.15
Expt	41.49[41]	698.6[41]		2.354[45]
Decane
Lipid14	49.34 ± 0.30	724.47 ± 0.07	1.44 ± 0.15	1.65 ± 0.15
Expt	51.42[41]	726.6[41]		1.39[45]
Tridecane
Lipid14	64.62 ± 0.27	756.19 ± 0.24	0.48 ± 0.04	0.57 ± 0.04
Expt	66.68[41]	756.4[41]		0.712[45]
Pentadecane
Lipid14	74.99 ± 0.39	770.67 ± 0.25	0.30 ± 0.02	0.36 ± 0.02
Expt	76.77[41]	768.5[41]		0.461[45]

All values at 298.15 K.

Thermodynamic properties for a selection of pan class="Chemical">alkenes calculated using the updated LJ and torsion class="Chemical">parameters are shown in Table 4. As previously stated, properties of unsaturated chains are not as well reproduced as for class="Chemical">pan class="Chemical">alkanes due to the difficulty in tuning LJ parameters to achieve the experimental heat of vaporization, resulting in an RMS error of 13.80% for ΔHvap when compared to experiment. The density values are again better reproduced with an RMS error of 2.35%.

Table 4

Thermodynamic Properties of a Selection of Alkene Chains Simulated Using the Updated Lipid14 Parameters and, Where Available, Comparison to Experiment

	ΔHvap (kJ mol-1)	ρ (kg m-3)
cis-2-Hexene
Lipid14	26.17 ± 0.21	656.23 ± 0.13
Expt	32.19[41]	683[44]
cis-5-Decene
Lipid14	45.27 ± 0.22	739.19 ± 0.16
Expt	42.9[41]	744.5[41]
cis-7-Pentadecene
Lipid14	69.5 ± 0.22	781.44 ± 0.24
Expt		775[44]

The fractions of trans, gauche, n class="Disease">end gauche, double gauche, and kinked gtg′+gtg conformers per molecule were computed for the selection of class="Chemical">pan class="Chemical">alkanes under study (see Table 5). Experimental data, estimated by FTIR, exists only for tridecane;[46] however, the updated Lipid14 parameters reproduce these results extremely well with an overestimation of the end gauche and double gauche conformations only. Furthermore, the population of gauche conformations per molecule falls close to the experimental value of 35% for all chains investigated (t/g ratio ∼1.86), meaning that the overpopulation of trans conformations which drives GAFF bilayer simulations into the gel phase is avoided.

Table 5

Average Number of trans, gauche, End gauche (eg), Double gauche (gg), and gtg′+gtg Conformers Per Alkane Molecule and Comparison to Experiment

	trans	gauche	t/g ratio	eg	gg	gtg′+gtg
Pentane
Lipid14	1.20	0.80	1.49	0.80	0.13	-
Hexane
Lipid14	1.83	1.17	1.57	0.83	0.22	0.14
Heptane
Lipid14	2.49	1.51	1.65	0.83	0.31	0.24
Octane
Lipid14	3.15	1.85	1.71	0.81	0.39	0.33
Decane
Lipid14	4.47	2.53	1.77	0.81	0.57	0.54
Tridecane
Lipid14	6.47	3.53	1.84	0.81	0.82	0.83
Expt[46]	6.5	3.5	1.86	0.68	0.64	0.77
Pentadecane
Lipid14	7.80	4.20	1.86	0.81	1.00	1.02

The quality of the n class="Chemical">hydrocarbon class="Chemical">parameters was further assessed by calculating the class="Chemical">pan class="Chemical">13C NMR T1 relaxation times for heptane, decane, tridecane, and pentadecane. Similar to diffusion, this is a dynamic property; however, it depends on the rotation of the CH vector. Results are shown in Figure 6. In general the simulation values follow the same profile as the experimental data. Simulation values tend to converge with experiment upon moving further from the end carbon of the alkane chains. Although the result for pentadecane is slightly high, the overall comparison between simulation and experiment is reasonable.

Figure 6

Calculated 13C NMR T1 relaxation times for selected alkane chains and comparison to experiment.[47] Values at 312 K.

Calculated pan class="Chemical">13C NMR T1 relaxation times for selected class="Chemical">pan class="Chemical">alkane chains and comparison to experiment.[47] Values at 312 K.

Lipid Bilayer Simulation

Initial Structures

Bilayers were constructed using the CHARMM Membrane Builder GUI[48] at the relevant experimental hydration level (see Table 6) and converted to n class="Chemical">Lipid14 PDB format using the charmmclass="Chemical">pan class="Chemical">lipid2amber.x script.[17] All systems used the TIP3P water model[49] and had 0.15 M KCl salt concentration added to the water layer, modeled using suitable AMBER parameters.[50]

Table 6

System Size, Hydration, Temperature, and Simulation Time for the Lipid Bilayer Systems

	no. of lipids	simulation time (ns)	temp (K)	waters/lipid nW
DLPC	128	5 × 125	303	31.3
DMPC	128	5 × 125	303	25.6
DPPC	128	5 × 125	323	30.1
DOPC	128	5 × 125	303	32.8
POPC	128	5 × 125	303	31
POPE	128	5 × 125	310	32

Equilibration Procedure

The full system was minimized for 10000 steps, of which the first 5000 steps used the steepest descent method and the remaining steps used the conjugate gradient method.[51]

The system was then heated from 0 K to 100 K using Langevin dynamics[52] for 5 ps at constant volume, with weak restraints on the pan class="Chemical">lipid (force constant 10 kcal mol–1 Å–2).

Following this, the volume was allowed to change freely and the temperature increased to a pan class="Chemical">lipid dependent value (see Table 6) with a Langevin collision frequency of γ = 1.0 ps–1, and anisotropic Berendsen regulation[37] (1 atm) with a time constant of 2 ps for 100 ps. The same weak restraint of 10 kcal mol–1 Å–2 was maintained on the class="Chemical">pan class="Chemical">lipid molecules.

Production Runs

Constant pressure and constant temperature (NPT) runs were performed on the six bilayers using the AMBER 12 package.[17] The GPU implementation of the AMBER 12 code (bugfix 21) was used to run the simulations on NVIDIA GPU cards, achieving approximately 30 ns per day for the 128-n class="Chemical">lipid bilayer systems.[17,19] Three dimensional periodic boundary conditions with the usual minimum image convention were employed. Bonds involving class="Chemical">pan class="Chemical">hydrogen were constrained using the SHAKE algorithm,[36] allowing a 2 fs time step. Structural data was recorded every 10 ps. 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. All simulations were performed at constant pressure of 1 atm and constant target temperature (Table 6). Temperature was controlled by the Langevin thermostat,[52] with a collision frequency of γ = 1.0 ps–1, as this method was identified as the most suitable in previous work.[25] Pressure was regulated by the anisotropic Berendsen method[37] (1 atm) with a pressure relaxation time of 1.0 ps.

Each pan class="Chemical">lipid type was simulated for 125 ns with five repeats. The first 25 ns of each run was removed for equilibration, resulting in a total of 500 ns of data per class="Chemical">pan class="Chemical">lipid system, an aggregate of 3 μs of data. Results are presented as block averages over the five repeats ± standard deviation. The majority of analysis in this paper used PTRAJ or CPPTRAJ analysis routines.[17,53]

To check stability over time of the pan class="Chemical">lipid bilayer systems, the simulations were extended from 125 ns to 250 ns. Additional GPU and CPU validations were also performed, and in all cases the GPU results were consistent with CPU results (see the Supporting Information).

Validation

Bilayer Structural Properties

Despite the degree of uncertainty in obtaining accurate experimental values,[72] the bilayer surface area each n class="Chemical">lipid occupies, or area per class="Chemical">pan class="Chemical">lipid, is easily calculated from membrane simulations and gives a quick indication of whether a bilayer is in the correct phase at a given temperature. The area per lipid for each system was calculated using the dimensions of the simulation box as per previous work.[22,26] The AL for each lipid type is reported in Table 7, with all simulation values within 3% of experimental values, indicating that all the bilayers are in the correct Lα-phase. The result for POPE is closer to the older experimental AL value of 56.6 Å2 than the more recent AL value 59–60 Å2. Nevertheless, the AL should be but one of a number of properties calculated to validate a lipid force field.[73,74]

Table 7

Average Structural Properties over Five Repeats of the Six Lipid Systems Simulated with Lipid14 and Comparison to Experiment

lipid system	area per lipid AL (Å2)	volume per lipid VL (Å3)	isothermal area compressibility modulus KA (mNm-1)	bilayer thickness DHH (Å)	bilayer Luzzati thickness DB (Å)	ΔDB-H = (DB–DHH)/2 (Å)	ratio r of terminal methyl to methylene volume
DLPC
Lipid14	63.0 ± 0.2	948.9 ± 0.3	281 ± 37	30.4 ± 0.4	30.2 ± 0.1	–0.1 ± 0.2	1.9
Expt	63.2,[54]60.8[55]	991[54]	-	30.8[54]	31.4[54]	0.8[56]	1.8–2.1[57]
DMPC
Lipid14	59.7 ± 0.7	1050.2 ± 1.5	264 ± 90	34.7 ± 0.6	35.2 ± 0.4	0.3 ± 0.2	2.2
Expt	60.6,[54]59.9[55]	1101[4,54]	234[58]	34.4,[59] 35.3[60]	36.3,[54] 36.7,[55]	0.8[56]	1.8–2.1[57]
					36.9[59]
DPPC
Lipid14	62.0 ± 0.3	1177.3 ± 0.5	244 ± 50	37.9 ± 0.5	38.0 ± 0.2	0.1 ± 0.2	2.1
Expt	63.1,[55]64.3[61]	1232[4]	231[4]	38,[62] 38.3[4]	39.0[55,62]	0.8[56]	1.8–2.1[57]
DOPC
Lipid14	69.0 ± 0.3	1249.6 ± 0.2	338 ± 31	37.0 ± 0.2	36.2 ± 0.2	–0.4 ± 0.1	2.1
Expt	67.4,[62] 72.5[4]	1303[4]	265,[58] 300,[63] 318[64]	35.3,[65] 36.7,[62,66] 36.9,[4] 37.1[67]	35.9,[4] 36.1,[65,67] 38.7[62]	1.0–1.7[57]	1.8–2.1[57]
POPC
Lipid14	65.6 ± 0.5	1205.4 ± 0.4	257 ± 47	36.9 ± 0.6	36.8 ± 0.3	–0.1 ± 0.2	1.9
Expt	64.3,[55] 68.3[68]	1256[68]	180–330[69]	37[68]	36.8,[68] 39.1[55]	0.8[56]	1.8–2.1[57]
POPE
Lipid14	55.5 ± 0.2	1138.7 ± 0.3	350 ± 81	42.4 ± 0.2	41.0 ± 0.1	–0.7 ± 0.1	2.0
Expt	56.6,[70] 59–60[71]	1180[71]	233[70]	39.5[71]	-	-	1.8–2.1[57]

Experimentally, the volume per n class="Chemical">lipid VL is more accurately measured and is thus a better comclass="Chemical">parison for simulation results than AL. The volume per class="Chemical">pan class="Chemical">lipid was calculated using the dimensions of the simulation box[26] and the volume of a water molecule as determined by simulating 1936 TIP3P waters in the NPT ensemble for 50 ns using an identical procedure to the bilayer simulations at the relevant temperature.

VL values for each pan class="Chemical">lipid are reported in Table 7, which although systematically underestimated, are within 5% of experimental values. It is likely that the headgroup LJ class="Chemical">parameters could be further tuned to remedy this discreclass="Chemical">pancy, as the thorough reclass="Chemical">paramaterization of the class="Chemical">pan class="Chemical">hydrocarbon chains makes it unlikely that the tail groups are causing this lack of agreement.

This intuition is confirmed when studying the n class="Chemical">lipid component volumes calculated with the SIMtoEXP software.[75] The headgroup volume of class="Chemical">pan class="Chemical">DOPC was found to be 305.41 Å3, which is below the experimental estimate of 319–331 Å3, while the hydrocarbon chain volume of 965.88 Å3 is closer to the experimental range of 972–984 Å3.[57]

The volume breakdown provided by SIMtoEXP was used to calculate the ratio r of terminal methyl to methylene volume. All pan class="Chemical">lipid systems report a value of r = 1.85–2.17, within or very close to the experimental range of 1.8–2.1.

Isothermal area compressibility modulus, KA, was calculated from the fluctuation in the area per n class="Chemical">lipid.[26] In general, KA values fall close to experiment, with experimental values falling within the standard deviation of class="Chemical">pan class="Chemical">DMPC, DPPC, DOPC and POPC simulation results; however the POPE value comes out high and there is a large standard deviation in all values. Although the DOPC value is above the published experimental value of 300 mN m–1,[63] a personal communication with E. Evans revealed that this KA value has recently been revised upward to 318 mN m–1,[64] closer to the Lipid14 simulation result. This was not known prior to the lipid simulations.

In this work, the Berendsen method was used for pressure coupling, given that it is the only barostat currently available in the AMBER MD package. It has recently been shown that Berendsen pressure control is not ideal for simulations in which volume fluctuations are important,[76] thus by implementing other barostats into AMBER better KA results may potentially be achieved. This is a work in progress, the results of which will be shown elsewhere. Furthermore, larger system sizes and longer simulation times could also be investigated, as such changes have been shown to speed up the convergence of KA values.[74]

The membrane thickness was examined by calculating DHH, the peak-to-peak distance, from electron density profiles of the membranes. Again, satisfactory agreement with experiment is achieved for all pan class="Chemical">lipids, though the POPE value is a little high, indicating that this system is slightly too ordered.

An alternative bilayer thickness, the Luzzati thickness DB, was calculated using the z-dimension of the simulation box and the integral of the probability distribution of the n class="Chemical">water density along the z-axis.[9,10]DB values are found to lie close to experimental values, though the DB thicknesses for the saturated class="Chemical">pan class="Chemical">lipids are slightly underestimated. Given that DB is the distance between the points along the membrane normal at which the water density is half of its bulk value, this suggests that water is penetrating slightly too far into the hydrophobic region of the bilayer, thereby lowering the value of DB.

Ordering and Conformation of Lipid Acyl Chains

The ordering of the n class="Chemical">lipid acyl chains may be determined by calculation of the order class="Chemical">parameter class="Chemical">pan class="Disease">SCD. This quantity can be directly compared to experimental SCD values determined by 2H NMR or 1H–13C NMR.[77−80]SCD is a measure of the relative orientation of the C–D bonds with respect to the bilayer normal and was calculated according towhere θ is the angle between the bilayer normal and the vector joining C to its deuterium atom, and < > represents an ensemble average.

Simulation NMR order parameters for the six pan class="Chemical">lipid systems and comclass="Chemical">parison to experiment.[77,78,80−84]

Figure 7 shows the n class="Chemical">Lipid14 order class="Chemical">parameters with comclass="Chemical">parison to experiment. All class="Chemical">pan class="Chemical">lipid systems follow the experimental order parameter trend. The carbon-2 atom of the sn-1 and sn-2 chains display markedly different order parameters owing to the different alignment of the acyl chains in this region. Experimentally, it has been found that the SCD order parameter of the C-D bonds near the headgroup in the sn-1 chains are greater than the sn-2 chains.[85,86] Splitting of the order parameter value of the sn-2 chain from the sn-1 chain is observed for the simulated lipid systems. The unsaturated chain of the DOPC, POPC and POPE lipids show a distinctive drop at the carbon-9 and -10 positions due to the cis double bond. The S values for the sn-1 chain of POPE are a little high, indicating that POPE may be slightly too ordered. In agreement with Jämbeck et al., the two chains of DOPC show differing behavior, the sn-1 chain having higher S values about the double bond than the sn-2 chain.[87]

Figure 7

Simulation NMR order parameters for the six lipid systems and comparison to experiment.[77,78,80−84]

The conformation of the acyl chains may be examined by analyzing the rotamers and rotamer sequences along the n class="Chemical">lipid tails and comclass="Chemical">paring results to experimental data collected by Fourier transform infrared (FTIR) spectroscopy.[88] FTIR can determine the number of trans (t) and gauche (g) conformers and sequences of t and g (class="Chemical">pan class="Disease">end gauche eg, double gauche gg, gtg, and kinks gtg′). The lipid bilayer simulations were analyzed by denoting torsion angles φ in the acyl chains as either t (φ< −150° or φ> 150°), g- (−90°≤ φ< −30°) or g+ (30°< φ≤90°).[89] The rotamer sequence gtg correspond to g+tg+ or g-tg- while the sequence (or kinks) gtg′ corresponds to g+tg- or g-tg+.

Results are shown in Table 8 and are in general satisfactorily close to available experimental values. The discrepancies observed between simulation and experimental values of gtg′ for DLPC and n class="Chemical">DPPC may result from the experimental ambiguity in assigning gtg and gtg′ wagging modes.[90] These results also confirm that the bilayers are in the correct phase, as the gel-to-liquid phase transition is associated with an increase in the number of gauche rotamers and kink rotamer sequences.[90−92] However the eg and gg results for POPE are not in accordance with the experimental values obtained by Senak et al. using FTIR,[88] who found a marked increase in eg, gg, and gtg′+gtg values going from class="Chemical">pan class="Chemical">DPPE to DPPC because of the tighter packing of the PE lipid in the Lα phase. The present simulation values for POPE and POPC, though, are similar.

Table 8

Analysis of Rotamers and Rotamer Sequences in the Acyl Chains of the Six Lipid Systems – End gauche (eg), Double gauche (gg), Kinks (gtg′), gtg′+gtg, and Number of gauche (ng)

lipid system	eg	gg	gtg′	gtg′+gtg	ng
DLPC
Lipid14	0.35	0.44	0.28	0.52	2.50
Expt[93]	0.45	0.32	0.88*	-	2.85
DMPC
Lipid14	0.34	0.48	0.35	0.62	2.82
Expt[94]	0.38	0.67	-	0.44	2.6
DPPC
Lipid14	0.36	0.66	0.47	0.83	3.58
Expt	0.38,[94] 0.4,[88] 0.54[93]	0.4,[88,93] 0.57[94]	1.19[93]a	0.46,[94] 1.0[88]	2.44,[94] 3.6–4.2,[95] 3.7,[93] 3.8[82]
DOPC
Lipid14	0.36	0.75	0.37	0.70	3.93
POPC
Lipid14	0.36	0.69	0.41	0.75	3.73
POPE
Lipid14	0.35	0.60	0.42	0.73	3.50
Expt[88]	0.05	0.2	-	0.8	-

The gtg′ sequence may be ascribed to a gtg′+gtg sequence.[90]

Electron Density Profiles

The electron density profiles (EDP) were calculated by assuming an electron charge equal to the atomic number minus the atomic partial charge, located at the center of each atom. Profiles have also been decomposed into contributions from the following groups: n class="Chemical">water, class="Chemical">pan class="Chemical">choline (CHOL), phosphate (PO4), glycerol (GLY), carbonyl (COO), methylene (CH2), unsaturated CH=CH and terminal methyls (CH3). These profiles, shown in Figure 8, are all symmetrical, with water penetrating up to the carbonyl groups, leaving the terminal methyl groups dehydrated in agreement with experimental findings.[54,68,61] The electron density profiles were then utilized for the calculation of scattering form factors using the SIMtoEXP software.[75]

Figure 8

The total and decomposed electron density profiles for each of the six lipid bilayer systems with contributions from water, choline (CHOL), phosphate (PO4), glycerol (GLY), carbonyl (COO), methylene (CH2), unsaturated CH=CH and terminal methyls (CH3).

The total and decomposed electron density profiles for each of the six n class="Chemical">lipid bilayer systems with contributions from class="Chemical">pan class="Chemical">water, choline (CHOL), phosphate (PO4), glycerol (GLY), carbonyl (COO), methylene (CH2), unsaturated CH=CH and terminal methyls (CH3).

Scattering Form Factors

Scattering data allow direct comparison between pan class="Chemical">lipid bilayer simulation and experiment, avoiding any intermediate modeling of experimental raw data.[57] X-ray and neutron scattering form factors can be computed by Fourier transformation of simulation electron density profiles and comclass="Chemical">pared to experimental scattering data.

Recent work determined the areas per n class="Chemical">lipid (AX and AN) at which class="Chemical">pan class="Chemical">DOPC bilayer simulations best replicate the experimental X-ray scattering and neutron scattering data by varying the area per lipid through application of a surface tension, with the ideal situation being AX = AN = ANPT (bilayer is run in the tensionless NPT ensemble).[57] In this work we were concerned with validating the Lipid14 parameters for tensionless bilayer simulation only; thus we report the X-ray and neutron scattering form factors for ANPT only.

Simulation X-ray scattering form factors for the six n class="Chemical">lipid systems (black line) and comclass="Chemical">parison to experiment[54,55,62,66,68] (cyan circles). Inset: Simulation neutron scattering form factors at 100% class="Chemical">pan class="Chemical">D2O (black line), 70% D2O (red line), and 50% D2O (blue line) and comparison to experiment[55,96] (black, red, and blue circles, respectively).

It can be seen from Figure 9 that there is general agreement between both the X-ray and neutron scattering form factors for all pan class="Chemical">lipids for which there is experimental scattering data available, indicating that the simulated bilayers have the correct structure. In general the minima of the experimental F(q) profiles are correctly reproduced, as are the relative lobe heights.

Figure 9

Simulation X-ray scattering form factors for the six lipid systems (black line) and comparison to experiment[54,55,62,66,68] (cyan circles). Inset: Simulation neutron scattering form factors at 100% D2O (black line), 70% D2O (red line), and 50% D2O (blue line) and comparison to experiment[55,96] (black, red, and blue circles, respectively).

The quantity ΔDB-H was computed from the membrane thickness values (see Table 7). Agreement with X-ray scattering data is sensitive to the value of DHH, while agreement with neutron scattering data is sensitive to the value of DB. Therefore it has been proposed that bilayer simulations should aim to replicate experimental ΔDB-H values to best achieve agreement with both types of scattering data, where ΔDB-H = (DB-DHH)/2.[57] The GROMOS united-atom n class="Chemical">lipid force field has been shown to match experiment for simulation ΔDB-H results.[56,57] As evidenced by Table 7, class="Chemical">pan class="Chemical">Lipid14 ΔDB-H values are lower than those found by experiment, though all simulation values do maintain a large standard deviation. In fact, analysis of ΔDB-H results for two other all-atom lipid force fields, CHARMM36[8] and Slipids[10,87] indicates that this quantity is difficult to reproduce using all-atom models, with only the Slipids POPC result falling close to experiment. Figure 10 plots ΔDB-H values against area per lipid for CHARMM36,[56] Slipids,[10,87] and AMBER Lipid14, displaying a downward trend in ΔDB-H with increasing area per lipid, in disagreement with the experimental trend. Results for Lipid14 are similar to CHARMM36 and most Slipids values. Improving this discrepancy with experiment for Lipid14 may further improve the comparison with scattering data; however, present simulation scattering profiles are still seen to be in satisfactory agreement with experiment.

Figure 10

Plot of ΔDB-H versus area per lipid AL for the three all-atom lipid force fields CHARMM36 (squares), Slipids (diamonds), and AMBER Lipid14 (circles). Values shown for DLPC (green), DMPC (magenta), DPPC (blue), DOPC (red), and POPC (orange).

Plot of ΔDB-H versus area per n class="Chemical">lipid AL for the three all-atom class="Chemical">pan class="Chemical">lipid force fields CHARMM36 (squares), Slipids (diamonds), and AMBER Lipid14 (circles). Values shown for DLPC (green), DMPC (magenta), DPPC (blue), DOPC (red), and POPC (orange).

Lipid Lateral Diffusion

To assess the ability of the n class="Chemical">Lipid14 class="Chemical">parameters to reproduce dynamic class="Chemical">pan class="Chemical">lipid properties, the lipid lateral diffusion coefficient D was calculated using the Einstein relation (eq 2) with two degrees of freedom (nf = 2). Diffusion coefficients were computed for each lipid as a block average over the five NPT production runs. The mean-square-displacement (MSD) curves were determined using window lengths spanning 20 ns and averaged over different time origins separated by 200 ps. The slope of this curve yields the diffusion coefficient using eq 2, with the linear 10–20 ns region used to perform the fit. Prior to the MSD calculation, the lipid coordinates were corrected to remove the artificial center of mass drift of each monolayer.[73]

Results are of the same order of magnitude as experimental values; although in general they are underestimated. Unlike the bulk n class="Chemical">alkane work there is no correction term to account for collective motion which cannot be sampled using a periodic box of limited size. Accordinclass="Chemical">pan class="Chemical">gly, the underestimation may be a result of size effects. As highlighted by Poger et al. there is a widespread in experimental lipid lateral diffusion values in the literature, with a range of experimental techniques applied to the calculation of diffusion values.[105] Even different groups applying the same experimental technique do not necessarily yield comparable diffusion coefficients. Our calculated diffusion coefficients are nonetheless found to be in good agreement with other simulation values.[10,23,87,105−107]

Given that the production runs were performed in the NPT ensemble and that temperature regulation methods such as Langevin dynamics, which randomizes particle velocities, may affect dynamic properties such as diffusion, the n class="Chemical">lipid lateral diffusion was also determined in the microcanonical (NVE) ensemble. A single production run of each class="Chemical">pan class="Chemical">lipid system was extended into the NVE ensemble for 100 ns using the same simulation settings as used for the alkane diffusion runs (see Parameterization Strategy). Resulting time averaged MSD curves are shown in Figure 11, and calculated diffusion coefficients are reported in Table 9. Although similar to the D values determined from the NPT runs, the diffusion coefficients from the NVE runs are slightly higher, with the results for DPPC and DOPC showing the largest differences. This supports a recent study on the effect of temperature control on dynamic properties by Basconi et al.,[108] who found that diffusion coefficients calculated from simulations applying Langevin dynamics approach the coefficients derived from NVE simulations, provided weak coupling is used for the temperature regulation.

Figure 11

Time averaged mean square displacement of the center of mass of the lipid molecules versus NVE simulation time.

Table 9

Lipid Lateral Diffusion Coefficients Calculated from NPT Runs, NVE Runs, and Experimental Values

lipid system	calcd NPT Dxy (10-8 cm2 s-1)	calcd NVE Dxy (10-8 cm2 s-1)	simulation temp (K)	exptl Dxy (10-8 cm2 s-1)	exptl temp (K)
DLPC	7.65	7.78	303	8.5[97]	298
DMPC	5.05	6.32	303	5.95,[98] 9[99,100]	303, 303
DPPC	9.21	11.94	323	12.5,[101] 15.2[102]	323, 323
DOPC	6.48	9.49	303	11.5,[100] 17[103]	303, 308
POPC	5.74	6.54	303	10.7[100]	303
POPE	4.67	4.85	310	5.2[104]a	305

Cell culture membrane containing 78% POPE at 305 K.

Time averaged mean square displacement of the center of mass of the pan class="Chemical">lipid molecules versus NVE simulation time.

n class="Chemical">Lipid lateral diffusion is known to occur via two regimes: fast ‘rattling’ of the class="Chemical">pan class="Chemical">lipid in the local solvation cage[109] and slower, long distance diffusion in the plane of the bilayer. The two regimes are clearly observed in the MSD versus time curves (Figure 11) and are also revealed by computing the diffusion coefficients using different time ranges of the MSD curve. Figure 12 plots the diffusion coefficient D against the starting time for fitting the MSD slope. Time windows were either 100 ps long (fit starts between 10 ps and 500 ps) or 500 ps long (fit starts between 500 ps and 20 ns).[23]D values decrease smoothly and then converge at a value representing the long time diffusion of lipids in the bilayer plane.

Figure 12

Lateral diffusion coefficients for the six lipid types calculated using different time ranges of the mean square displacement curve for the linear fit.

Lateral diffusion coefficients for the six pan class="Chemical">lipid types calculated using different time ranges of the mean square displacement curve for the linear fit.

Conclusions

The n class="Chemical">Lipid14 force field represents a significant advancement for the simulation of phospohoclass="Chemical">pan class="Chemical">lipids in the AMBER MD package. Hydrocarbon parameters have been refined, resulting in good reproduction of thermodynamic and dynamic properties for a number of simple carbon chains, thus we can be confident that the hydrocarbon region of the lipid membrane is correctly represented. Head group parameters have also been updated, with the final parameter set finding good agreement with experiment for a range of properties, including the area per lipid, volume per lipid, bilayer thickness, NMR order parameters, scattering data, and lipid lateral diffusion, without applying a surface tension in the simulations. Crucially, the experimental raw data that requires no empirical input to derive, namely the NMR order parameters and scattering data, are well reproduced. Results for POPE however indicate that PE lipids may require further attention, as the order parameter results for POPE indicate that it remains somewhat artificially ordered in comparison to experiment. Results from five GPU repeats and CPU runs are seen to be consistent (these additional results are provided in the Supporting Information), with a number of tests performed on GPUs using both different starting structures and extending production runs to 250 ns. Future improvements may involve further refinement of parameters in order to address the underestimation of the volume per lipid and bilayer Luzzati thickness values in addition to PE lipid types.

Although the present study only concerns the validation of n class="Chemical">Lipid14 for the simulation of three saturated and three unsaturated class="Chemical">pan class="Chemical">lipids, the modular nature of the Lipid14 force field allows for a number of different lipids to be constructed from headgroup and tail group ‘building blocks’ and for the easy insertion of new lipid species into the force field. The Lipid14 charge derivation follows the usual AMBER convention, making this force field compatible with other AMBER potentials, such as the General Amber Force Field[16] and the ff99SB protein force field.[14] As such, the interaction of other species, such as small molecules or proteins, with lipid membranes can be studied in AMBER using the Lipid14 force field.