Compatibility of advanced water models
with a united atom model of lipid in lipid
bilayer simulation
Cite as: J. Chem. Phys. 151, 065104 (2019); https://doi.org/10.1063/1.5108830
Submitted: 02 May 2019 . Accepted: 18 July 2019 . Published Online: 13 August 2019
Shakkira Erimban, and Snehasis Daschakraborty

J. Chem. Phys. 151, 065104 (2019); https://doi.org/10.1063/1.5108830© 2019 Author(s).

151, 065104

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

Compatibility of advanced water models
with a united atom model of lipid in lipid
bilayer simulation
Cite as: J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Submitted: 2 May 2019 • Accepted: 18 July 2019 •Published Online: 13 August 2019
Shakkira Erimban and Snehasis Daschakrabortya)
AFFILIATIONS
Department of Chemistry, Indian Institute of Technology Patna, Bihar 801106, Indiaa)Author to whom correspondence should be addressed: snehasis@iitp.ac.in

ABSTRACT

Molecular dynamics simulation of lipid bilayers generally uses all-atom, united-atom, and coarse-grained models of lipid molecules. The
GROMOS united-atom model of lipid constructs a balance between accuracy and computational cost. The above-mentioned model satisfactorily reproduces many of the structural and dynamical properties of different lipid bilayers. However, the GROMOS force field is parameterized
only with the SPC model of water. Unfortunately, SPC is not an excellent model of water for predicting the structure and dynamics of theinterfacial water near the lipid bilayer. More advanced water models, such as TIP3P-FB and TIP4P-FB, outperform the SPC model in predicting different thermodynamic and microscopic properties of bulk water. This motivates us to check the compatibility of five different watermodels, including SPC, with the GROMOS96 53A6L united atom model of two different lipid bilayers, DPPC and POPC. A systematic comparison of the bilayer structure and dynamics, resulting from the simulations with different water models, is done. We find general agreementof the results for different water models with the experiment. In fact, the more advanced water models provide better agreement with theexperiment. This study, therefore, helps in widening the range of choice of water models in simulating the lipid bilayer using the GROMOS
united-atom model of the lipid molecules.
Published under license by AIP Publishing. https://doi.org/10.1063/1.5108830., s

I. INTRODUCTION
Molecular dynamics (MD) simulation studies have immensely
contributed to understanding the structural and dynamical properties of lipid bilayers. This has been made possible by the continuous
development of more and more accurate force fields of various lipid
molecules. Three major classes of force fields are used in the MD
simulation of the lipid bilayers. They are all-atom, united atom, andcoarse-grained models. CHARMM,1–5 AMBER,6–11 OPLS-AA,12–14
and SLIPIDS15,16 are three major all-atom force fields of lipid, among
which the CHARMM force field has been extensively used in many
MD simulations. Being all-atom, these models accurately describe
the structure of a lipid molecule. However, the agreement of different simulated parameters with experimental quantities is not always
seen.17 The coarse-grained model of the lipid molecule, such as
MARTINI force field,18,19 considers mainly four-to-one mapping,where four heavy atoms are represented by only one coarse-grained

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishingsite. This allows choosing a longer simulation time (microsecond)
and a much bigger system size (millions of lipid molecules) thanthose in all-atom models. The coarse-grained model can be used
to study much more complex phenomena, such as undulations ofmembrane surfaces, formation of aggregates, phase transformation,
behavior of membrane proteins, and lipid-protein interaction.20
However, the use of this model is associated with several limitations, such as freezing of the coarse-grained lipid model at normal
temperature.18,19
Apart from the all-atom and the coarse-grained model, there
is also the united-atom model of lipid, which constructs a balance
between the complexity of the lipid molecules and the computational overhead. GROMOS21,22 is a computationally efficient unitedatom force field, which considers nonpolar CH, CH2 , and CH3 assingle carbon atoms. There are mainly two classes of force fieldparameters: one with original GROMOS parameters21–23 and theother with Berger modifications.24 Major versions of the GROMOS

151, 065104-1

The Journalof Chemical Physics

force field are GROMOS8725 and GROMOS96.26 The GROMOS
force field reproduces the properties of the system almost accurately. The latest version of GROMOS96 includes 53A5 and 53A6.27
GROMOS96 53A6L ,28 a modified version of 53A6, attempts to maintain the fluidity of the membrane, which was lacking in the 53A6
version.
The structural and dynamical properties of the lipid bilayer
are sensitive to the interaction between the hydrophilic head groupof lipids with water. Use of an arbitrary water model in a lipid
bilayer simulation is generally not recommended. In principle, one
should use that particular water model, with which the lipid modelis parameterized. For example, the CHARMM force field of lipid
was initially been parameterized with the TIP3P model29 of waterand later reparameterized with a flexible SPC30 model of water.31
However, recently,32,33 it has been found that the CHARMM36 force
field of lipid is also compatible with the TIP4P/2005, TIP3P-FB,
and TIP4P-FB water models34 as the above combinations of force
field do not lead to severe flaws in the structure and dynamics ofthe lipid bilayer. Unfortunately, the GROMOS96 53A6L model oflipid has been parameterized only with the SPC model of water.28
Even though the above combination of force fields reproduces theexperimentally observed properties with reasonable accuracy, it isstill essential to work with more advanced water models—such as
TIP4P/200535 and TIP4P-FB—mainly when a study focuses on the
structure and dynamics of water at the lipid bilayer interface. Compatibility of water models, such as TIP4P/2005 and TIP4P-FB, would
also allow one to simulate the lipid bilayer at low temperatures where
the SPC model does not reproduce the experimental properties ofliquid water as well as the former models. Note that the MD simulation study of lipid membrane at low temperatures is significant
as it may shed some light on the mechanism of cold adaptation of
psychrophiles, which is a class of extremophilic organisms capableof growth and reproduction at low-temperature conditions (−20 ○ C
to +10 ○ C).36
Widely used three site water models, TIP3P and SPC, reproduce many properties of bulk water under normal conditions.37
However, the self-diffusion coefficient and static dielectric constant
are poorly described.38 In addition, the solvation entropies are notwell reproduced by these water models.39 The freezing temperatures

ARTICLE

scitation.org/journal/jcp

of these models are less than 200 K.40 The density maximum arises
at a much lower temperature for these two models.41 Four-pointwater models perform better than the three- and five-point models
in reproducing many of the thermodynamic properties of water at
a broader range of temperatures. TIP4P/2005 accurately reproduces
the density of water at a wide range of temperatures (240 K–370 K).35
However, it deviates from the observed values for some other properties, such as dielectric constant34 . The relatively recently developed
water models, TIP3P-FB and TIP4P-FB, are parameterized to well
reproduce the enthalpy of vaporization, density, dielectric constant,isothermal compressibility, heat capacity, and the thermal expansion
coefficient of liquid water.34 They are found to describe many of the
key physical properties of water, such as the dielectric constant andthe viscosity, with very good accuracy.34 Some of these parameters
are listed in Table S1 of the supplementary material for different
water models.
The above discussion may motivate one to reparameterize theexisting model of the lipid membrane, which can be compatible
with more advanced models of water. However, it is far easier tocheck the compatibility between the GROMO96 53A6L and different water models in the first place. Here, our goal is to checkthe reproducibility of a large number of experimentally observed
quantities using the GROMO96 53A6L 28 force field for two lipids
POPC and DPPC, in combination with five different models ofwater: SPC, TIP3P, TIP4P/2005, TIP3P-FB, and TIP4P-FB. Note,
POPC contains a saturated sn1 chain and an unsaturated sn2 chain,
and it is observed in most naturally occurring lipid membranes(see Fig. 1 for the labeling of the chains). It is mostly present
in the outer leaflet where it functions as part of the permeability barrier. DPPC contains two saturated chains and is a welldocumented lipid. This study may guide the best possible watermodel(s) to consider with the GROMOS G53A6L force field. In addition, the present study may indicate the requirement of reparameterization of the lipid force field to combine with a particular watermodel.
The following is the outline for the remainder of this paper.
The molecular models and simulation details are detailed in Sec. II.
Section III details the simulation results and discussions. Concluding
remarks are offered in Sec. IV.

FIG. 1. Scheme for the two lipid
molecules with the relevant atoms andthe two chains labeled.

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-2

The Journalof Chemical Physics

ARTICLE

II. SIMULATION DETAILSEach of the two bilayers, POPC and DPPC, is simulated withfive different water models: SPC, TIP4P/2005, TIP3P, TIP3P-FB,
and TIP4P-FB. These aggregate a total of ten different simulation
systems. The lipid molecules are modeled as a united-atom modelwith GROMOS96 G53A6L 28 force field parameters.
Initial configurations of the lipid bilayers are built using the
PACKMOL software.42 Each system includes 200 lipids with 7000
water molecules in a simulation box of dimension 8 nm × 8 nm
× 8.5 nm. The above composition ensures a well-hydrated lipid
bilayer structure with the hydration number of 35. The above areaof the xy plane (64 nm2 ) at the initial stage is set close to the experimentally observed area per lipid 0.643 nm2 for the POPC43 and0.631 nm2 for the DPPC lipid bilayer.44 The MD simulation for thesesystems is done using the GROMACS software package.45 First, we
use the steepest-descent algorithm for energy minimization of thesesystems. The bonds are constrained with the LINCS algorithm46 forboth lipid and water. The systems are then equilibrated using the
NPT ensemble for 50 ns. The temperature is maintained at 300 K
for POPC and 323 K for the DPPC lipid layer—which is above theirrespective melting temperatures (for POPC T m = 270.5 K47 and for
DPPC T m = 314 K48 )—using the Nose-Hoover thermostat49,50 (coupling constant 0.1 ps). The pressure is fixed at 1 bar in the lateral andnormal directions (semi-isotropically) using the Parrinello-Rahman
barostat.51 The semi-isotropic pressure bath has an isothermal compressibility of 4.6 × 10−5 bars−1 and a coupling constant of 1 ps.
These simulation parameters were also used elsewhere.28 The equations of motion are always solved using the Verlet leapfrog algorithm
every after 2 fs. The Lennard-Jones interactions are scaled to zeroafter a distance of 14 Å for SPC and 20 Å for all other water models. The cutoff distance for short-range force is taken as 8 Å for
SPC and 21 Å for other models. Long-range electrostatic interactions are handled by the Particle Mesh Ewald (PME) summation
method. The interactions are updated every 5 simulations steps.
The equilibration of each system up to 50 ns time is followed by
the final production run by the NPT ensemble for another 150 ns,
during which we save the configurations every after 250 fs for the
analyses.

A. Structural properties
1. Average area per lipid
The average area per lipid (⟨AL ⟩) is calculated by dividing
the time-averaged lateral dimension of the simulation box by the
number of lipid molecules present in each leaflet nl ,(1)

Some studies also probe the curved area per lipid, which takes
into account the curved lipid membrane surface. However, a recent
simulation study17 has indicated that for PC lipids, such as DPPC
and POPC, the standard area per lipid shows closer agreement withthe experimental value than the curved area. Therefore, we avoid
the calculation of the curved area per lipid. Note that the reliability of the experimental ⟨AL ⟩ values is always questionable. The⟨AL ⟩ value of a lipid bilayer cannot be measured directly from theJ. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishingexperiment. In addition, the thermal fluctuations of the fluid phaseof the lipid bilayer do not allow an accurate determination of the⟨AL ⟩ values from the experimental studies.52,53 Due to the above fact,
a wide range of ⟨AL ⟩ values for a particular lipid bilayer is reported
in the literature.54 Although a direct comparison with experimental results is not very much reliable, it is always important for asimulation study to predict the ⟨AL ⟩ value in agreement with theexperimental range.
The simulated ⟨AL ⟩ values are reported in Table I for different models of water and compared with the experimental43,55–57 andsome previously reported simulated data.24,58,59 The time dependencies of ⟨Axy ⟩ for all ten systems are shown in Fig. S1 of the
supplementary material, where we see the convergence of the lateraldimension of the box after 50 ns simulation time. General agreement
between the simulated ⟨AL ⟩ values for different water models andthe experimental values is evident in Table I. While the TIP4P/2005
water model gives relatively lower values of ⟨AL ⟩, a closer agreement is observed for the other water models, TIP3P, TIP3PFB, and TIP4P-FB, between the simulated and experimental ⟨AL ⟩values55 for both the lipid membranes. The value of ⟨AL ⟩ is lowest in the case of TIP4P/2005 water model probably because of thelowest dielectric constant (∼60) of the latter water model comparedto those of all other models. The dielectric constant values for different water models are listed in Table S1 of the supplementarymaterial.
2. Isothermal area compressibility modulus
The isothermal area compressibility modulus K A is calculated
from the fluctuation of average area per lipid in an equilibrium
simulation by the following equation:
KA =kB T⟨AL ⟩.
nL ⟨δA2L ⟩(2)Here, kB is the Boltzmann constant, T is the average temperature,⟨AL ⟩ is the average area per lipid, and δA2L is the variance associated with AL . The simulated average area compressibility moduli are

TABLE I. The simulated ⟨AL ⟩ values of the two lipid bilayers—each of which is simulated with different water models—and their comparison with earlier simulated andexperimental data. The values in the parentheses are the standard error calculatedby block averaging.

III. RESULTS AND DISCUSSION⟨AL ⟩ = ⟨Axy ⟩/nl .

scitation.org/journal/jcp

Lipid

POPC

DPPC

Water model⟨AL ⟩ (nm2 )

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

0.61(0.003)
0.66(0.003)
0.67(0.002)
0.60(0.002)
0.68(0.002)

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

0.60(0.002)
0.65(0.001)
0.67(0.003)
0.57(0.002)
0.67(0.002)

Literature (nm2 )Expt.:43,55 0.630–0.683
Simulation:58 0.655

Expt.:56,57 0.57–0.717
Simulation:24,59 0.61–0.637

151, 065104-3

The Journalof Chemical Physics

ARTICLE

given in Table II. The simulated values are consistent with the experimental values for both the POPC and DPPC except when SPC watermodel is used. This is because ⟨AL ⟩ is also not correctly predicted by
the SPC water model. The higher values of K A , predicted by the SPC
water, were also reported in a previous study.17 Therefore, the agreement between the experimental60,61 and simulated17,62,63 K A values
is significantly improved by the use of more advanced water models.
However, it is also true that the accuracy of the experimental K A values is limited by the accuracy of the measurement of ⟨AL ⟩ value inthe experiment.

TABLE III. Average volume per lipid molecule V L for DPPC and POPC lipid bilayers
simulated with different water models and the comparison between the present resultswith earlier simulated and experimental data. The values in the parentheses are thestandard error calculated by block averaging.

SystemWater model

V L (nm3 )Literature (nm3 )

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

1.23(0.0003)
1.20(0.0002)
1.20(0.0005)
1.26(0.0003)
1.20(0.0001)

Expt.:43,74 1.223–1.256
Simulation:28 1.232

DPPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

1.20(0.0003)
1.17(0.0002)
1.18(0.0003)
1.20(0.0004)
1.18(0.0002)

Expt.:43,61 1.229–1.232
Simulation:28 1.226

3. Average volume per lipid
The volume occupied by lipid molecules is one of the essential
parameters characterizing the lipid bilayer structure. In the simulation, the following equation calculates the volume per lipid molecule
V L:
VL = (V − nW VW )/nL .

(3)Here, nL and nw are the numbers of lipids and water molecules,
respectively. V w is the volume per water molecule, which has been
determined from a 15 ns long independent simulation of pure bulk
water. V L can be measured by different experimental techniques,
such as neutral buoyancy,64–66 dilatometry,67,68 pyknometry,69 andvibrating tube densitometry.70,71 Unfortunately, any of these techniques do not provide very accurate V L values.72 Therefore, MD
simulations73 better predict the molecular volume and overcome theerror bars associated with the experimental estimation of V L . Simulated V L for ten different systems are reported in Table III. In thecase of POPC, water models do not significantly influence the simulated values, which are in good agreement with the experimental
values.43,61,74 However, we see some deviation from the experimental
values for the DPPC lipid bilayer.
4. Bilayer electron density and scattering form factors
Detailed structural information of the lipid bilayer is obtained
from electron density profiles (EDPs) along the bilayer normal.
However, the experiment cannot measure EDPs directly. Therefore, we Fourier transform the simulated electron and neutron
TABLE II. The simulated K A values of the two lipid bilayers—each of which is simulated with different water models—and their comparison with earlier simulated andexperimental data. The values in the parentheses are the standard error calculatedby block averaging.

Lipid

Water model

K A (N/m)

Literature (N/m)

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

0.54(0.08)
0.23(0.01)
0.35(0.03)
0.34(0.07)
0.31(0.04)

Expt.:60 0.18–0.33
Simulation:62,63 0.2–0.6

DPPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

0.48(0.05)
0.28(0.01)
0.27(0.03)
0.23(0.07)
0.28(0.03)

Expt.:61 0.231
Simulation:62,63 0.2–0.6

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishingscitation.org/journal/jcp

density profile to obtain the simulated X-ray and neutron scattering form factors, which are directly compared with those from theexperiment.
Figures 2(a) and 2(b) present the total EDPs for DPPC and
POPC bilayers, each with five different water models. The EDPs forwater are also exhibited in the same figures. The EDPs for individual
groups are shown in Fig. S2 of the supplementary material. Comparison of the EDP for the total system with those of the individual
groups reveals the origin of different peaks and dips of the total EDPs
for the overall system. The nearly symmetrical profiles for all thesystems indicate proper equilibration of the lipid bilayers. The twomost pronounced peaks of the EDPs correspond to the phosphategroup in the head group region of lipid. Water molecules penetrate
the bilayers up to the CO groups, but the terminal methyl groups in
the acyl chains stay dehydrated. This is in agreement with the experiment.43 A dip is observed in the middle of the EDPs, which indicatesa reduction in the interaction between the lipid layers at the center.
The EDPs for different water models are overall consistent with eachother. However, we see that the dip of the bilayer center is slightly
more intense in the case of SPC and TIP4P/2005 water models than
the remaining three.
The lipid bilayer thickness is an important quantity to identify the flexibility and packing of the lipid bilayer. There are multiple interpretations for the bilayer thickness. One of them is the
distance between the two phosphate peaks DHH in EDPs. Experimental scattering form factors can provide the DHH value61 only
indirectly. This is because the scattering form factors only probe thetransverse lipid organization of the bilayer. Therefore, the Fouriertransformed electron density profiles do not accurately give the headgroup-head group distance for the calculation of bilayer distance
DHH .61 However, being motivated by many simulation works,77–81
we have compared the obtained thickness with the experimentally
reported values. The simulated DHH values of the two lipid bilayers,
simulated with different water models, are presented in Table IV. All
water models result in the bilayer thickness almost consistent withthe experimental values.43,61,75
Now, we compare the simulated structures of the lipid membranes with those from the experiment.43 The X-ray and neutron

151, 065104-4

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

FIG. 2. Electron densities [(a) and (b)],X-ray [(c) and (d)], and neutron [inset of(c) and (d)] scattering form factors for
POPC [(a) and (c)] and DPPC [(b) and(d)] bilayers. Total electron densities areshown in solid lines, while the electron
densities for water are shown in broken lines. The results for different watermodels are shown in color-coded lines(SPC: red; TIP3P: blue; TIP4P/2005:
pink; TIP3P-FB: dark red; and TIP4P-FB:
dark gray). The experimental X-ray andneutron scattering form factor (for different D2 O concentration) data43 are presented in open circles.

scattering experiments probe the lipid bilayer structure by scatteringform factors, which determine the scattering length density along thebilayer normal. The form factor is connected to the atom number
density through the following equation:
F(q) = (∫

D/2

−D/2

∑α fα (qz )nα (z) − ρS ) exp(izqz )dz.

F(q) = ∫(4)Here, nα (z) is the number density of the αth atom as a function of the
distance from the bilayer center across the bilayer normal. f α (qz ) is
TABLE IV. The average thickness of the DPPC and POPC lipid bilayers simulatedwith different water models and the comparison between the present results with earlier simulated and experimental data. The values in the parentheses are the standard
error calculated by block averaging.

Lipid

Water model

DHH (nm)

Literature (nm)

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

3.42(0.04)
3.46(0.04)
3.41(0.02)
3.47(0.01)
3.35(0.02)

Expt.:43 3.70
Simulation:28 3.46

DPPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

3.59(0.02)
3.43(0.04)
3.34(0.01)
3.52(0.03)
3.44(0.05)

Expt.:61,75 3.42–3.83
Simulation:28 3.57

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

the scattering length density, ρS is the solvent scattering length density, and D is the bilayer thickness. For a symmetric bilayer, Eq. (4)reduces to the following form:
D/2

−D/2

Δρe (z) cos(zqz )dz.

(5)Here, Δρe (z) is the scattering length density difference between
the solvent and the bilayer. We calculate the atom number density
profile from the simulations and determine the form factor using
the SIMtoEXP software.76 X-ray scattering form factors for the lipid
bilayer systems are studied in water, while the neutron scatteringform factors are studied in three different concentrations of D2 O
(100%, 70%, and 50%) for POPC and two different D2 O (100%and 50%) concentrations for the DPPC lipid bilayer. The simulated
and experimental form factors for different systems are shown in
Figs. 2(c) and 2(d). Usually, qualitative comparisons are made withthe experimentally obtained scattering profiles. Some quantitative
comparisons, such as the peak heights,16,80,81 the minima and maxima positions,78 and the root mean square deviation (RMSD) ora reduced χ 2 value,76,82 are also made between the simulated andmeasured scattering profiles.
Overall, we see a good agreement with the experiment for boththe DPPC and POPC lipid bilayers simulated with all five watermodels. A closer comparison reveals that the TIP3P, TIP3P-FB,
and TIP4P-FB water models reproduce the positions and the heights
of all three experimental lobes much better than the SPC and
TIP4P/2005 water models. Importantly, in the case of SPC, thesimulation does predict well the experimental peak heights andthe peak positions. A better agreement is observed at lower q

151, 065104-5

The Journalof Chemical Physics

ARTICLEvalues for all the water models. Note that the first lobe of the X-ray
scattering profile represents the structural properties of the bilayer
and therefore better agreement of the first simulated lobe with theexperimental one indicates the better mimicking of the lipid bilayer
structure.
Neutron scattering length density (NSLD) profiles of protonand deuterium are sharply different. Thus, we calculate the NSLD
for the lipid bilayer relative to the deuterons of heavy water and thereciprocal space neutron scattering curves give the neutron scattering form factors. As expected, the neutron form factor increases withincreasing concentration of deuterium. The simulated form factors
match well with the experimental values for q < 0.2. We see better
agreement with experimental form factors for TIP3P, TIP3P-FB, and
TIP4P-FB water models. SPC and TIP4P/2005 water models slightly
overestimate the NSLD value at q = 0.
5. Deuterium order parameter of the acyl chain
and the head group
The ordering of the acyl chains is determined using the deuterium order parameter SCD , which can be calculated using thefollowing equation:
1
SCD = ⟨3cos2 θ − 1⟩.
2(6)Here, θ is the angle between a C−−D bond of the methylene groupof the hydrophobic tail and the bilayer normal. The angular bracket
indicates the ensemble average. Since no hydrogen/deuterium atom
is present in the united-atom model of the lipid molecule, one
needs to construct the positions of the methylene bound hydrogen/deuterium using the neighboring carbon atoms’ positions andassuming a tetrahedral geometry around the carbon.28,83
There is another way to obtain SCD , which we have adopted
in our calculation. We have used the “gmx_order” analysis toolscitation.org/journal/jcp(as incorporated in Gromacs 2019), which calculates SCD using thefollowing equation:2
1
SCD = Sxx + Syy ,
(7)
3
3
2
1
where Sxx = 2 ⟨3cos θx − 1⟩ and θx is the angle between the x axis
and the bilayer normal z axis. For this, the order of atoms should becorrectly provided. Note that Eq. (7) assumes the same SCD values forthe two hydrogen atoms of the methylene carbon atoms. To rectify
this issue, the equations are modified as follows:
√
2
1
2 2
SCD1 = Sxx + Syy −
Sxy ,
(8)
3
3
3
√
1
2 2
2
Sxy .
(9)
SCD2 = Sxx + Syy +3
3
3
Equations (8) and (9) work well for saturated carbons, but for theunsaturated site, they are slightly modified as follows:√1
3
3
Sn = Szz + Syy −
Syz ,
(10)
4
4
2
√
1
3
3
Sn+1 = Szz + Syy +
Syz .
(11)
4
4
2
These equations assume an angle of 120○ around the double bond, which has been taken care of in the fixed version of
the “gmx_order” program (“fixed-gmx_order”) provided by Van
Lehn et al.84 However, this modified tool still produces wrong
results for the second carbon of the double bond. Further modification has been done by Piggot et al.85 The SCD values of theunsaturated methylene positions of POPC are recalculated using the“fixed-gmx_order” tool.
The simulated SCD values for the 9th and the 10th methylene
carbon are in better agreement with the experimental values. (Thetool provides the SCD values same for the two hydrogen atoms.)

FIG. 3. Deuterium order parameters forthe sn1 and sn2 chains of the POPC [(a)and (b)] and DPPC [(c) and (d)], simulated with different water models: SPC(red), TIP3P (blue), TIP4P/2005 (pink),
TIP3P-FB (dark red), and TIP4P-FB(dark gray). Black filled squares present
the experimental values. Error bars ofsimulated data are the standard error,
calculated from block averaging. See
Fig. 1 for the labels of the two chains.

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-6

The Journalof Chemical Physics

ARTICLETABLE V. Order parameter SCD for the head group carbon atoms of POPC andDPPC lipid bilayers simulated in different water models and their comparison withexperimental data.

Lipid

Water model

α

β

γ1

γ2

γ3

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

0.08
0.09
0.09
0.09
0.10

0.02
0.03
0.03
0.01
0.02

0.05
0.10
0.13
0.07
0.12

0.09
0.13
0.15
0.11
0.15

0.15
0.21
0.23
0.17
0.23

Expt.

0.04

0.05

0.13

0.2

0.21

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

0.09
0.11
0.10
0.11
0.09

0.004
0.03
0.02
0.01
0.02

0.04
0.09
0.11
0.03
0.09

0.10
0.16
0.16
0.12
0.16

0.17
0.22
0.23
0.18
0.23

Expt.

0.05

0.02

0.07

0.21

0.22

Literature

DPPC

Literaturescitation.org/journal/jcp

However, for the saturated chains of both DPPC and POPC, theresults are found to be invariant for the original “gmx_order” andthe fixed version. Thus, we have provided the results for the sn1
chain of POPC and for both the sn1 and sn2 chains of DPPC as
obtained from the original “gmx_order” tool in Fig. 3. It is worth
mentioning that the original “gmx_order” tool can still provide SCD
values of the unsaturated position, which are in semiquantitative
agreement with the experimental values. The method of getting
semiquantitative results using the original “gmx_order” tool and theperformance analysis of both the original and the fixed “gmx_order”
tools have been detailed in Sec. S3 of the supplementary material.
SCD gives an insight into the overall ordering of the bilayer. For
example, a higher value indicates an increase in the overall orderedarrangement. A lower SCD value is generally observed for unsaturated lipids, such as POPC.86 The presence of the double bond introduces kinks in the lipid molecule and makes the system less ordered.
Figure 3 presents SCD for both the sn1 and sn2 chains of DPPC
and POPC in different water models. These simulated values arecompared with the experimental values, obtained from NMR based
measurements.87–89 For both the DPPC and POPC bilayers, the SPC

FIG. 4. Probability distribution of theangles between the bilayer normal anddifferent residues: the head group vector P-N [(a) and (b)], the C==O vectorof the sn1 chain [(c) and (d)], and the
C==O vector of the sn2 chain [(e) and (f)]for the DPPC and POPC lipid bilayers in
different water models. The distributions
are color-coded for different water models: SPC: red; TIP3P: blue; TIP4P/2005:
pink; TIP3P-FB: dark red; and TIP4P-FB:
dark gray. The experimental tilt value90
for vectors P-N and the simulated values
are shown in Table VI.

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-7

The Journalof Chemical Physics

and TIP4P/2005 water models overestimate the experimental values. TIP3P, TIP3P-FB, and TIP4P-FB reproduce the experimental
order parameter more accurately for both the sn1 and sn2 chainsof DPPC and POPC bilayers. Consistent with the experimental values for the first few carbons in DPPC, the order parameter remains
almost the same. The sharp drop for the double bond of oleoyl chain(sn2) of POPC is more accurately reproduced by TIP3P, TIP3P-FB,
and TIP4P-FB water models.
The ordering of the head group carbon atoms of choline andglycerol groups are calculated in different water models. The methylene carbons in the choline head group are termed α and β. Twomethylene carbons and one methanetriyl carbon in the glycerol
group are termed γ1, γ3, and γ2, respectively. The different groups
are labeled in Fig. 1. The SCD values in different water model show
characteristic differences in their ordering. The results are comparedwith experimental values for assessing the obtained structures in different water models. The SCD values are listed in Table V. Similar
to the previous simulation, we observed that the SCD for α carbon
is always higher than the β carbon for both the DPPC and POPC
bilayers in all the water models. For α carbons, the simulated values are not in good agreement with the experimental values for boththe lipid bilayer. Better agreement is obtained for β carbon atoms
with values in the range of 0.02–0.03, which is less than that ofthe α carbon. For γ1, γ3, and γ2 carbons, TIP3P, TIP3P-FB, and
TIP4P-FB water models better predict the values. The groups areless ordered in SPC and TIP4P/2005 water models. We will see in
the next paragraph that the SCD values for the head group carbons
closely follow the angular distribution of the head groups to the lipid
bilayer normal.
6. The orientation of the head groups
and the carbonyls
The fluctuation of the head group of a lipid provides a qualitative insight into the flexibility of the lipid bilayer. This is analyzed by
the angular distribution of characteristic vector(s) in the head group
region to the bilayer normal. Here, we consider the P-N and C==O

ARTICLE

scitation.org/journal/jcp

TABLE VI. Most probable values of ζ for the head group atom vector P-N, ω1 for
the C==O vector of sn1 chain, and ω2 for the C==O vector of sn2 chain of POPC andDPPC lipid bilayers in different water models.

Lipid

Water model⟨ζPN ⟩ (deg)⟨ω1 ⟩ (deg)⟨ω2 ⟩ (deg)

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

77
75
71
78
72

110
111
110
113
110

135
130
129
132
130

Expt.

70○ –80○

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

79
77
72
80
73

110
109
109
112
109

131
125
126
129
127

Expt.

70○ –80○

Literature

DPPC

Literaturebond vectors of the sn1 and sn2 chains of the two lipids, DPPC and
POPC. See Fig. 1 for the atom labels. Figure 4 presents the following
distributions: (i) probability distribution p(ζ) of angle ζ between the
P-N vector and the bilayer normal, (ii) p(ω1 ) of angle ω1 between the
C==O vector of the sn1 chain and the bilayer normal, and (iii) p(ω2 )of angle ω2 between the C==O vector of the sn2 chain and the bilayer
normal.
For SPC, TIP3P-FB, and TIP4P-FB water models, p(ζ) is centered around 70○ –80○ for both the lipid bilayers. This is in agreement with the experimental value.90 Therefore, on average, thelipid head groups lie parallel with the bilayer surface for both thelipids. The most probable values are given in Table VI. Characteristic differences in the shape of the distribution in different water

FIG. 5. The radial distribution function
g(r) between the oxygen atom of the carbonyl group (OC ) with water oxygen (OW )[(a) and (b)], and between the phosphorus atom of the PO43- group (P) andthe water oxygen (OW ) [(c) and (d)] for
the DPPC and the POPC bilayers, simulated with different water models: SPC(red), TIP3P (blue), TIP4P/2005 (pink),
TIP3P-FB (dark red), and TIP4P-FB(dark gray). See Fig. 1 for the atomlabels.

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-8

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

TABLE VII. Characteristics of the lipid hydration calculated as a water coordination number to a phosphorus (P) and carbonyl
oxygen (OC ) atom in its first solvation shell along with the solvation shell width (FSW) in different water models for POPC andDPPC lipid bilayers. The values in the parentheses are the standard error calculated by block averaging.

CN

FSW (nm)

SystemWater model

P−−OW

OC −−OW

P−−OW

OC −−OW

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

4.36(0.02)
5.86(0.01)
5.48(0.01)
4.35(0.03)
5.65(0.02)

1.18(0.008)
1.44(0.005)
1.43(0.01)
1.12(0.009)
1.42(0.02)

0.46(0.0002)
0.44(0.006)
0.43(0)0.45(0.0007)
0.43(0)0.34(0.0005)
0.33(0.0004)
0.32(0.0005)
0.33(0.0005)
0.32(0.0004)

DPPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

4.41(0.03)
5.95(0.02)
5.50(0.01)
4.39(0.01)
5.62(0.02)

1.10(0.009)
1.37(0.004)
1.40(0.003)
1.05(0.004)
1.42(0.003)

0.47(0.0008)
0.44(0.0004)
0.43(0.0006)
0.46(0)0.43(0.0004)

0.34(0.0004)
0.33(0.0004)
0.33(0.0005)
0.33(0.0008)
0.33(0.0004)

models are observed. For TIP3P, TIP3P-FB, and TIP4P-FB, sharp
peaks are obtained for the Gaussian type distributions. For SPC
and TIP4P/2005, the distributions are broader with the indistinct
peak position. In particular, for SPC, the distribution is unsymmetrical with a hump at ∼110○ , which is 20○ below the bilayer
plane. For TIP4P/2005 water model, the distribution is bimodal
with equal populations at ∼70○ (20○ above the bilayer plane) and∼110○ (20○ underneath the bilayer plane). Although the bimodal
distributions may emanate from a fluid lamellar phase,91 a previous study92 has observed that molecules with preferential headgroup orientations of ±20○ to the membrane plane are quite stable. The occurrence of the additional shoulder/peak at 110○ may
be linked to the dielectric constant of the water model. Both SPC
and TIP4P/2005 water models have the dielectric constant less than
65 (see Table S1 of the supplementary material), while the dielectric constant for other water models is more than 77. The connection between the dielectric constant of the water solvent and thehead group orientation can be understood in terms of the different
ion-dipole or dipole-dipole interaction between the head group andwater.
The above results suggest that the head group atoms of boththe lipids are more ordered in the presence of any of the three watermodels TIP3P, TIP3P-FB, and TIP4P-FB. We note that the SCD

values for the head group carbons, presented in Table V, strongly
follow the angular distribution trend. Higher SCD values of the headgroup come from the wide angular distribution and vice versa.
Both the distributions, p(ω1 ) and p(ω2 ), are found to be almost
insensitive to the water model. This is expected since the C==O
groups are more buried inside the bilayer than the head group. While
the angular distribution p(ω1 ) is centered at ω1 = 110○ , the distribution p(ω2 ) shows the peak at ω1 = 135○ for all the water models inboth the DPPC and POPC lipid membrane.
7. Radial distribution function
In order to understand the hydration structure of the lipid headgroup, we calculate the radial distribution functions g(r) between thelipid head-group atoms and water. Figure 5 presents g(r) between
the P/OC atom of the lipid (see Fig. 1 for the atom labels) and water
OW for different models of water. Interestingly, the hydration structures of POPC and DPPC lipid head groups are very similar to eachother, suggesting a negligible influence of slight modifications in
the hydrophobic chain(s) of the lipid molecule on the hydrationof the lipid head group. We see that the force field of water influences both the position of the first peak and the first solvation shellwidth [FSW, the position of the first minima of g(r)]. For TIP3P,

FIG. 6. Lateral mean square displacement (MSD) of the DPPC (a) and POPC(b) molecules along the xy plane (perpendicular to the lipid bilayer normal) in
different water models SPC (red), TIP3P(blue), TIP3P-FB (dark red), TIP4P/2005(pink), and TIP4P-FB(dark gray).

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-9

The Journalof Chemical Physics

TIP3P-FB, and TIP4P-FB water models, the intensities of the firstpeaks of g(r) between OC −−OW and P-OW for both the lipids arehigher than those for SPC and TIP4P/2005 water models. The larger
first hydration shell width for both the SPC and the TIP4P/2005
water models suggest significantly low penetration of these watermodels in lipid bilayers. The mass density distributions of the phosphate group and water, plotted in Fig. S5 of the supplementarymaterial, clearly show that the SPC and TIP4P/2005 water models penetrate the lipid bilayer much less than the remaining threewater models. This is probably the reason behind the poorer hydration of the lipid head group by the SPC and TIP4P/2005 watermodels.
Now, we calculate the first-hydration-shell coordination number around the phosphate group by integrating the g(r) up to the
distance where g(r) has a minimum. The hydration numbers of thephosphate group and the OC atom, and their first-hydration-shellwidths are listed in Table VII. The first hydration shell coordination
numbers of the phosphate group are less for SPC and TIP4P/2005
water models than those for the remaining three water models byapproximately 1–1.5 water molecules. On the other hand, the firsthydration shell coordination numbers of the OC atom are less for
SPC and TIP4P/2005 water than those for the remaining three watermodels by approximately 0.2–0.3 water molecule. Therefore, we see
comparatively weak penetration of the SPC and TIP4P/2005 waterin the lipid bilayers. This is consistent with the earlier study,17 where
similar low hydration of the phosphate group was observed for the
GROMOS united atom model of lipid simulated with the SPC force
field of water. The relatively higher first solvation shell coordination

ARTICLE

scitation.org/journal/jcp

TABLE VIII. The simulated Dxy values of the two lipid bilayers—each of which is simulated with different water models—and their comparison with earlier simulated andexperimental data. The values in the parentheses are the standard error calculatedby block averaging.

SystemWater model

Dxy /10−8(cm2 /s)

POPC

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

13.07(0.4)
14.1(0.9)9.39(0.5)
4.77(0.4)
9.81(0.5)

Expt.:95,96 13.5 × 10−8
Simulation:86 9 ± 3 × 10−8

SPC
TIP3P
TIP3P-FB
TIP4P/2005
TIP4P-FB

17.21(0.8)
22.97(1.0)
15.24(0.6)
8.34(1.6)
16.86(0.9)

Expt.:95,96 14 × 10−8
Simulation:86 12 ± 3 × 10−8

DPPC

Literature value (cm2 /s)number (∼5.5–6.0)17 around the phosphate groups in all-atom models of the lipid molecule is well reproduced by the GROMOS united
atom model of lipid molecules when more advanced water models(TIP3P or TIP3P-FB or TIP4P-FB) are used. Therefore, the current
results suggested that these water models describe the hydration of
the lipid head group better than the SPC or TIP4P/2005 water modeldo.

FIG. 7. First-order reorientation autocorrelation functions for the head groups (PN) [(a) and (b)], sn1 chain [(c) and (d)],and sn2 chain [(e) and (f)] of the two
lipids, POPC and DPPC, respectively,
simulated with different water models:
SPC (red), TIP3P (blue), TIP3P-FB (darkred), TIP4P/2005 (pink), and TIP4P-FB(dark gray). See Fig. 1 for the atom labels
and the label of two chains.

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-10

The Journalof Chemical Physics

ARTICLEIn addition, the first-hydration-shell widths (FSWs) around thephosphate group for SPC and TIP4P/2005 water models are slightly
higher than for the rest of the water models. The predictions of FSW
with the remaining three water models are closer to those for allatom models of lipids.17 Therefore, while the SPC and TIP4P/2005
water models result in the reasonable structure of lipid membrane,
the TIP3P or TIP3P-FB or TIP4P-FB water models increase theaccuracy of the lipid bilayer structure, which is even comparable to
the results obtained from an all-atom model of lipid.
B. Dynamical properties
1. Lateral diffusion coefficient
DPPC and POPC have a distinct structure in their hydrophobic region and thus have different packing and fluidity. It is known

scitation.org/journal/jcp

that the introduction of unsaturation produces characteristic kinks
and thereby makes less ordered packing. To have an understandingof the fluidity of the lipids, we determine the lateral diffusion coefficient of lipids in the membrane. The lateral diffusion coefficient oflipids in the xy direction Dxy is obtained by fitting the mean squaredisplacement (MSD) using the following equation:lim ⟨∣r(t) − r(0)∣2 ⟩ = 4Dxy t.

t→∞

(12)

Here, Dxy is the lateral self-diffusion coefficient, r(t) and r(0)are the center of mass positions of the lipid molecules at time t and
t = 0, respectively. The single brackets indicate the ensemble average
which is an average over all lipid molecules and many time origins(t = 0). Figure 6 plots the MSD profiles of the POPC and DPPC
lipid molecules in different water models as functions of time. We

TABLE IX. The fitting parameters and the average rotational time for the head group, sn1 and sn2 regions of the POPC andDPPC lipid bilayers in different water models.

Lipid

Water model

Regiona1τ 1 (ps)a2

τ 2 (ns)a3

τ 3 (ns)⟨τ⟩ (ns)

SPC

Head
sn1
sn2

0.37
0.08
0.06

90.12
63.87
40.69

0.39
0.12
0.06

3.27
4.06
2.99

0.24
0.8
0.88

35.1
258
430

9.73
207
379

TIP3P

Head
sn1
sn2

0.36
0.11
0.08

45.30
101.8
76.32

0.28
0.23
0.12

0.54
3.91
3.95

0.36
0.66
0.8

7.20
180
242

2.76
120
194

TIP3P-FB

Head
sn1
sn2

0.32
0.09
0.06

66.26
100.51
74.25

0.31
0.21
0.11

0.77
4.95
4.75

0.37
0.7
0.83

10.8
198
304

4.26
140
253

TIP4P/2005

Head
sn1
sn2

0.33
0.07
0.05

113.18
60.35
38.60

0.46
0.11
0.06

3.46
4.50
3.74

0.21
0.82
0.89

40.7
283
538

10.2
233
479

TIP4P-FB

Head
sn1
sn2

0.31
0.09
0.07

58.92
107.87
75.77

0.31
0.22
0.11

0.63
5.04
4.61

0.38
0.69
0.82

8.04
178
288

3.27
124
237

SPC

Head
sn1
sn2

0.63
0.10
0.07

42.91
72.24
45.34

0.28
0.19
0.09

0.41
4.30
3.99

0.09
0.71
0.84

36.22
275
248

3.40
196
209

TIP3P

Head
sn1
sn2

0.43
0.17
0.12

18.22
120.49114.07

0.49
0.34
0.23

0.11
3.74
5.64

0.08
0.49
0.65

3.07
201
144

0.31
100
94.9

TIP3P-FB

Head
sn1
sn2

0.44
0.13
0.09

25.73
114.19
92.26

0.49
0.32
0.18

0.15
4.11
4.96

0.07
0.55
0.73

5.04
167
154

0.44
93.2
113

TIP4P/2005

Head
sn1
sn2

0.61
0.09
0.07

46.44
68.84
43.80

0.30
0.17
0.08

0.42
4.66
4.51

0.09
0.74
0.85

47.9
307
362

4.47
228
308

TIP4P-FB

Head
sn1
sn2

0.42
0.14
0.09

23.44
119.96
97.92

0.51
0.31
0.19

0.13
3.96
5.13

0.07
0.55
0.72

3.89
169
160

0.35
94.2
116

POPC

DPPC

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

151, 065104-11

The Journalof Chemical Physics

ARTICLE

see three distinct regions in the MSD profiles. In between, the ballistic regime at the short time and the diffusive regime at the longtime, there exists a subdiffusive regime. This reminds of the MSD
profiles obtained for the diffusion of supercooled water.93,94 The Dxyvalues are obtained by fitting Eq. (12) only in the diffusive regime(>1 ns). As the lipid molecules encounter large friction while diffusing, the uncertainty in the MSD analysis is expected to be very high.
Therefore, we may obtain less accurate Dxy values. To encounter
this problem, we analyze two independent trajectories (each of 200
ns duration) for calculating the diffusion coefficients. Table VIII
presents the Dxy values (averaged over two independent trajectories)
for the lipid molecules in different water models. The values are in
reasonable agreement with experimental95,96 and the earlier simulation data.86 Very similar diffusivities are observed for the lipids
in TIP3P or SPC water models. In addition, the diffusion coefficients of the lipids in these water models are higher than those in
the rest models. TIP4P-FB, TIP3P-FB, and TIP4P/2005 water models give almost similar diffusivities of the lipids. Therefore, we canarrange the water models according to the decreasing order of diffusion coefficients of the lipids as follows: TIP3P ∼ SPC > TIP4P-FB∼ TIP3P-FB ∼ TIP4P/2005. Interestingly, the above trend closely follows the following viscosity trend of the water models at the bulk
regime. ηTIP3P ∼ ηSPC < ηTIP4P-FB ∼ ηTIP3P-FB ∼ ηTIP4P/2005 (see Table S1of the supplementary material).
It is worthwhile mentioning that there is a large range in experimental and MD simulation predictions of Dxy values of the lipid. Dxyof a lipid molecule can be estimated experimentally using a numberof methods, including fluorescence techniques, electron paramagnetic resonance (EPR) spectroscopy, quasielastic neutron scattering(QENS), and NMR spectroscopy. A list of values of Dxy for DPPC
lipid, obtained from different experiment and MD simulations, isavailable elsewhere,54 where it is seen that the Dxy values span overthree orders of magnitude. This large range of the experimental Dxyvalues makes it hard to judge a force field of a lipid in predicting thelateral diffusion coefficient.
2. Reorientation autocorrelation function
Reorientation time autocorrelation function C(t) is calculated
for estimating the rotational time of a molecule. C(t)s of lipid
molecules in the bilayer give insight into the fluidity of the lipid
bilayer. Here, we calculate C(t) of 1st order for the head group andthe two chains (sn1 and sn2) of the lipid. For studying the rotationof the head group, we consider the N-CN1 vector. The vectors C1AC1B and C2A-C2B are considered for the rotation of the sn1 and thesn2 chains, respectively (see Fig. 1 for the atom labels).
Figure 7 presents the above three C(t)s as functions of time
for different systems. Evidently, the correlation functions (particularly for the lipid tails) do not decay to zero over an available rangeof time. This makes it challenging to obtain an average rotationaltime by a simple integration of the correlation functions. Thus, we
fit the curves by multiexponential equations. The correlation decays
are nicely fit with tri-exponential functions (biexponential fitting is
not proper),C(t) = a1 e(−t/ τ1 ) + a2 e(−t/ τ2 ) + a3 e(−t/ τ3 ) ; a1 + a2 + a3 = 1.

(13)Here, a1 , a2 , and a3 are the coefficients and τ 1 , τ 2 , and τ 3 are therelaxation times. Therefore, the average rotational time is given by

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishingscitation.org/journal/jcp⟨τ⟩ = a1 τ1 + a2 τ2 + a3 τ3 .
The fitting parameters and the average rotational times fordifferent systems are listed in Table IX.
The average rotational time for the head group is the lowest. The sn1 chains rotate a bit faster than the sn2 chain (particularly for POPC). Other observations are summarized as follows: (i)The rotational time for the head group and the sn2 chain is much
faster for DPPC compared to those for POPC. (ii) The sn1 chain
rotates at comparable time scale for both the lipid in the same watermodel. Of course, the values are not the same, but the difference isthe lowest for the sn1 chain. (iii) The average rotational times for
all regions follow the following trend for different water models:
TIP3P < TIP4P-FB ∼ TIP4P-FB < SPC < TIP4P/2005.
IV. CONCLUDING REMARKSIn this work, we have checked the compatibility of five different water models with the GROMOS96 53A6L united atom modelof two different lipids, DPPC and POPC. A systematic comparisonof the structure and dynamics of the lipid bilayers, simulated withdifferent water models, has been done. We have calculated properties, such as area per lipid, isothermal compressibility modulus,
average volume per lipid, electron density profile, X-ray and neutron
scattering form factors, and bilayer thickness, and deuterium orderparameters of the hydrocarbon chains and the head groups, angular distribution, and radial distribution function for understanding
the structure of the lipid bilayer. Lateral diffusion and reorientation correlation function of the lipid molecules are calculated for
elucidating the dynamics of the lipid bilayer. The broad comparison of these parameters for different water models reveals thatthe simulations better predict the structure and dynamics of thelipid bilayers when either TIP3P-FB or TIP4P-FB water model isused.
Simulations with SPC and TIP4P/2005 water models predict
more densely packed structures as understood from the lower area
per lipid and slower rotation of the lipid molecules. Importantly,
the simulated deuterium order parameters and the scattering form
factors of the lipids are in quantitative agreement with the experimentally measured values when any of the three water models,
TIP3P, TIP3P-FB, and TIP4P-FB, is considered in the simulation.
The worst agreements are found for the SPC and TIP4P/2005 watermodels. This is a crucial observation since the above two quantities
are measured directly from experiments. The rotational time of thelipid molecules in different water models closely follows the area perlipid for the two bilayers in different water models. For example, thelipid bilayers, simulated with SPC or TIP4P/2005 water, are of minimum area, and therefore, the rotation of the lipids is slower in thesewater models. The lateral diffusion coefficient of the lipid molecules,
on the other hand, strongly follows the simulated viscosities of different models of water. Since both the SPC and TIP3P water models
have a lower viscosity coefficient (η ≤ 0.4) than the other models, thelateral diffusion coefficients of both the lipid molecules are higher.
These results indicate that the bulk properties of different models
of water may strongly influence the structure and dynamics of thelipid bilayer. Therefore, it is essential to work with a better modelof water (such as TIP3P-FB and TIP4P-FB), which better predicts
various experimentally measured quantities.

151, 065104-12

The Journalof Chemical Physics

Finally, we can conclude that water models, such as TIP3PFB, TIP4P-FB, and TIP3P, are also very much compatible with the
GROMOS96 53A6L united-atom model of lipid for the simulationof lipid bilayers. More importantly, these water models (particularly TIP3P-FB and TIP4P-FB) overall improve the agreement of
the simulated quantities with the experimental ones. However, theperformance of the TIP4P/2005 water model for quantitative predictions of lipid bilayer properties is not satisfactory, particularly
for properties like deuterium order parameter. These results, therefore, actively support using TIP3P-FB and TIP4P-FB water models(in combination with the GROMOS96 53A6L united-atom model oflipid) for quantitative studies of the structure and dynamics of thelipid bilayer and interfacial water. In addition, our study does notparticularly recommend the use of TIP4P/2005 model of water in
combination with the GROMOS96 53A6L force field of lipid for thesimulation of the lipid bilayer in water.
SUPPLEMENTARY MATERIAL
The supplementary material contains a table listing the physical properties of different water models, the time dependences of theaverage area per lipid of the two lipid bilayers simulated with fivedifferent water models, decomposed electron density profiles for thedifferent groups of the lipid DPPC and POPC, performance comparison of “gmx_order” and “fixed-gmx_order” for the predictionof the SCD values, and the mass density profiles for the head group of
the lipid and water.
ACKNOWLEDGMENTS
We acknowledge IIT Patna for the computational facility. We
acknowledge Professor Norbert Kučerka for providing the experimental X-ray and neutron scattering form factor data for thetwo lipid bilayers. Shakkira acknowledges DST Inspire Fellowship(INSPIRE No. IF180090) for the research fellowships.
REFERENCES
1

A. D. MacKerell, J. Wiorkiewicz-Kuczera, and M. Karplus, J. Am. Chem. Soc.
117, 11946 (1995).
2
A. D. MacKerell et al., J. Phys. Chem. B 102, 3586 (1998).
3
S. E. Feller et al., Biophys. J. 73, 2269 (1997).
4
J. B. Klauda et al., J. Phys. Chem. B 114, 7830 (2010).
5
J. B. Klauda et al., J. Phys. Chem. B 109, 5300 (2005).
6
K. V. Damodaran, K. M. Merz, and B. P. Gaber, Biochemistry 31, 7656 (1992).
7
U. Essmann, L. Perera, and M. L. Berkowitz, Langmuir 11, 4519 (1995).
8
P. B. Moore, C. F. Lopez, and M. L. Klein, Biophys. J. 81, 2484 (2001).
9
J. Wang et al., J. Comput. Chem. 25, 1157 (2004).
10
B. Jojart and T. A. Martinek, J. Comput. Chem. 28, 2051 (2007).
11
C. F. Lopez et al., J. Phys. Chem. B 108, 6603 (2004).
12
W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J. Am. Chem. Soc. 118,
11225 (1996).
13
A. Maciejewski et al., J. Phys. Chem. B 118, 4571 (2014).
14
W. Kulig, M. Pasenkiewicz-Gierula, and T. Rog, Chem. Phys. Lipids 195, 12(2016).
15
J. P. M. Jämbeck and A. P. Lyubartsev, J. Phys. Chem. B 116, 3164 (2012).
16
J. P. M. Jambeck and A. P. Lyubartsev, J. Chem. Theory Comput. 8, 2938 (2012).
17
K. Pluhackova et al., J. Phys. Chem. B 120, 3888 (2016).
18
S. J. Marrink, A. H. de Vries, and A. E. Mark, J. Phys. Chem. B 108, 750 (2004).
19
S. J. Marrink et al., J. Phys. Chem. B 119, 6184 (2019).

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

ARTICLE

scitation.org/journal/jcp

20

S. J. Marrink et al., Chem. Rev. 119, 6184 (2019).
J. Hermans et al., Biopolymers 23, 1513 (1984).
22
L. D. Schuler, X. Daura, and W. F. Van Gunsteren, J. Comput. Chem. 22, 1205
(2001).
23
I. Chandrasekhar et al., Eur. Biophys. J. 32, 67 (2003).
24
O. Berger, O. Edholm, and F. Jähnig, Biophys. J. 72, 2002 (1997).
25
W. V. Gunsteren and H. Berendsen, Netherlands: Biomos (1987), http://www.
gromos.net/gromos87/GROMOS87_manual.pdf.
26
X. Daura, A. E. Mark, and W. F. Van Gunsteren, J. Comput. Chem. 19, 535(1998).
27
C. Oostenbrink et al., J. Comput. Chem. 25, 1656 (2004).
28
D. Poger and A. E. Mark, J. Chem. Theory Comput. 6, 325 (2010).
29
D. J. Price and C. L. Brooks III, J. Chem. Phys. 121, 10096 (2004).
30
I. G. Tironi, R. M. Brunne, and W. F. van Gunsteren, Chem. Phys. Lett. 250, 19(1996).
31
C. J. Högberg, A. M. Nikitin, and A. P. Lyubartsev, J. Comput. Chem. 29, 2359
(2008).
32
F. Sajadi and C. N. Rowley, PeerJ 6, e5472 (2018).
33
A. Srivastava, S. Malik, and A. Debnath, Chem. Phys. 525, 110396 (2019).
34
L.-P. Wang, T. J. Martinez, and V. S. Pande, J. Phys. Chem. Lett. 5, 1885 (2014).
35
J. L. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).
36
C. L. Moyer and R. Y. Morita, Encyclopedia of Life Sciences (Wiley, 2007),Vol. 1.
37
D. van der Spoel, P. J. van Maaren, and H. J. Berendsen, J. Chem. Phys. 108,10220 (1998).
38
A. Glättli, X. Daura, and W. F. van Gunsteren, J. Chem. Phys. 116, 9811 (2002).
39
H. S. Ashbaugh et al., J. Chem. Phys. 132, 124504 (2010).
40
C. Vega, E. Sanz, and J. Abascal, J. Chem. Phys. 122, 114507 (2005).
41
C. Vega and J. Abascal, J. Chem. Phys. 123, 144504 (2005).
42
L. Martínez et al., J. Comput. Chem. 30, 2157 (2009).
43
N. Kučerka, S. Tristram-Nagle, and J. F. Nagle, The J. Membr. Biol. 208, 193(2006).
44
N. Kucerka et al., Biophys. J. 95, 2356 (2008).
45
D. Van Der Spoel et al., J. Comput. Chem. 26, 1701 (2005).
46
B. Hess et al., J. Comput. Chem. 18, 1463 (1997).
47
R. N. Lewis, B. D. Sykes, and R. N. McElhaney, Biochemistry 27, 880 (1988).
48
S. Mabrey and J. M. Sturtevant, Proc. Natl. Acad. Sci. U. S. A. 73, 3862 (1976).
49
S. Nosé, J. Chem. Phys. 81, 511 (1984).
50
W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
51
M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
52
Y. Liu and J. F. Nagle, Phys. Rev. E 69, 040901 (2004).
53
J. B. Fournier, A. Ajdari, and L. Peliti, Phys. Rev. Lett. 86, 4970 (2001).
54
D. Poger, B. Caron, and A. E. Mark, Biochim. Biophys. Acta, Biomembr. 1858,
1556 (2016).
55
J. M. Smaby et al., Biophys. J. 73, 1492 (1997).
56
G. Büldt et al., J. Mol. Biol. 134, 673 (1979).
57
R. L. Thurmond, S. W. Dodd, and M. F. Brown, Biophys. J. 59, 108 (1991).
58
H. Heller, M. Schaefer, and K. Schulten, J. Phys. Chem. 97, 8343 (1993).
59
S.-W. Chiu et al., J. Phys. Chem. B 113, 2748 (2009).
60
H. Binder and K. Gawrisch, J. Phys. Chem. B 105, 12378 (2001).
61
J. F. Nagle and S. Tristram-Nagle, Biochim. Biophys. Acta, Biomembr. 1469,
159 (2000).
62
E. Lindahl and O. Edholm, Biophys. J. 79, 426 (2000).
63
C. Anézo et al., J. Phys. Chem. B 107, 9424 (2003).
64
C.-H. Huang and J. P. Charlton, J. Biol. Chem. 246, 2555 (1971),
http://www.jbc.org/content/246/8/2555.
65
M. C. Wiener et al., Biochim. Biophys. Acta, Biomembr. 938, 135 (1988),
https://doi.org/10.1016/0005-2736(88)90153-8.
66
S. Tristram-Nagle, H. I. Petrache, and J. F. Nagle, Biophys. J. 75, 917 (1998).
67
D. L. Melchior and H. J. Morowitz, Biochemistry 11, 4558 (1972).
68
D. Wilkinson and J. Nagle, Anal. Biochem. 84, 263 (1978).
21151, 065104-13

The Journalof Chemical Physics

69

S. H. White, R. E. Jacobs, and G. I. King, Biophys. J. 52, 663 (1987).
M. Posch et al., J. Biol. Chem. 258, 1761 (1983), http://www.jbc.org/content/
258/3/1761.short.
71
R. Koynova et al., Chem. Phys. Lipids 48, 205 (1988).
72
D. Uhríková et al., Chem. Phys. Lipids 145, 97 (2007).
73
R. S. Armen, O. D. Uitto, and S. E. Feller, Biophys. J. 75, 734 (1998).
74
G. Pabst et al., Phys. Rev. E 62, 4000 (2000).
75
L. J. Lis et al., Biophys. J. 37, 657 (1982).
76
N. Kučerka, J. Katsaras, and J. F. Nagle, J. Membr. Biol. 235, 43 (2010).
77
R. W. Benz et al., Biophys. J. 88, 805 (2005).
78
N. Kučerka et al., Biophys. J. 95, 2792 (2008).
79
M. Roark and S. E. Feller, J. Phys. Chem. B 113, 13229 (2009).
80
N. Kucerka et al., J. Phys. Chem. B 116, 232 (2012).
81
J. P. Jämbeck and A. P. Lyubartsev, J. Chem. Theory Comput. 9, 774 (2012).
82
J. B. Klauda et al., Biophys. J. 90, 2796 (2006).
70

J. Chem. Phys. 151, 065104 (2019); doi: 10.1063/1.5108830
Published under license by AIP Publishing

ARTICLE

scitation.org/journal/jcp

83

P. Niemelä, M. T. Hyvönen, and I. Vattulainen, Biophys. J. 87, 2976 (2004).
R. C. Van Lehn and A. Alexander-Katz, J. Phys. Chem. B 118, 12586 (2014).
85
T. J. Piggot et al., J. Chem. Theory Comput. 13, 5683 (2017).
86
S. Ollila, M. T. Hyvönen, and I. Vattulainen, J. Phys. Chem. B 111, 3139 (2007).
87
T. Huber et al., J. Am. Chem. Soc. 124, 298 (2002).
88
H. I. Petrache, S. W. Dodd, and M. F. Brown, Biophys. J. 79, 3172 (2000).
89
J. Seelig and N. Waespe-Sarcevic, Biochemistry 17, 3310 (1978).
90
J. P. Doux, B. A. Hall, and J. A. Killian, Biophys. J. 103, 1245 (2012).
91
K. Tu et al., Biophys. J. 70, 595 (1996).
92
L. Saiz and M. L. Klein, J. Chem. Phys. 116, 3052 (2002).
93
S. Indra and S. Daschakraborty, Chem. Phys. Lett. 685, 322 (2017).
94
V. Dubey, N. Kumar, and S. Daschakraborty, J. Phys. Chem. B 122, 7569
(2018).
95
A. Filippov, G. Orädd, and G. Lindblom, Biophys. J. 84, 3079 (2003).
96
G. Lindblom, G. Orädd, and A. Filippov, Chem. Phys. Lipids 141, 179 (2006).
84

151, 065104-14