Derivation and systematic validation of a refined all-atom force field for phosphatidylcholine lipids.

An all-atomistic force field (FF) has been developed for fully saturated phospholipids. The parametrization has been largely based on high-level ab initio calculations in order to keep the empirical input to a minimum. Parameters for the lipid chains have been developed based on knowledge about bulk alkane liquids, for which thermodynamic and dynamic data are excellently reproduced. The FFs ability to simulate lipid bilayers in the liquid crystalline phase in a tensionless ensemble was tested in simulations of three lipids: 1,2-diauroyl-sn-glycero-3-phospocholine (DLPC), 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), and 1,2-dipalmitoyl-sn-glycero-3-phospcholine (DPPC). Computed areas and volumes per lipid, and three different kinds of bilayer thicknesses, have been investigated. Most importantly NMR order parameters and scattering form factors agree in an excellent manner with experimental data under a range of temperatures. Further, the compatibility with the AMBER FF for biomolecules as well as the ability to simulate bilayers in gel phase was demonstrated. Overall, the FF presented here provides the important balance between the hydrophilic and hydrophobic forces present in lipid bilayers and therefore can be used for more complicated studies of realistic biological membranes with protein insertions.

Introduction

Biological membranes are crucial components of the cell. They divide the inner part of the cell from the outer environment as well as dictate the intercellular transport and cell fusion/division. Membranes also provide the matrix in which important proteins (transmembrane pumps, channels, receptors etc.) reside. In nature, biological membranes are very complex systems, consisting of a mixture of different lipids, sterols, and proteins. These constituting parts are involved in a number of processes making the cell membrane a very dynamic structure. It has been proposed that protein–lipid interactions regulate the functionality of proteins, making the lipid bilayer more than a passive hydrophobic slab.[1−7] Most studies are, however, performed on pure phospholipid bilayers since phospholipids are a major part of biological membranes. Before one can study biological membranes in their true essence, the properties of simpler membranes must be fully understood. Experimental techniques used to study lipid bilayers span from X-ray and neutron scattering to IR/Raman and NMR.[8−10]

Lipid bilayers have a very pronounced temperature-dependent phase behavior. The biologically most relevant phase for membranes is the liquid-crystalline phase (Lα-phase), which means that the aliphatic chains of the lipid have a certain degree of freedom which result in a disordered structure. Although a large number of experimental studies have been performed on membrane systems and already shed some light on the dynamical and structural properties, a more detailed understanding of these systems is desirable. Due to the difficulties of studying these systems in the correct phase, computer simulation methods are an ideal complement to experiments. Molecular dynamics (MD) computer simulations have become an essential tool for studying membrane on an atomistic level during the past two decades.[11−17]

If a sufficiently large system is chosen to be simulated, the quality of the simulations is limited by two factors: (1) sampling or time scale and (2) the empirical potential energy function or force field (FF) employed. The sampling issue is always present and it has been found that fluctuations in simulations of pure lipid bilayers happen on a multinanosecond time scale.[18] As a consequence of this, long simulations (>100 ns) are needed in order to study equilibrium properties of lipid bilayers. For systems with proteins embedded, much longer simulations may be needed.[19,20] The other issue, i.e., the ability of the potential energy function to provide realistic dynamics and thermodynamic properties of the system, is extremely crucial. Two major types of FFs are available for lipid bilayers: united atom (UA) and all-atomistic (AA). In UA FFs all nonpolar hydrogens have been included in the heavier atoms, creating pseudoatoms. The popular Berger lipids,[13] the G53A6L FF,[21] the 43A1-S3 FF,[22] and the potential developed by Ulmschneider and Ulmschneider[23] and Kukol[24] are examples of these FFs. The AA FFs have explicit hydrogens included and the most common ones for lipid simulations are CHARMM[25−27] and the general AMBER FF (GAFF).[28] Until now, these have been the most widely used AA FFs, although several problems have been reported regarding older versions of CHARMM[29−31] and GAFF FFs.[32−34] Other FFs are nonadditive (polarizable)[35−37] or coarse-grained.[38−40] In some cases, an AA description is preferred and recent studies have shown that there exists a unique designated hydrogen bond network around the choline head group.[41,42] In order to describe and study this in full detail, an accurate AA model is crucial. AA description may also be crucial for describing order parameters of the glycerol group and in the description of double bonds in unsaturated lipids. An inaccurate FF often results in a lipid bilayer in a too ordered phase, a gel phase (Lβ′), at a higher temperature than observed in experiments and this leads to a contracted surface area which results in a too low area per lipid (AL).[30,31,43] Most often, this has been a problem for the models including explicit hydrogens. In order to solve this problem, a positive surface tension can be applied, so the simulations are performed in a NγPT ensemble. Another solution is to fix the total simulation box area leading to a NPAT ensemble. Neither of these two solutions is appealing since they restrict the membrane from stretching laterally and suppress membrane undulations. Therefore, the appropriate ensemble for studying lipid bilayers as well as insertions to the bilayer is the isothermal–isobaric ensemble (NPT) which has been explained in more depth earlier.[13,44,45]

It has proven to be very difficult to obtain accurate FFs for lipid bilayers and often specialized FFs have to be used; e.g., CHARMM uses specialized parameters for lipids that differ from the parameters used for proteins to describe the same type of atoms. One possible explanation for this is the anisotropic nature of these systems. Lipid bilayers are very peculiar systems with a polar head group in contact with the aqueous environment and a strong attraction of the tails driven by the hydrophobic effect and van der Waals (vdW) forces. A complexity like this poses a challenge for the modeling community. Especially, the fact that in order to study complexes including both lipids and proteins we need FFs that work well together implies further complications. Previous work has shown that the interactions between amino acids and lipid molecules need to be carefully evaluated before embarking on larger simulations.[46]

In the present paper, we present an AA FF for fully saturated phosphocholine lipids based primary on high-level ab initio calculations together with some limited empirical input. Partial atomic charges for the whole lipid, Lennard-Jones (LJ) parameters, and torsion potentials for the aliphatic chains are recalculated. The introduction of the scaling factor of 0.8333 for 1–4 Coulombic interactions and 0.5 for 1–4 vdW interactions makes the parameters set presented here compatible with the AMBER FFs for proteins and biomolecules. Computing of partial atomic charges from the electrostatic potential (which is a quantum mechanical observable wheres the charges are not) makes these very essential FF parameters physically grounded. The fact that all charges are computed at the same level of theory and with the same procedure together with the LJ parameters and torsional potential results in a FF with balanced terms that have been developed in a consistent manner. This is important since it makes it easy to improve the FF in the future by using even more precise ab initio methods and this also creates a good ground for transferability of the computed parameters. With a FF well balanced and capable of describing the physics and chemistry of a lipid bilayer, computer simulations can be extended to do more than just reproduce experimental findings. Different properties and phenomena can be predicted by accurate models and the fundamentals of molecular biology can be better understood with the help of modeling. Naturally, verification against experimental findings plays a key role in development of these models, so a large number of properties for the lipid bilayers have been computed and compared with experimental data. Area per lipid is one important structural property of lipid bilayers. However, due to experimental inconsistencies it has been proposed that this should not be the only property investigated.[18,47,48] Therefore, different thicknesses, scattering form factors, volumes, and order parameters have been calculated and compared to experimental data. Due to the distinct thermal dependency of the aliphatic chains of the lipids, simulations at different temperatures were performed in order to verify that the set of parameters proposed here captures this effect. The final comparisons show that the FFs presented here reproduce the experimental findings very well and therefore is highly suited for simulations of lipid bilayers in the correct phase.

Methods and Models

Parameterization Strategy

The main focus of the parametrization procedure has been on the hydrophobic tails of the lipids, since these parameters were previously shown to have strong influence on simulated properties of bilayers, and the partial atomic charges for the whole lipid. As model compounds for the hydrophobic tails, a range of n-alkanes has been chosen (hexane to hexadecane). The assumption that the hydrophobic core behaves similar to a bulk alkane liquid has proven to be valid a number of times.[13,22,23,31,35,49] If the fitting is done simultaneously for a series of alkanes, the results are generally better.[49] Furthermore, studies on NMR relaxation rates have shown that the hydrophobic part of the lipid bilayer behaves similar to bulk alkane solutions.[50]

The following potential energy function was chosen for this workwhich is a standard form for FFs like AMBER, CHARMM, and GROMOS. As a starting point for our parametrization, we took parameters of the CHARMM36 (C36) FF.[27] Parameters for all covalent bonds and angles as well as LJ and torsional parameters for the lipid head group were taken from the FF described by Klauda et al.[27] We then derived new parameters (as described in detail below): partial atomic charges for the whole lipid and LJ and torsion parameters describing the lipid tails.

We started our parametrization by computing partial atomic charges for typical configurations of alkanes and averaging over these configurations. Because of the known difficulties of reproducing the vdW dispersion interaction by ab initio methods, experimental heats of vaporization and densities were used during the fitting of the LJ parameters. First, the new charges were used with the original parameters from C36 (LJ and torsional) and the LJ parameters were then altered until satisfactory agreement between simulations and experiments was obtained. After this, the torsional potentials were fitted from ab initio computations for the model compound. After one round in the parametrization scheme, it was necessary to refit the LJ parameters again and the torsional potentials until self-consistency was obtained. The Lorentz–Berthelot combination rules were used for the vdW interactions.[51] The introduction of scaling factors for the 1–4 interactions with the CHARMM FFs for lipids together with a new set of charges have earlier been proven to be successful.[31] Below, we present more detailed descriptions of each parametrization step.

Partial Atomic Charges

For the lipid tails hexadecane was chosen as a model compound for computing the charges. It is well-known that partial atomic charges are conformation dependent[52] but previous FFs have been parametrized from optimized geometries. In order to address this issue, we performed a 10 ns long MD simulation with pure hexadecane with FF parameters earlier derived by our group.[31] After the simulation, 54 random conformations were extracted and used for computing the charges which were then averaged over all conformations in order to obtain a final set. In this way, we obtained Boltzmann-averaged charges over an ensemble of conformations in a procedure equivalent to the one used by Sonne et al.[30] We hope that by averaging over an ensemble of conformations the effects of the conformational dependence of partial charges are minimized. In computation of atomic charges, a dielectric constant of 2.04 was used to mimic the dielectric environment of the membrane’s hydrophobic part.

Atomic charges for the lipid head group were obtained in a similar fashion where 26 random conformations were chosen from a 20 ns long simulation of an equilibrated bilayer (DMPC) with the same FF parameters used in the initial simulation of hexadecane. A large part of the hydrophobic parts of the lipids were then cut off in order to save CPU time and the cropped lipids were placed in dielectric continuum with ε = 78.4 in order to mimic the aqueous environment. Inclusion of solvent effects results in a FF with implicitly polarized charges optimized for condensed phase simulations. This has been proven to give reliable results without any performance loss.[53]

For each molecular conformation, the charges were computed using the restricted electrostatic potential approach[54] (RESP) with the DFT method using the B3LYP exchange-correlation functional[55−58] and the cc-pVTZ basis set.[59] The electrostatic potential was sampled with the Merz–Singh–Kollman scheme[60] by single-point calculations and fitted during the two-stage procedure developed by Cornell et al.[61] All solvent effects were modeled by placing the molecule in a polarizable continuum with different dielectric constant (see above) with the IEFPCM continuum solvent model.[62,63] The quantum mechanical calculations were performed with the Gaussian09 software package,[64] and the RESP calculations were performed with the Red software.[65] In subsequent molecular dynamics simulations, Coulombic 1–4 interactions were scaled by a factor of 0.8333.

The way the atomic charges have been calculated and used in MD simulations makes them compatible with the AMBER03 FF[53] and since the charges in all AMBER FFs are derived from the RESP the lipid FF presented here is compatible with most members of the AMBER FF family. This is of importance since there is a growing interest in simulating membrane proteins in their native environment[19,66] and also peptide partitioning in biological membranes.[67,68] Ongoing work aims to clarify which AMBER biomolecular FFs that work sufficiently well together with the current parameters. A preliminary test of the compatibility of the lipid parameters and the AMBER03 FF is presented further down.

Boltzmann averaging over charges introduces temperature dependency on the charges and in order to see the impact of temperature, simulations with different temperatures (298, 303, 310, 318, and 325 K) were performed with hexadecane using the methodology described above. No explicit temperature dependence could be found over this range of temperatures, making the charges reliable and robust with respect to temperature, at least within the interval tested here (data not shown).

Lennard-Jones Parameters

By reproducing experimental densities correct σ-values for the LJ potential can be assigned. Heats of vaporization are important since they describe the interactions between the molecules which is primarily determined by the well-depth (ε) of the LJ potential. The average density is straightforward to obtain from MD simulations and the heat of vaporization was computed according towhere ⟨E⟩gas was obtained as an average of eight separate 50 ns long, in vacuo simulations with different starting conformations taken from a condensed-phase simulation. A time step of 0.25 fs was used for the in vacuo simulations. Test simulations in condensed phase lasted for 5 ns and production simulations for 100 ns at 298.15 K with the partial atomic charges obtained in the previous step. The simulation length of the condensed phase is longer than usually considered but it is necessary since a variation of parameters during the fitting procedure can freeze the system and this transition can occur on a relatively long time scale.[23] The LJ 1-4 interactions between the hydrogens of the lipid carbon chains were scaled by a factor of 0.5 and for the remaining interactions special 1-4 interactions were specified. In order to make sure that the alkanes intermolecular interactions are of the right magnitude and that the molecules have the correct size, the LJ parameters ε and σ were fitted simultaneously.

Torsion Potential

The potential energy surface (PES) for octane was scanned in order to obtain an accurate torsion potential for the hydrophobic tails of the lipids. The torsion angles of interest were CH2–CH2–CH2–CH2 and CH2–CH2–CH2–CH3. The relative energies between different conformations were obtained by employing the hybrid method for interaction energies (HM-IE).[69] Applying this method makes it possible to estimate intermolecular and intramolecular interactions with the precision of very high level ab initio methods, namely CCSD(T)[70] with a large basis set (LBS). Assuming that the effect of using a larger basis set is additive, the following expression can be used to estimate the intermolecular energieswhere SBS is the abbreviation for the small basis set. This energy was estimated at a CCSD(T)/cc-pVQZ level of theory; i.e., LBS and SBS are equivalent to cc-pVQZ and cc-pVDZ, respectively. The geometry of the molecules was optimized with the second order of Møller–Plesset perturbation theory (MP2)[71] with the cc-pVDZ basis set and then single-point calculations were used to estimate the relative energies with the HM-IE scheme. It has been proven that MP2/cc-pVDZ gives reasonable geometries whereas it is necessary to use more exact methods for the relative energies.[72]

Since longer carbon chains give lower rotational barriers[73,74] and due to the fact that the number of medium-range 1,3-interactions increases the stability of the compounds,[75] a slightly larger model compound (octane) was used than that used in parametrization of CHARMM.[26] Different conformations were used during the scanning of the PES and used in the fitting procedure, resulting in a total number of roughly 200 conformations estimated at a CCSD(T)/cc-pVQZ level of theory. More detailed information regarding the fitting procedure can be found in the Supporting Information along with all FF parameters derived in this work.

Testing of Newly Derived Parameters

A number or properties were computed for the bulk alkane systems and several lipid bilayers (DLPC, DMPC, and DPPC) in order to test the quality of the newly derived parameters. For bulk alkane liquids, uncorrected (DPBC) and corrected (D) self-diffusion coefficients[76] were determined as well as 13C longitudinal relaxation times (T1) and the number of gauche bonds per molecule. The ability of the newly derived FF to reproduce experimental data for lipid bilayers was verified by computing a series of structural and dynamic properties: AL, volume per lipid (VL), isothermal area compressibility modulus (KA), thermal area expansivity (αAT), and lateral diffusion coefficient. Three types of bilayer thicknesses were computed and compared to experiments: head-to-head distance (DHH), the Luzzati thickness (DB),[8] the hydrophobic core thickness (2DC), and the two corresponding thermal contractivities, αDBT and αDCT. Electron density profiles along the bilayer normal were obtained from simulations with the contributions from different groups (see Figure 1). From the total electron density, the scattering form factors, F(q), were obtained via Fourier transformation. Deuterium order parameters (|SCD|) were also computed and compared to experimental findings. More information regarding the computed properties for both alkane liquids and lipid bilayers can be found in the Supporting Information.

Figure 1

Atomic charges calculated via the RESP method for the whole lipid. Charges within parentheses are the charges on the hydrogen(s) attached to the heavier atom. The second column corresponds to the groups used in Figure 7. α5, α4, and α3 are torsional angles that are discussed in the Results section.

Simulations Details

All simulations were performed in the isothermal–isobaric ensemble, NPT, at a pressure of 1 atm. The pressure was held constant by using the Parrinello–Rahman barostat[77] with a coupling constant of 10.0 ps with an isothermal compressibility of 4.5 × 10–5 bar–1. For the bulk liquids an isotropic pressure coupling was used and for the bilayer simulations a semi-isotropic pressure coupling scheme was used. The temperature was kept constant by the Nosé–Hoover thermostat[78,79] with a coupling constant of 0.5 ps. The lipid bilayer and water were coupled separately to the thermostat. Long-range electrostatic interactions were treated by a particle-mesh Ewald scheme[80,81] with a real-space cutoff at 1.4 nm with a Fourier spacing of 0.10 nm and a fourth-order interpolation to the Ewald mesh. Single-atom charge groups were used. van der Waals interactions were truncated at 1.5 nm and treated with a switch function from 1.4 nm. Long-range corrections for the potential and pressure were added.[51] The inclusion of long-range corrections should eliminate the LJ cutoff dependency in the simulations. Due to the fact that lipid bilayers are inhomogeneous systems the method introduced by Lagüe et al.[82] to add long-range corrections could be applied instead. Periodic boundary conditions were imposed in every dimension. A time step of 2 fs was used with a Leap-Frog integrator. The LINCS algorithm[83] was used to freeze all covalent bonds in the lipid, and the analytical SETTLE[84] method was used to hold the bonds and angle in water constant. The TIP3P water model[85] was the water model of choice. The choice of water model can be explained by the fact that TIP3P is the default water model in major FFs such as AMBER and CHARMM and since one of the aims of the work presented here was to create a lipid FF compatible with AMBER this was a natural choice. Further, earlier work of Högberg et al.[31] has shown that there is flexibility in the choice of water model for AA simulations of lipid bilayers. Atomic coordinates were saved every 1 ps and the neighbor list was updated every 10th step.

Bulk liquids were simulated with a simulation box consisting of 128 molecules for the larger alkanes and 256 for the smaller alkanes (hexane and heptane) at a temperature of 298.15 K. The lipid bilayer systems were prepared using the CHARMM-GUI[86,87] with 128 lipids in total, 64 in each leaflet. In order to achieve proper hydration, 30 TIP3P water molecules were added per lipid. Three different lipid types were simulated, DLPC (12:0/12:0), DMPC (14:0/14:0), and DPPC (16:0/16:0). These system were investigated under a range of temperatures; see Table 1 for an overview of all simulations performed. All lipid bilayer systems were equilibrated for 40 ns before production runs were initiated which lasted for 300–500 ns. All MD simulations were performed with the Gromacs[88] software package (versions 4.5.3 and 4.5.4). All analysis were made with the analysis tools that come with the MDynaMix software package.[89] System snapshots were rendered and analyzed with VMD.[90] Neutron scattering form factors were computed with the SIMtoEXP software.[91]

Table 1

List of Lipid Bilayer Simulations Performed in the Present Studya

lipid	temp (K)	time (ns)
DLPC	303	400
	323	400
	333	400
DMPC	303	400
	323	400
	338	400
DPPC	293	300
	323	500
	333	500
	353	500

All simulations were performed with 128 lipids and 30 water molecules per lipid. Total simulation time: 4.2 μs.

The calculations of free energies of solvation in water and cyclohexane were performed by using thermodynamic integration over 35 λ values in the range between 0 and 1. A soft core potential (SCP) was used to avoid singularities when the solute is almost decoupled from the solvent. The α-parameters used for the SCP and the simulation workflow were set following the methodology described by Sapay and Tieleman.[92] The amino acid analogues were solvated with 512 and 1536 molecules of cyclohexane and water, respectively.

Results and Discussion

Bulk Alkane Liquids

Our initial charge calculations implied that the major part of the alkanes have almost zero (less than 0.005) charge of carbons and hydrogens with small but noticeable dipole moments at the end of the chains. Therefore, charges were set to zero for the inner methylene groups (both for carbons and hydrogens). As a result of this, the performance gap between AA and UA FFs can be reduced significantly. The total charge on the outermost methylene groups was set to +0.033 which is balanced by the charge of −0.033 on the methyl groups. The partial atomic charges used during the simulations of bulk alkane solutions were the same as presented for the lipid tail in Figure 1. The major parts of the charges in the aliphatic tails are set to zero which is similar to most UA FFs[13,22,23,47] but differs from all version of CHARMM FFs[25−27] having a charge of +0.09 on all aliphatic hydrogens. The present set of charges differ from a majority of other FFs due to the presence of the small dipole moments at the end of the lipid tails.

Figure 7

Electron density profiles from simulations of DLPC (top), DMPC (middle), and DPPC (bottom) lipid bilayers. Contributions from individual components are shown with colors: N, choline groups, dark blue; P, phosphate groups, cyan; CO, ester groups, red; CH2:, methylene groups, light gray; CH3, methyl groups, dark gray.

The final LJ parameters for methylene and methyl groups are shown in Table 2 and the final torsional parameters in Table 3. The same number of terms in the Fourier sum with the same multiplicity and phase shift as used in C36 gave the best result.

Table 2

Final Values for the Optimized Lennard-Jones Parameters in This Work

atom type	σ (nm)	ε (kJ mol–1)
C (CH2)	0.358	0.228
H (CH2)	0.240	0.112
C (CH3)	0.350	0.340
H (CH3)	0.220	0.095

Table 3

Values for the Optimized Torsional Potentials

torsion	force const, kϕ(kJ mol–1)	mult., n	phase shift, δ (deg)
CH2–CH2–CH2–CH2	0.0603	2	0
	0.2036	3	180
	0.5541	4	0
	0.4020	5	0
CH2–CH2–CH2–CH3	–0.0835	2	0
	0.2659	3	180
	–0.3460	4	0
	0.4006	5	0

The results of computations with the final optimized parameter set for heat of vaporization, density, and diffusion are summarized in Table 4. The linear trend in ΔHvap is reproduced perfectly by the simulations with only a slight overestimation for the longest alkanes resulting in a rms error of only 0.9%. Further, the rms error for the density is 0.8% which means that the correct balance between size and interactions in the FF is obtained. During the parametrization procedure the alkane system can easily freeze due to incorrect parameters in a multinanosecond simulation.[23] A too high value of ε combined with too small atomic sizes (σ) gave rise to a stacked lamellae phase which was also observed in the present investigation. This behavior will affect the FFs ability to describe a lipid bilayer in the Lα-phase at the correct temperature in a tensionless ensemble. The freezing sometimes occurred after only a couple of nanoseconds, making it straightforward to avoid this problem. In some cases, longer simulations were needed just as it was for Ulmschneider and Ulmschneider.[23] Due to the possibilities of these phase transitions, all simulations were run for 100 ns.

Table 4

Heat of Vaporization (ΔHvap), Densities (ρ), and Uncorrected (DPBC) and Corrected Self-Diffusion Coefficients (D) from Simulations after the Final Parameterization Compared to Experimental Values

	ΔHvap(kJ mol–3)		ρ (kg m–3)		D (10–5cm2 s–1)
alkane	sim	expt[93]	sim	expt[93]	DPBC	D	expt[94,95]
hexane	31.86 ± 0.17	31.56	668.9 ± 0.22	660.6	3.18	3.72 ± 0.31	4.21
heptane	36.45 ± 0.14	36.57	687.1 ± 0.13	679.5	2.76	3.16 ± 0.13	3.12
octane	41.11 ± 0.11	41.49	702.3 ± 0.32	698.6	1.81	2.18 ± 0.08	2.36
nonane	46.98 ± 0.13	46.55	724.7 ± 0.15	719.2	1.51	1.78 ± 0.10	1.78
decane	51.04 ± 0.19	51.42	732.2 ± 0.11	726.6	1.10	1.31 ± 0.04	1.39
undecane	57.12 ± 0.10	56.58	746.6 ± 0.23	740.2	0.86	1.02 ± 0.06	1.11
dodecane	60.93 ± 0.15	61.52	754.4 ± 0.10	749.5	0.65	0.77 ± 0.06	0.87
tridecane	67.13 ± 0.24	66.68	761.0 ± 0.18	756.4	0.49	0.59 ± 0.05	0.71
tetradecane	71.10 ± 0.18	71.73	767.1 ± 0.21	759.6	0.51	0.59 ± 0.04	0.55
pentadecane	77.56 ± 0.26	76.77	772.0 ± 0.14	768.5	0.41	0.47 ± 0.08	0.46
hexadecane	82.02 ± 0.21	81.35	776.6 ± 0.10	770.1	0.31	0.36 ± 0.03	0.39
rms error	0.9%		0.8%		20%	8.6%

As the thermodynamic properties for the alkanes converged and resulted in good agreement with experimental data, the FFs ability to reproduce dynamic properties such as diffusion and NMR relaxation times was investigated. Since the FF only has been parametrized against thermodynamic data, it is of utmost interest to test the dynamics of the alkanes. The uncorrected diffusion coefficients, DPBC, in Table 4 are all underestimated due to the finite box size in the simulations. Once corrections for the hydrodynamics and long-range effects[76] were added, the rms error decreases from 20% to 8.6%. Similar trends have been reported by Klauda et al.[26] and Davis et al.[35] Figure 2 shows the 13C NMR T1 relaxation times obtained from simulations compared to experiments[96] (in the Supporting Information more details on how to obtain 13C NMR T1 from simulations are discussed). The dynamics of the chain is important since it will affect how the inner part of the bilayer behaves when solutes and proteins are inserted. For the shortest alkane chain (heptane) the relaxation time is overestimated which indicates that the local chain dynamics is too rapid. The agreement between simulations and experiments becomes better when the chain size is increased. The trend seen here follows the general trends found in the recent literature.[26,35,97] Albeit the fact that the longer chains have relaxation times closer to the experimental ones, the dynamics around the torsional angles are still slightly too slow. As can be seen further on, this does not affect the lipid bilayer properties but purposes a future improvement of the FF. Schemes where high-level ab initio data are used together with available experimental data in order to fine-tune parameters have been suggested before.[98] This method could be used to eventually reach the limit of accuracy of fixed point charge FFs.

Figure 2

Computed 13C NMR T1 relaxation times for alkanes for selected carbons compared to experimental data.[96] Note the difference in scale.

The statistics of different conformations of alkane chains were also investigated by computing the fraction of gauche conformations in condensed phase. In Table 5 the average number of gauche and trans bonds per alkane molecule are presented. The experimental estimation of 35% gauche conformations per molecule[99] is reproduced; however, this comparison should be done with some care since the interpretation of the experimental data is not trivial. Still, it is reassuring that a high enough number of gauche conformations occur since in some previous FFs the trans conformations have been too pronounced, resulting in lipid bilayers that are too contracted when simulated under zero surface tension.[31]

Table 5

Average Number of Gauche and Trans Bonds Per Alkane Molecule from Simulations and Experiment

alkane	gauche	trans	t/g ratio
hexane	1.07	1.93	1.80
heptane	1.58	2.42	1.53
octane	1.98	3.02	1.52
nonane	2.22	3.78	1.70
decane	2.54	4.46	1.76
undecane	2.85	5.15	1.81
dodecane	3.16	5.84	1.85
tridecane	3.47	6.53	1.88
tetradecane	3.78	7.22	1.91
pentadecane	4.10	7.90	1.92
hexadecane	4.39	8.61	1.95
experimental[99],a	3.50	6.50	1.86

For tridecane, estimated by FTIR.

Transfer free energies from a polar environment to an apolar environment have been used earlier to investigate the compatibility of lipid–protein FF combinations.[92,100] It has been shown that this scale, which is referred to as the Radzicka scale,[101] correlates with studies where the same amino acid analogues have been transferred from aqueous solution to an explicit membrane environment.[102−104] Similar approaches have also been used in the parametrization of the GROMOS FF.[105] The computed free energy of transfer, ΔΔG, together with free energy of solvation in water and cyclohexane are presented in Figure 3. The free energies of hydration indicate that the AMBER03 FF (used to describe amino acids) is slightly too hydrophobic with the MD simulations performed here when compared to experimental data.[106,107] This is a general feature of most major FFs developed for proteins.[92] The rms in ΔGhyd is 4.86 kJ mol–1 which does not differ much from the rms values obtained by Shirts et al.[108] and Sapay and Tieleman.[92] The largest errors are for HIE (7.77 kJ mol–1), ASN (6.61 kJ mol–1), and MET (6.93 kJ mol–1) which also correlates with earlier work.[92,108,109] Computed free energies of solvation of the amino acid analogues in cyclohexane are in better agreement with experimental data[101] and have a rms error of 2.76 kJ mol–1. Largest errors in ΔGcyc are obtained for TRP (5.68 kJ mol–1), CYS (4.86 kJ mol–1), and HIE (4.34 kJ mol–1). The transfer free energies also are in good agreement with experimental findings with a rms error of 3.40 kJ mol–1. Some error cancellation also occurs, e.g. for MET, HIE and TRP. The residues that contributes mostly to the error are GLN (6.75 kJ mol–1), ASN (6.26 kJ mol–1), and PHE (5.14 kJ mol–1). Overall, the free energies agree well with experiments and the partitioning of amino acid analogues between the polar and apolar solvents suggests that the LJ parameters presented here are compatible with the AMBER03 FF. Since there are no intermolecular interactions of Coulombic nature present in the cyclohexane simulations, it is likely that other AMBER FFs will work together with the FF presented here. The data presented here shows that despite the fact that a different cutoff scheme was applied than the one used in the derivation of AMBER03, the two parameter sets are very compatible. For future work, we will test explore this compatibility further by performing simulations of membrane-bound proteins in their native environment. However, one should bear in mind that this hydrophobicity scale is missing some crucial points as protein backbone contributions etc. and therefore is only a crude measurement.

Figure 3

Free energy of solvation for aminoacid analogues in water, ΔGhyd (left), cyclohexane, ΔGcyc (center), and free energy of transfer from water to cyclohexane, ΔΔG (right). The light gray area corresponds to an offset of ±kBT and the dark gray to an offset of ±2kBT. Errors in the computed values are within the dot size in the plots.

Lipid Bilayer Simulations

Atomic Charges for the Lipid Head Group

Partial atomic charges for the head group of the lipid were Boltzmann-averaged over a number of conformations generated from a 20 ns MD simulation and are presented in Figure 1. The sum of the individual atoms in the groups agree very well with the work of Högberg et al.[31] and Sonne et al.[30] which both describe AA FFs. The charge distribution in the head group, due to the long-range character of electrostatic interactions, has a major impact on the fluidity of the lipid bilayer.[30,31,43] Compared to most UA FFs[13,21−23,92] and the all-atomistic C36[27] FF the charges presented here differ quite a lot. For the former, certain charge groups are usually constrained to a certain charge, as +1 for the choline group and −1 for the phosphate group. In the present paper, no constraints were imposed, resulting in a slightly different charge distribution in the head group, +0.75 for the choline group and −0.86 for the phosphate group. This will affect the P–N dipole moment of the head group, making it a weaker dipole in the case of our charges and a stronger dipole when the charge groups have predetermined values. A decrease in dipole moment reduces the attraction of the lipid head group with its neighbors in the membrane plane resulting in a more upright position.[43,110] In Figure 4, the probability distribution of the angle between the bilayer normal and the P–N vector is shown for DLPC, DMPC, and DPPC. The most probable angle is around 71° for all lipids compared to the experimental estimate of 72°.[111] This means that the dipole lies almost completely in the lipid bilayer plane. Due the larger dipole moment in certain FFs, an angle of around 90° has been reported for lipid bilayer simulations,[23,47] wheres other reported values have been in within the range of 60–86°.[34,43,112−115] The torsional angles shown in Figure 1 (α5, α4, and α3) are all in gauche conformations which agrees with the experimental findings of Akutsu and Nagamori.[111] This means that the intramolecular interactions in the lipid head group are described well and that its conformation is not fully extended. This might be important when studying membrane-bound proteins and the interactions between small solutes and lipid bilayers since the conformations of the head group alter the dipole moment. When using the popular Berger lipids, the head group is in a more extended conformation which is in disagreement with experiments[111,116]

Figure 4

Average orientation of the P–N dipole with respect to the normal of the membrane.

Other parts of the head group that are important are the charges describing the electron distribution of the glycerol/carbonyl region. Previous studies have shown that the dipole moments of the ester groups correlate with the area per lipid, where a stronger dipole results in a higher area per lipid.[23] A more polar ester group results in a higher water density in this region which in turn increases the area per lipid. In the C36 FF the polarity of this group was increased and it has been proposed that this polarity is one of the reasons for the success of the popular Berger FF.[23] The charges presented in Figure 1 are similar to previously suggested charges.[30,31] Interestingly, the ester oxygen bears almost the same charge in the present FF as the sn-1 ester oxygen in the charge set proposed by Sonne et al.[30]

Area and Volume Per Lipid and Isothermal Area Compressibility Modulus

Area per lipid for DLPC, DMPC, and DPPC at various temperatures is presented in Table 6. The simulations using the FF described here give values that are in excellent agreement with experiments for these lipids. It should be mentioned that experimental data have been reported in a wide range, see, e.g., Table 1 in ref (47). In Table 6 the simulated AL are compared to the most recent experimental data.[117,118] Generally, area per lipid is one of most computed properties from simulated lipid bilayers but the difficulties in obtaining accurate experimental data[10] make these comparisons not conclusive. It has been proposed that the area per lipid alone cannot be an indicator of the quality of a FF[18] and the summary of published experimental and simulation values needs to be used for FF validation.[10,47] Therefore, other properties than just area per lipid must be computed and verified against experiments and only a comparison with an ensemble of different experimental values can really verify the FF and the methodology used. Below, such a comparison for a range of temperatures for the simulated lipid bilayers is presented. Despite the above arguments, area per lipid can still be used to indicate whether or not the lipid bilayer ends up in the correct phase. Incorrect parametrization may lead to underestimated areas indicating that the lipid bilayer is in the Lβ′-phase.[23,30,31] This is not the case in the present investigation, which can be seen from Table 6 showing area per lipid corresponding to Lα-phase for the range of temperatures, in good agreement with experimental results.

Table 6

Average Area Per Lipid (in nm2) Obtained from Simulations at Different Temperatures Compared to Experiments

lipid	temp (k)	simulation	experiment
DLPC	303	0.624 ± 0.004	0.608 ± 0.012[117]
	323	0.646 ± 0.005	0.648 ± 0.013[117]
	333	0.657 ± 0.005	0.659 ± 0.013[117]
DMPC	303	0.608 ± 0.005	0.599 ± 0.012[117]
	323	0.631 ± 0.006	0.633 ± 0.013[117]
	333	0.649 ± 0.005	0.657 ± 0.013[117]
DPPC	323	0.626 ± 0.005	0.631 ± 0.013[117]
	333	0.640 ± 0.005	0.650 ± 0.013[117]
	353	0.672 ± 0.004	0.719[118]

Volume per lipid is a property which can be more rigorously defined, and the interpretation of these properties from experiments is more straightforward compared to area per lipid. In Table 8 these values are compared to experimental data. Values obtained from simulations are systematically underestimated when compared to experiments. It is likely that this lack of consensus is due to the parameters for the head group since the agreement between simulations and experiments for bulk alkanes in Table 4 is very good. This indicates future improvements of the parameter set. The isothermal area compressibility is a property that can be difficult to obtain accurate values for in MD simulations. Often the variance in area per lipid is underestimated resulting in a too high value of KA.[14,45] As can be seen from Table 8, the agreement between simulations and experiments is very good. Since the fluctuations in area per lipid occur on a long time scale,[18,47] relatively long simulations are needed and therefore the simulations were run for 400–500 ns. As a consequence of this, fluctuations on several time scales can occur and contribute to KA. This results in computed values of KA that agree well with experiments. However, in a small system thermal fluctuations could be suppressed resulting in a size dependency[119] since in a larger system undulations make the bilayer more compressible.[45] Poger and Mark studied the impact of system size by using lipid bilayers consisting of 128 and 722 lipids.[47] It was found that KA for both systems converged to the same value, but roughly double the sampling time was required for the smaller system to reach the value of the larger. Possibly, longer simulations could have made the difference between the computed KA and the experimental data even smaller. Based on these conclusions, the system-size dependencies of the systems studied here were not explored.

Table 8

Structural Properties of Lipid Bilayers from Simulations and Experiments: Volume Per Lipid (VL), Isothermal Area Compressibility Modulus (KA), Distance between the the Head Groups (DHH), Luzzati Thickness (DB), and the Hydrophobic Thickness (2DC)

		VL (nm3)		KA (mN m–1)		DHH (nm)		DB (nm)		2DC (nm)
lipid	T (K)	sim	expt	sim	expt	sim	expt	sim	expt	sim	expt
DLPC	303	0.951	0.991[122]	268 ± 24		3.01	3.08[122]	3.04	3.14[122]	2.06	2.09[122]
											2.10[118]
											2.17[117]
	323	0.974		223 ± 31		2.95		3.00	3.10[117]	2.00	2.02[118]
											2.08[117]
	333	0.978		226 ± 33		2.92		2.99	3.07 m	1.98	2.06[117]
DMPC	303	1.060	1.096[123]	250 ± 29[124]	234 ± 23[124]	3.45	3.44[125]	3.53	3.63[122]	2.50	2.54[122]
			1.101[122]		257[122]		3.53[122]		3.67[117]		2.56[118]
									3.69[125]		2.57[117]
	323	1.110		255 ± 32		3.35		3.41	3.52[117]	2.39	2.40[118]
											2.48[117]
	333	1.120		230 ± 20		3.29		3.31	3.42[117]	2.32	2.41[117]
DPPC	323	1.201	1.229[120]	238 ± 35	231 ± 20[8]	3.77	3.80[120]	3.83	3.90[117,120]	2.80	2.84[118,120]
			1.232[8]								2.85[117]
	333	1.209		234 ± 42		3.71		3.70	3.81[117]	2.74	2.72[118]
											2.79[117]
	353	1.223		211 ± 22		3.63		3.59		2.61	2.60[118]

Computed thermal area expansivities are shown in Table 7 and compared to experimental findings. The values obtained from simulations are systematically underestimating this property, which means that the lipid bilayer does not have the same tendency to expand as the temperature is increased. However, the disagreement is not significant and the relative values are reproduced.

Table 7

Thermal Area Expansivity (αAT) and Thermal Contractivities (αDBT, αDCT) from Simulations Compared to Experiments[117] a

	DLPC			DMPC			DPPC
temp (K)	αAT	αDBT	αDCT	αAT	αDBT	αDCT	αAT	αDBT	αDCT
303
sim	1.8	0.6	1.3	2.2	2.0	2.4	–	–	–
expt	2.8	1.9	1.6	3.2	2.2	2.0	–	–	–
323
sim	1.7	0.6	1.4	2.1	2.1	2.5	2.5	2.0	2.3
expt	2.6	2.0	1.7	3.0	2.3	2.1	3.0	2.3	2.1
333
sim	1.7	0.6	1.4	2.1	2.2	2.6	2.4	2.1	2.3
expt	2.6	2.0	1.7	2.9	2.4	2.2	2.9	2.4	2.2

Reported values are expressed in the unit 103 K–1.

Bilayer Thickness

Three types of membrane thicknesses were analyzed: the distance between the phosphate peaks in the electron density, the Luzzati thickness[8] (DB), and the hydrophobic thickness (2DC). DB is calculated according to[22,120]where d is the repeat z-spacing along the the bilayer normal of the simulation box, i.e., the z-dimension of the simulation box. ρw(z) is the probability distribution of water in the z-direction and is calculated from the time-averaged histograms of the water distribution along the z-axis with a thickness of dzwhere nw is the time average of water molecules per slice and dV is the time-averaged volume of each slice. 2DC is computed in a similar way to DBwhere ρCH(z) = ρCH + ρCH.[22,120] The probability distributions are calculated according toAs can been seen from Table 8, the agreement between simulations and experiments is excellent for the phosphate–phosphate distance for all three lipids. Consistent with previous work, the Luzzati thickness is slightly underestimated.[21,22,47] Since the DB is a measure of how deep water is allowed to penetrate the lipid bilayer, it is important to have good agreement between the simulation and experiments in order to reassure that the models used are describing the physics of the system in a correct manner. The polarity of the ester group will determine how deep the water is willing enter the lipid bilayer in order to form hydrogen bonds with the carbonyl dipole. As discussed previously, this is strongly correlated with the area per lipid. The hydrophobic thickness is an important property of a lipid bilayer since this property has been proposed to dictate a number of cellular processes via, e.g., hydrophobic mismatch.[121] The differences between simulations and experiments for 2DC are small for all temperatures, which is consistent with the computed values for DB.

The decrease in bilayer thickness with increasing temperature can be explained by the increase in chain distortion.[117] Due to the pronounced thermal dependency of the trans–gauche isomerization, this phenomenon gives rise to the major part of the differences observed. As the chain disorder increases, the lipid bilayer becomes more compact in the direction of the membrane normal (the z-direction). In an expansion of the lipid bilayer in the xy-plane, the area is increased and the thickness has to be decreased. In Table 9 the number of gauche bonds for the DPPC lipid simulated in a bilayer at various temperatures are presented.

Table 9

Number of Gauche Bonds per DPPC Lipid Molecule from Simulations at Various Temperatures Compared to Experimental Estimates

	293 K	323 K	333 K	353 K
simulation	0.8	3.9	4.0	4.2
experiement	<1.0 (300 K)[133]	3.6–4.2[130]
		3.7[131]
		3.8[126]

Thermal contractivities from simulations are compared to experiments in Table 7. The small temperature dependency is captured by the simulations and for DMPC and DPPC the hydrophobic thickness thermal contractivity is slightly overestimated whereas for DLPC it is underestimated to a small degree. The agreement for the bilayer thickness contractivity between simulations and experiments for the former two lipids is better.

Taking all these properties into consideration, we can conclude that the current FF balance the hydrophilic interactions of the lipids well with the hydrophobic effect. Water is not allowed to penetrate the structure too deeply, and the hydrophobic core of the membrane is kept intact, which is important. This might be shown important in future simulations of realistic models of biological membranes in the presence of proteins and membrane-associated molecules.

Ordering and Conformations of Lipid Tails

The deuterium order parameter is a property that is determined from experiments without model assumptions, and the methods used are robust which results in measured vales that are reproducible. As a result of this, |SCD| is one of the most vital properties of a set of FF parameters for which to get accurate values. Computed deuterium order parameters (|SCD|) together with experimental data for the aliphatic tails of DLPC, DMPC, and DPPC are presented in Figure 5. Values obtained from simulations agree very well with experiments.[118,126,127] The splitting of order parameter for the second carbon in the sn-2 tail observed in experiments is also present in the simulations for each lipid type, however, not to the same extent. Experimental studies have shown that the upper part of the two lipid tails are aligned differently with the sn-1 tail in a more perpendicular position relative the bilayer normal than the sn-2 tail.[128,129] As can be seen in Figure 5, this difference is also present in the simulations for all lipids. The agreement in order parameters indicates that the simulated lipid bilayers have the correct internal (disordered) structure and that the bilayers are all in the Lα-phase. A large discrepancy between experimental findings and simulations where |SCD| would have been overestimated by the latter would indicate a lipid bilayer in a too ordered phase, i.e., the Lβ′-phase.

Figure 5

Deuterium order parameter, |SCD|, from simulations and experiments[118,126,127] for the lipid tails of DLPC (top), DMPC (middle), and DPPC (bottom) at 303, 303, and 323 K, respectively.

The number of gauche bonds in the aliphatic chains of the lipids can be computed in order to further test if the lipid bilayer is in the correct phase during the simulations. In Table 9 the number of gauche bonds from simulations of a DPPC lipid bilayer is compared to experiments.[130,131] Simulations performed at 323 K agree well with the experimental estimates, and this together with the calculated deuterium order parameter and area per lipid indicates that the lipid bilayer is in the correct phase. A higher percentage of trans conformations (and low number of gauche bonds) induces order in the system resulting in a gel phase.[132] In Figure 6 three simulation snapshots are shown at different temperatures where the increase in gauche bonds per lipid molecule and general disorder with increasing temperature are well illustrated.

Figure 6

Snapshots from simulations of a DPPC lipid bilayer at different temperatures: 293 K (left), 323 K (center), and 353 K (right). The choline groups are rendered in blue, the phosphorus groups in cyan, the glycerol/carbonyl groups in red, and the hydrophobic tails in white.

Electron Density Profiles and Scattering Form Factors

In MD simulations it is straightforward to obtain electron density profiles and it is possible to divide the total electron density into the constituting groups of the lipid molecule as defined in Figure 1. The overall electron densities for the studied lipids are presented in Figure 7 together with the contributions from each group. From these plots the head-to-head distance (DHH) can be obtained (see previous sections). Other valuable information such as the fact that water molecules penetrate the lipid bilayer down to the carbonyls of the aliphatic tails can also be extracted, which is in agreement with experiments.[134] Furthermore, it is evident that water does not penetrate the hydrophobic core of the lipid bilayer, in agreement with experimental findings.[122,135−137] For all lipids, the leaflets are interdigitated[138] and the peaks corresponding to the head groups are broad, indicating that they are all in the Lα-phase since in gel phase the leaflets are less interdigitated and the head group’s electron density peaks are divided into two sharp peaks instead of one broad peak.[139,140]

One inherent problem with experimental electron density curves is that they are all based on certain assumptions and models which makes direct comparison between simulations and experiments difficult.[48,117,141−143] For this reason, the densities in Figure 7 were not compared to experiments. Electron density profiles are obtained from an inverse Fourier transformation of the scattering form factors. A better way of comparing simulations and experiments is to compute the form factors from the total electron density profiles obtained from the simulations. These form factors can then be directly compared with the experimental data. This comparison together with the order parameters should serve as the major part of the testing procedure when the quality of a FF is established. F(q) and |SCD| (at least for the lipid tail) can be obtained from experiments in a unambiguous manner and are reproducible and therefore should receive the most focus.[48,117,143] In Figure 8 form factors obtained from MD simulations are compared to experimental form factors[117,120,122] at various temperatures. The most common temperatures to perform experiments and simulations for DLPC, DMPC, and DPPC are 303 and 323 K. For these temperatures, the agreement between simulations and experiments is very good, providing a clear proof that the FF presented here gives the correct structure for a fully hydrated lipid bilayer. The position where the form factor goes to zero is reproduced in an excellent manner, and the relative peak heights agree very well with the available experimental data. A direct indication of this is that all the lipid bilayers have the correct thickness and structure. A small offset can be seen for the third peak of DMPC, but the difference is very small. When the temperature is increased, the agreement is slightly decreased. Still, the simulations give more than reasonable agreement with experiments. The largest discrepancy with experimental data is found for DLPC, for the third peak, when the temperature is increased from 303 to 323 and 333 K. The insets in Figure 8 show the neutron scattering FFs obtained from simulations compared to experiments.[117] The SDP model introduced by Kučerka et al.[120] was used to compute the FFs. As can be seen, the agreement is good and the points where the form factors have their minima are reproduced by the simulations. These findings extend the evidence that the FF presented here is highly suitable for simulating lipid bilayers.

Figure 8

Form factors for the simulated lipid bilayers (lines) compared to available experimental data (points)[117,120,122] at various temperatures. In the insets, neutron scattering form factors are shown. Experimental F(0) values are shown as black squares.[122]

Lateral Diffusion

In membranes the lipids are free to diffuse in the plane that is perpendicular to the bilayer normal. To test the parameter set proposed here, lateral diffusion coefficients, Dl, were therefore computed. The Einstein relation was used to calculate Dl via the mean-squared displacement (MSD). In order to correct for the artificial center of mass motion of each monolayer,[18,144] this contribution was subtracted during the course of the MSD calculations. The diffusion of lipids in a bilayer occurs on two different time scales: on a short time scale (<5 ns) where the center of mass of the lipid is in motion due to the conformational freedom of the acyl chains and the head group is relatively fixed in space,[145] and on the longer time scale the whole lipids move in the monolayer more freely. In the present investigation, only the long time diffusion was investigated. The linear fitting was performed from 50 to 350/450 ns. The final results are shown in Table 10. For DLPC, DMPC, and DPPC Dl is underestimated with respect to the experimental data.[146−151] The values obtained for DMPC at 303 K are close to other values obtained from MD simulations.[145,152,153] The drastic decrease in Dl observed for DPPC when the temperature is lowered is expected due to the phase transition. As in the latter case the lipids are in a gel phase, the configurational freedom for the acyl tails is very limited, resulting in a very low diffusion coefficient. A large drop in lateral diffusion (about 2 orders of magnitude) when going from the Lβ′-phase to the Lα-phase has been reported earlier by Marrink and co-workers in coarse-grained simulations.[154] The increase in the lateral mobility of the lipids with an increase in temperature can be observed which is related to the area per lipid. As the temperature is increased, the area per lipid is increased, making it easier for the lipid to move in the two-dimensional liquid. For DMPC, the lateral diffusion coefficients obtained from simulations when the temperature is increased are about half of the experimental values as can be seen in the fifth and sixth columns of Table 10. In order to enhance Dl accelerated MD simulations[155] could be applied where the sampling of the conformational space is improved. This approach has been used earlier in lipid bilayer simulations.[153] One possible explanation for the systematic underestimation of the lateral diffusion could be the system size. In a lipid bilayer patch of 128 lipids, no efficient collective movement can occur that could result in a decrease in Dl. This is similar to the case of bulk alkanes when the diffusion coefficients computed in a box of limited size were underestimated, but after correction on collective motion[76] they occur in very good agreement with experimental results (Table 4). The agreement between simulations and experiments for the normalized values (with respect to the highest temperature) in Table 10 supports this, indication that there could be a systematic underestimation which could originate from this effect.

Table 10

Lateral Diffusion Coefficients for DLPC, DMPC, and DPPC in a Range of Temperatures Compared to Experimental Findings

		Dl (10–8cm2 s–1)		normalized Dl
lipid	temp (K)	sim	expt	sim	expt
DLPC	303	6.66 ± 0.45	8.5[146],a	0.249
	323	14.3 ± 3.50		0.536
	333	26.7 ± 4.70		1.000
DMPC	303	5.22 ± 0.49	5.95,[147] 9.00[148,151]	0.332	0.291
	323	11.8 ± 2.10	22.3[151]	0.752	0.722
	333	15.7 ± 3.10	30.9[151]	1.000	1.000
DPPC	293	0.15 ± 0.08		0.005
	323	9.84 ± 1.30	12.5,[149] 15.2[150]	0.316
	333	13.6 ± 0.90	9.7[156]	0.437
	353	31.1 ± 3.50		1.000

298 K.

Fully Hydrated Gel Phase

To further test the FFs ability to describe the nature of a lipid bilayer under various thermal conditions, simulations of a DPPC bilayer in gel phase were conducted. The phase transition from the Lα-phase to the Lβ′-phase occurs at 314 K for DPPC,[157] so the simulations were performed at a temperature of 293 K in order to ensure that the lipid bilayer entered the correct phase. In order to prepare the gel phase, the size of the carbon atoms for the special 1-4 LJ interactions was increased so that the bilayers entered the gel phase without any major problems of metastable configurations. For the production run, the 1-4 LJ interactions were set to the original values with a equilibration time of 40 ns and a production simulation of 300 ns with the same procedure as for the bilayers in the Lα-phase. This approach to force a lipid bilayer into the Lβ′-phase is similar to the methodology described by Schubert et al.[158] Simulations where the temperature was slowly decreased to 293 K without the biased 1-4 interactions resulted in a lipid bilayer in the rippled phase (Pβ′). As this is a metastable state at 293 K the system would eventually reach the more stable Lβ′-phase. Since the dynamics of the system is slowed down by roughly 2 orders of magnitude, biased intramolecular interactions were introduced merely to force the lipid bilayer into the correct phase on a shorter time scale. Thus, we can conclude that the correct phase behavior is to be expected. We believe that if the parameter set would have given rise to a more fluid phase than Lβ′ at 293 K this would have been observed during the 300 ns production simulation since melting is a process that generally occurs on a smaller time scale than freezing. In Figure 6 the obtained gel phase, Lβ′, is illustrated where the lipids in each monolayer are parallel.

In Table 11 the results of the simulation are summarized. It is evident that the bilayer entered the gel phase and that the hydration number (30 water molecules per lipid) did not play a significant role in accordance with the findings of Marrink et al.[154] The snapshot from the simulation presented in Figure 6 verifies this graphically. The area per lipid is overestimated by roughly 0.02 nm–2, which is the only discrepancy in the comparison between simulations and experiment. The tilt of the lipids with respect to the normal of the bilayer is characteristic of the Lβ′-phase and the tilt angle obtained from simulations is in very good agreement with experimental data.[159,160] The head-to-head distance also agrees well with experimental findings.[159] The significant decrease in number of gauche defect per chain (Table 5) illustrates that the bilayer has entered the gel phase and that the conformations of the chains is of importance in controlling the bilayers fluidity.

Table 11

Structural Properties from Simulations Compared to Experiments for the Fully Hydrated Gel Phase of DPPC

	AL (nm–2)	DHH (nm)	title angle (deg)
sim (293 K)	0.492 ± 0.006	4.37	33.6
expt (292 K)[159]	0.472 ± 0.005		32.0 ± 0.5
expt (298 K)[160]	0.473 ± 0.003	4.28 ± 0.02	31.6 ± 0.4

These findings strengthen the present FFs ability to model lipid bilayers under a range of temperatures and shows that the correct phase behavior can be expected. There is a growing interest in simulating membranes in other phases than the liquid crystalline one and the data presented here makes this FF a good candidate for these simulations.

Conclusion

It is critical for computer simulations of lipid bilayers to have sets of parameters obtained in a clearly defined and physically grounded procedure, in order to produce reliable and meaningful results. Here, the development of new AA lipid FF has resulted in a consistent FF where the hydrophilic and hydrophobic forces balance each other, which is crucial. This balance is what drives the formation of the different structures that can be formed by lipids. As model compounds for the aliphatic chains of saturated lipids, a series of alkanes was chosen and many experimentally accessible properties for these bulk liquids were reproduced in the simulations. Furthermore, a large number of important properties for lipid bilayer are in very good agreement: area and volume per lipid, membrane thicknesses, isothermal area compressibility, and trans/gauche ratio. But, most importantly, the key properties, NMR order parameters and scattering form factors, are excellently reproduced which is a clear proof of the abilities of the current FF. Simulations performed with DPPC in the gel phase show that the FF gives the correct phase behavior as well. To our knowledge, this is the first reported set of parameters that reproduces all these properties consistently under a range of temperatures. The ability of the FF to reproduce a lipid bilayer in a liquid-crystalline phase in a tensionless ensemble and describe the physics of lipid in a membrane accurately from more than just the point of view of one or two properties makes it an excellent choice for more advanced simulations, including several lipid types and important membrane-inserted molecules. Most importantly, |SCD| and F(q), which are properties that can be determined without any empirical input from an experimental point view, are in very good agreement with available experimental data.

As protein simulations in membranes become more popular and since the FF presented here is compatible with the AMBER, FFs for biomolecule simulations can easily be extended to include transmembrane proteins. Ongoing work in our group aims to expand the FF to unsaturated and polyunsaturated lipids, lipids with different head groups, and cholesterol. A FF that has been developed in a consistent manner for a large number of lipids will be important in future modeling of biological membranes.