J. Phys. Chem. 1995, 99, 10035-10042
10035
Constant Pressure and Temperature Molecular Dynamics Simulations of Crystals of the Lecithin Fragments: Glycerylphosphorylcholine and Dilauroylglycerol Kechuan Tu, Douglas J. Tobias, and Michael L. Klein* Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323
Received: August 19, 1994; In Final Form: April 7, 199.5@
We have performed constant temperature and pressure classical molecular dynamics simulations of crystals of the lecithin fragments glycerylphosphorylcholine (GPC), a model for the headgroup and glycerol regions, and dilauroylglycerol (DLG), a model for the glyceroYacyl ester and hydrocarbon regions. By comparing the simulation results to experimental X-ray crystal structures, we have evaluated the quality of all-atom force fields for phospholipids. For GPC we examined two existing force fields: the Stouch et al. (J. Comput. Chem. 1991, 12, 1033) and CHARMM 22 (MacKerell, A.; Schlenkrick, M.; Brickmann, J.; and Karplus, M. Unpublished) potentials. For DLG we employed a force field that is a combination of the Stouch et al. parameters for the glyceroVacy1 ester linkage and a compilation of new and existing potentials for the hydrocarbon chains. The average unit cell parameters from the GPC simulation with the Stouch et al. potential were all within 2% of the experimental values, and the density was less than 4% too low. The agreement between the simulation and experimental cell parameters was better, and the density error was 3%, for DLG. The agreement of the CHARMM 22 model for GPC was not as good: the deviation of the unit cell lengths was as much as 7%, and the density was 11% too high. Overall, the molecular structures and hydrogenbonding networks were well-reproduced by the DLG simulation and the GPC simulation employing the Stouch et al. potential. The most significant deviations were in the torsion angles in the acyl ester linkages of DLG, but the deviations in adjacent torsions offset each other, leaving the hydrocarbon chain packing unchanged. In general, the results of our simulations, in which the Ewald method was used to sum the electrostatic interactions, agree even better with the experimental data than the constant energy and volume simulations of Stouch et al.. in which the electrostatic interactions were truncated.
1. Introduction In principle, molecular dynamics (MD) simulations can provide atomic-scale structural and dynamical information for the discussion of the structure and function of biological membranes. Such information has not been accessible experimentally because the disorder in biological and model membranes, especially in the biologically relevant liquid crystal phase, has made quantitative structure determination impossible with present techniques. The lack of high-resolution experimental data on membranes is simultaneously good and bad for modelers: the void leaves a lot of room for predictions from simulations but, at the same time, makes it difficult to assess their credibility. Before carrying out extensive simulations of membrane systems where high-resolution experimental data is just emerging,’-3 a reliable simulation methodology (potential energy functions, algorithms) should be established. Several reports of MD simulations of model biological membranes (lipid monolayer^^,^ and bi1ayers6-l5) have appeared recently in the literature. For the most part, these studies employed potentials that were derived, after small modifications, from protein and nucleic acid simulations. However, there is no guarantee that these potentials will work well in simulations of membranes where the relative importance of different types of interactions (hydrophobichydrophilic, electrostatic, etc.) is likely quite different. In the one case where careful tests of a potential were reported, Stouch and colleagues studied crystals of lecithin fragments.I6.l7 They tested their potential by using constant volume and energy (NVE) MD simulations of flexible molecules16 and performed optimizations of lattice parameters by using energy minimization of rigid molecules in a flexible @
Abstract published in Advance ACS Abstracts, May 15, 1995.
lattice.” One question that remains unanswered from this work is how well does the potential perform when both the molecular and crystal degrees of freedom are allowed to relax at finite temperature? Furthermore, all of the membrane simulations to date have used serious, untested approximations, Le. in the handling of long-ranged electrostatic forces4 and/or boundary conditions’ for computational convenience and/or expediency. It is likely that the results of lipid simulations depend sensitively on how the electrostatic interactions are treated. All but onell of the lipid simulations to date have employed truncation of electrostatic forces. However, there is no well-prescribed procedure for determining the optimal truncation scheme, and it has been shown that poorly chosen schemes can introduce serious artifacts into the simulation result^.^ Finally, most of the membrane simulations have been carried out at constant volume, using experimental values for the system dimensions (surface are/lipid, bilayer thickness, interlamellar spacing), although there is considerable disagreement on the correct values for these parameters even for the most studied lipid.I8 It is difficult to tell whether or not the chosen dimensions are appropriate for the simulated systems because, in most of the published simulation results, the pressure in the system was not reported. In general, where comparison is possible, the structural properties of phospholipid bilayers computed from constantvolume MD simulations, for example, atomic density profiles, bilayer thicknesses, and interlamellar spacings,*.’ and hydrocarbon chain order parameter profiless~lO~’ agree qualitatively with experimental measurements. The surface areaslmolecule and interlamellar repeat distances from a short (80 ps) constant pressure simulation of a hydrated dipalmitoylphosphatidylcholine (DPPC) bilayerI5 fall within the range of experimental
0022-365419512099-10035$09.0010 0 1995 American Chemical Society
’
10036 J. Phys. Chem., Vol. 99, No. 24, I995 estimates, although the slow drift in the volume observed at the end of a 200 ps simulation of a dipalmitoylphosphatidylethanolamine bilayerI4 suggests that much longer simulation times are needed for convergence at constant pressure. Furthermore, short-time dynamical properties, such as diffusion constants for “rattling” motion, exhibit order-of-magnitude agreement with experimental data,lO~’l and acyl chain CH reorientation times agree within 10-30%.10 It is natural to ask how much the results (especially regarding structure) of the constant volume simulations would change if the calculations were carried out at constant pressure and the interactions between the molecules (i.e. the potentials) were allowed to determine the system dimensions. Moreover, would systems containing charged or zwitterionic lipids in which the long-ranged electrostatic forces are truncated remain stable on order-of-magnitude greater time scales than are currently being explored, where interesting dynamical phenomena like headgroup rotation and lipid lateral diffusion can be observed and where acyl chain order parameters are likely to start converging?I9 Recently, a bilayer of an uncharged, nonzwitterionic lipid, glycerol monooleate, was shown to be stable on the nanosecond time scale in a simulation that employed electrostatic truncation,I3 but this is likely a very different situation from the case of a charged or zwitterionic lipid. To further address some of these important issues, we have performed simulations of systems for which detailed structural and, in some cases, dynamical and thermodynamic properties are known experimentally, i.e. fragments of phospholipid molecules: crystalline glycerylphosphorylcholine, a model for the headgroup/glycerol portion; crystalline dilaurylglycerol, the glyceroUacy1 ester part; and solid and liquid alkanes, models for the hydrocarbon tails. These simulations were carried out under the same conditions as the experiments by using a novel constant temperature and pressure (NPT) MD algorithm developed recently in our 1aboratory2Oand an Ewald summation2Iof the long-range forces. The use of fully flexible cells in constant NPT MD simulations provides a stringent test, as the unit cell sizes and shapes are sensitive to the details of the potentials. By combining existing and new force fields, we arrived at a set of potentials that reproduce the experimental densities and crystal structures of the molecules investigated, as well as the energies and diffusion constants in cases where they are known. In this paper we report results for the lecithin fragments, glycerylphosphorylcholine(GPC) and dilauroylglycerol (DLG). The investigations of alkanes will be reported elsewhere.**
2. Methodology 2.1. Potential Energy Function. The potential energy function (PEF), which gives the potential energy as a function of the positions of the atoms in the system, is the primary ingredient used in molecular mechanics and molecular dynamics or Monte Carlo simulation studies. It is the function from which the forces are derived and used, either to direct an energy minimization or to propagate the equations of motion in a dynamics simulation. The development of a PEF takes place in three stages: an empirical (usually analytical) functional form is chosen, parametrized, and checked and, if necessary, refined. In our work, the PEF is of a standard macromolecular mechanics f ~ r m , ’containing ~ . ~ ~ six types of terms: harmonic bond stretching and bond angle bending terms, periodic and harmonic dihedral angle (torsion) potentials, and intra- and intermolecular van der Waals and Coulomb interactions. The parameters of the PEF include the force constants and equilibrium geometrical values in the bond and angle potentials, the van der Waals radii and interaction strengths, and the partial atomic charges in polar moieties.
Tu et al. Here, the PEF is parametrized for all of the atoms in the system. This is in contrast to the model known as the “polarH ’model, commonly employed in simulations of biological molecules,24 including l i p i d ~ , ~ 1,14,15 - ~ $ ’ where the nonpolar hydrogen atoms are not explicitly represented but rather are included implicitly through modification of the heavy atoms to which they are bonded. Studies of systems with densely packed hydrocarbon chains have shown that an all-atom representation is required for a proper description of the chain structure and dynamics.25 We chose to employ an all-atom representation because we are interested in studying various phases of lipids (crystal, gel, and liquid crystal) using the same model. Due to the lack of high-resolution structural data on membranes, the best way to check the potentials is to carry out simulations of lipid crystals. It is desirable to simulate a crystal of an intact lecithin molecule, but since it is difficult to obtain well-ordered crystals of these molecules, no highly resolved crystal structures are available. Although the crystal structure of dimyristoylphosphatidylcholine DMPC has been reported,26 the resolution is too low to be used in a quantitative comparison. However, high-resolution crystal structures are available for GPC27and DLG,28two fragments of a lecithin molecule. GPC represents the polar headgroup and glycerol regions, and DLG represents the glycerollacy1 ester and hydrocarbon regions. Together these two molecules contain all of the functional groups, structural features, and interactions that occur in the lecithin molecules. Thus, as in previous work,16 we used simulationsof the GPC and DLG crystals to check the potentials, and we assume that the results can be extrapolated to the lecithins. At present, there are two all-atom potentials available for lecithin^.'^^^^ The first was developed by Stouch et a1.I6 and has been used recently in simulations of a DMPC monolaye# and bilayer.’* In this model, the parameters for the intemal degrees of freedom were taken from the Discover30 and AMBER3’ protein and nucleic acid force fields. The van der Waals interactions are treated as Lennard-Jones potentials with the Discover parameters. The partial charges were determined by fitting the electrostatic potentials from ab initio calculations on DLG and GPC. Although the model of Stouch et al. has been tested by energy minimization and constant NVE simulations of lipid ~rystals,’~-’’ one of our objectives in this paper is to subject the headgroup portion to a more stringent test, namely, a constant NPT simulation of the GPC crystal with a more rigorous treatment of the long-range forces. We also evaluate the headgroup potential from a similar all-atom lipid model, the CHARMM 22 for comparison with the Stouch model. The CHARMM 22 model has been recently employed in a MD simulation of a DPPC bilayer.’O In recent work, there have been two quite different approaches to the design of all-atom models for hydrocarbons. One a p p r ~ a c h ’ ~is, ~a* fully flexible representation with partial charges on the C and H atoms. The overall torsional potential is made up of contributionsfrom a single cosine “intrinsic” term for each torsion in the molecule and from 1-4 Lennard-Jones and Coulomb interaction^.^^ Another approach33is a partially rigid model in which all of the intemal degrees of freedom involving H atoms (except for terminal methyl group rotations) are frozen by the use of holonomic constraints, enabling a relatively large MD time step to be used. This model has no charges and uses an exponential-6 potential due to Williams34 for the intra- and intermolecular nonbonded interactions and cosine power series torsional potential^.^^ It was verified by a constant pressure simulation of solid octane that showed good agreement with the experimental crystal structure. Given the
Glycerylphosphorylcholine and Dilauroylglycerol TABLE 1: Fitted Torsion Parameter@ dihedral angles! d k (kcavmol) C’CT-0-0’ CT-C’-0-CT
11.60
-0.800 2.150 0.00 0.00 0.10 0.10 0.16 0.581 0.1 10 0.341 0.410 0.096 0.234 0.177 0.165
HA-CT-C’-0 HA -CT -C’ -0’ CT-CT-C’-0’ CT-CT-C’-0 C’-CT-CT-HA C’-CT-CT-CT CT-CT-CT-CT HA-CT-CT-CT HA-CT-CT-HA
J. Phys. Chem., Vol. 99, No. 24, 1995 10037
n
6 (deal
2 1 2 3 3 2 2 3 1
180.0 180.0 180.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 3 3
2 1 3 3
+
The form of the torsion energy is E, k,[1 cos(n,@,- 6J]. Atom types are defined as follows: C’ and 0’are carbonyl carbon and oxygen, 0 is an ether oxygen, and CT and HA are aliphatic carbon and hydrogen atoms.
TABLE 2: Fitted Lennard-Jones Parameter@ atom HA CT
E
(kcalhol) 0.0097 1 0.0946
(4 2.797 3.439
“ For two atoms separated by the distance r, the Lennard-Jones energy is 4c[(u/r)’* - ( u / ~ ) ~ I . proven success of the exponential-6 potential, as well as computational considerations, we recently combined the two approaches to construct another hydrocarbon potential.22 Anticipating the use of multiple time step MD methods,36 we preferred to develop a potentially fully flexible model over the cumbersome constraints. The use of the exponential-6potential, as opposed to the Lennard-Jones plus Coulomb approach, eliminates a large number of long-ranged electrostatic interactions. Accordingly, we took the molecular geometry and nonbonded parameters from Williams34but the force constants from Smith and K a r p l ~ s . ~As * the Smith-Karplus torsional potentials appear to be realistic,37we fit a cosine Fourier series to them for use in our alkane model (see Table 1). We recently demonstrated by using constant NPT MD simulations that this model exhibits good agreement with experimental results for the crystal structures and/or densities and energies of both solid and liquid alkanes over a range of chain lengths.22 For use in the lipid molecules, where most of the van der Waals interactions are modeled by Lennard-Jones potentials, we fit LennardJones potentials to the Williams potentials. We use these Lennard-Jones parameters (Table 2) with the traditional geometric mean combining rules for the interactions of the hydrocarbon chains with the remaining atoms. Finally, in order to “connect” the new hydrocarbon chains to the remainder of the lipid molecules through the acyl ester linkages, we used the bond stretching and angle bending potentials of Stouch et al.I6 However, since the model of Stouch et al. employs 1-4 van der Waals and electrostatic contributions to produce the overall torsional potential and we have avoided such interactions in our hydrocarbon chains, we discarded the 1-4 interactions and fit the original Stouch et al. torsional potentials to a Fourier cosine series. Table 1 lists the resulting torsion parameters for the acyl ester region. 2.2. Constant Pressure and Temperature Molecular Dynamics. We have employed the so-called “hybrid” scheme of Martyna et aL20 for carrying out MD calculations in fully flexible simulation cells under the conditions of constant temperature and pressure. Technical details of the method and
its implementation are given elsewhere,20so we give here only a brief description of the salient features. The “hybrid” scheme is an “extended system” (ES) in which the equations of motion for the positions and momenta are supplemented by equations of motion for thermostat and barostat degrees of freedom (and their momenta). The name “hybrid” comes from the fact that both isotropic and anisotropic volume fluctuations are allowed by the introduction of separate equations of motion for the volume and the cell variables, together with the constraint that the cell variables are consistent with the volume at any instant in time. This ES method is attractive because it produces continuous dynamics with a welldefined conserved quantity and, assuming the dynamics is ergodic, it generates the isothermal-isobaric ensemble and satisfies the appropriate tensorial virial theorem. To help ensure ergodicity, the particles and the barostat variables are coupled to separate NosB-Hoover chain thermostat^.^^ Finally, as the “hybrid” method is modularly in~ariant,~’ the results do not depend on the choice of the basis for the lattice vectors. 2.3. Details of the Simulations. The starting structures for the simulations were the X-ray crystal structures of GPC27and DLG.28 Space group operators and lattice translations were used to generate all of the atomic coordinates from the asymmetric units. The GPC simulation cell contained 96 molecules in a 3a x 4b x 2c lattice, and the DLG cell contained 48 molecules in a 6a x 4b x c lattice. The simulations were carried out at the temperature 300 K and zero pressure. The ES equations of motion were integrated by using an iterative Verlet-like algorithm42with a time step of 1 fs. The SHAKE algorithm43was used to constrain the lengths of bonds involving hydrogen atoms. The fictitious masses of the ES variables were chosen according to the prescription given by Martyna et al.,20with time scales of 0.5 ps for the thermostats and 1 ps for the volume and cell variables. The NosB-Hoover chain length was five. The Ewald method2I was used to compute all of the electrostatic interactions in the crystals. Minimum image periodic boundary conditions2I were employed to calculate the van der Waals interactions and the real-space part of the Ewald sum with simple truncation at 10 A. Long-range corrections to account for the truncated van der Waals interactions2I were included in the energies and pressures. The long-range correction to the energy is computed from the following formula:
where Uv(r) is the long-ranged part of the pair potential, gv(r) is the pair correlation function between sites i andj, r, is the cutoff, and V is the volume. Taking gv(r>rc) = 1 and U , = -Cjjr-6, i.e. the attractive part of the Lennard-Jones or exponential-6 potentials, we obtain
The corresponding correction to the pressure is simply C/v. We used the CHARMM program,23version 23, as modified by us to implement the ES dynamics, minimum image periodic boundary conditions, and Ewald summation for all of the calculations. Each simulation consisted of 30 ps constant NVE and 10 ps constant NPT equilibration, followed by a 40 ps constant NPT data collection period. These modest simulation lengths were sufficient, as there is not much molecular flexibility in the GPC and DLG crystals and, hence, all of the computed
Tu et al.
10038 J. Phys. Chem., Vol. 99, No. 24, 1995
y chain 31.2
91 h
h
", 6
8 kl
31
3
30.8
U
30.6
0
10
20 30 40 time (ps)
90.5 90 89.5
0
50
10
20 30 time (ps)
40
50
102.5
31.6 h
3
E?
31.4
a
d
Y
31.2
~
102
p chain
101.5 101
Figure 3. Atom and torsion angle labeling in a lecithin molecule (phosphatidylcholine), after Hauser et al.44
100.5
31
0
10
20 30 40 time (ps)
0
50
10 20 30 time (ps)
50
40
91
33.4 h
8
33.2
e
", 0
33
90.5
3
90
r-
89.5
v
32.8
89
0
10
20 30 time (ps)
40
LI , , , , I , , , , I , ,, , I , , , , d 0
50
10 20 30 time (ps)
40
50
Figure 1. Time evolution of the simulation cell parameters (3a x 4b x 2c lattice) during the constant pressure simulation of the Stouch model of GPC.
-
91 34
h
8
0
", 1
90
$
33.5
Y
e
89
33 0
10
20 30 40 time (ps)
50
0
FL""I""I""I""I""I~
8 t! ? v
30.4
a
v
30.2
~
30
""
I
"
50
' I' ' 0 ' q
94 93.5 93
0
10
20 30 40 time (ps)
0
50
10
20 30 time (ps) I' " ' I i ' " I
40
10 20 30 time (ps)
40
50
"m
n.
YI
34.8
3
'I
" "
40
94.5
30.6
3
10 20 30 time (ps)
h
90.5
34.6
#
0
v
34.4
r-
90 89.5
34.2
0
10
20 30 time (ps)
40
50
0
50
Figure 2. Time evolution of the simulation cell parameters (6a x 4b x c lattice) during the constant pressure simulation of DLG. properties coverged during the 10 ps constant NPT equilibration (see Figures 1 and 2). 3. Results and Discussion In this section we compare simulation results to experimental X-ray diffraction data for several properties of L-2-glycerylphosphorylcholine (GPC) and 2,3-dilauroyl-~-glycerol (DLG) crystals. The properties of interest include the dimensions of the crystal lattices, the structures and orientations of the molecules,
and the intra- and intermolecular interactions. We will employ the atom and torsion angle nomenclature of ref 44 (see Figure 3). 3.1. L-2-Glycerylphosphorylcholine(GPC). The simulation of GPC was initiated from the crystal structure solved by Abrahamsson and Pascher.*' The GPC molecule includes the glycerol moiety and the choline headgroup that comprise the polar region of a full lecithin molecule. The crystal is monoclinic with four molecules per unit cell and two molecules in the asymmetric unit in the space group P21. The crystal structure can be described as bilayers of one conformer alternating with bilayers of the other which are parallel to the ab plane. Both independent conformers of the asymmetric unit form the same hydrogen-bonding network within their bilayer, but there are no intramolecular hydrogen bonds and no hydrogen bonds between the different conformers. For both independent molecules, the glycerol terminal hydroxyl group is hydrogen bonded to one unesterified oxygen of the phosphate group of a symmetry-related molecule to form the interlayer hydrogen bond. The other glycerol hydroxyl group is hydrogen bonded to the other unesterified oxygen of the phosphate group of a b-translated molecule. This intralayer hydrogen bond links the molecules to produce spirals parallel to the b axis, The geometry of the lattice reflects the balance of the intermolecular forces and the intramolecular geometry and, in a constant NPT simulation, is very sensitive to the details of the force field. Table 3 compares the X-ray data with the average lattice parameters and density from the GPC simulation using the Stouch model. The a and b unit cell lengths are both about 2% too long in the simulation, and c is about 1% too short. The simulation maintained perfectly the monoclinic symmetry of the crystal, although the /3 angle is almost 2% too small and the density is 3.8% too low. Overall, the Stouch model appears to reproduce very well the lattice of the GPC crystal. The data in Table 3 show that, with the exception of the b axis, the lattice parameters and density of the GPC crystal are not as well reproduced by the CHAFWM 22 model. Perhaps the most notable discrepancies in the CHARMM 22 results are that the c axis is about 7% too short, there is a slight distortion from monoclinic symmetry, and the density is 11% too high. It is conceivable that the performance of the CHARMM 22 potential is sensitive to the method used to calculate the electrostatic interactions. We have used the Ewald method in our calculations while the CHAFWM 22 parameters were optimized for use with a 12 8, spherical truncation.29 Regardless of the reason for the relatively poor performance of the CHARMM 22 model observed here, we conclude from the comparison that, for our purposes, the Stouch model is superior for GPC (and lecithin headgroups). Henceforth, therefore, we
J. Phys. Chem., Vol. 99, No. 24, 1995 10039
Glycerylphosphorylcholine and Dilauroylglycerol TABLE 3: Lattice Parameters and Density of the GPC Crystal" a
(A)
b (A) 7.71 7.84(3) 7.65(2)
10.10 10.33(3) 9.69(2)
expt MD Stouchb MD C22'
c (A) 16.62 16.46(5) 15.43(5)
a (deg)
P (de&
Y (deg)
d (g/cm3)
90.0 90.0(4) 90.1(4)
102.7 101.1(3) 99.7(4)
90.0 90.0(4) 90.4(3)
1.32 1.27 1.47
The MD values are averages over 40 ps. The MD root-mean-squared fluctuations are given as the last significant digits in parentheses. This work, Stouch et al. model.I6 This work, CHARMM 22
TABLE 4: Hydrogen Bond Lengths in the GPC Crystalu mol. B
mol. A interlayer intralayer
expt
MD
expt
2.70 2.77
2.77(3) 2.69(3)
2.70 2.69
MD
2.7l(2) 2.71(2)
( I Units are A. The MD root-mean-squared fluctuations are given as the last significant digits in parentheses.
TABLE 5: Root-Mean-Squared Deviation from the Crystal Structure of the Internal Geometric Parameters of GPC bonds MD (this work) MD16
mol. A mol. B mol. A mol. B
expt 30
(A)
angles (den) dihedrals (deg)
0.0282 0.0182 0.0181 0.0283 0.038
2.10 2.64 3.00 2.91 2.01
3.86 4.60 8.01 4.90
will restrict our attention to the results of the Stouch model GPC simulation in the remainder of this section. Hydrogen bonding reflects the short-ranged intermolecular interactions and is likely to play an important role in determining the structure and properties of the membrane surface. As described above, there are two kinds of hydrogen bonds in the GPC crystal: the interlayer hydrogen bond connecting the OH3 atom in each molecule to the 0 1 4 of a symmetry replicate and the intralayer hydrogen bond between the OH2 atom in each molecule and the 0 1 3 atom of a b-translated molecule. Table 4 shows the comparison of the average hydrogen bond distances from the Stouch model simulation to the experimental values. The hydrogen bonds involving molecule B are quantitatively reproduced, while those involving molecule A are slightly off the interlayer hydrogen bond is 0.07 8, too long, and the intralayer bond is 0.08 8, too short. Nonetheless, overall, the hydrogen bonding network is well-preserved by the simulation. Now we turn to an analysis of the internal structure of the GPC molecules. In Table 5 we list the root-mean-squared (rms) deviations, between the 40 ps time-averaged structure and the X-ray structure, of the internal geometric parameters (bond lengths, bond angles, and dihedral angles) involving heavy (nonhydrogen) atoms. We also include for comparison in the table the experimental error (30 for 95% confidence) and the results from the constant NVE simulation of Stouch et a1.l6 The average bond deviation from our constant NFT simulation is well within, and the angle deviation is slightly greater than, the experimental error. The largest deviations, about 4.5",occur in the P-012-Cll and C11-C12-N angles. The dihedral angles show small deviations (about 4-5" on average) from the experimental structure in our simulations, the largest being in the torsion about the P-012 bond. Excluding the bond TABLE 6: GPC Headgroup Torsion Angles" a1 expt A expt B MD A (this work) MD B (this work) MD16 a
165 - I72 166(15) -174(9) -174(9)
a2 -7 1 64 -73( 12) 63(9) W8)
lengths in molecule A, the deviations from the X-ray structure are smaller in our simulation compared to the simulation of Stouch et al. The rms deviation of the average heavy atom positions with respect to the crystal structure was 0.38 8, in our constant NPT simulation. The headgroups (a-chains) of the two molecules in the asymmetric unit exhibit a mirror image symmetry; i.e., their torsion angles have almost the same value but opposite signs. These preferred conformations are similar to those of the parent compound, phosphorylcholine. Restrictions in the headgroup conformation are imposed by the tendency of the zwitterionic headgroups to optimize their intra- and intermolecular interactions and by steric hindrance, which affects the conformations about the P-0 ester bonds. Thus, a detailed analysis of the a-chain torsion angles and the distances between the negatively charged nitrogen atom and the phosphate moiety provides us with a sensitive test of the precision of the intramolecular and short-ranged intermolecular interactions in the force field. Table 6 compares the average a-chain torsion angles from the simulations with the crystal structure data. All of the experimental torsion angles are within the simulation fluctuations except the 011-P-012-Cll torsion angle (a3). The ranges of the torsional fluctuations during our simulation are reflected in the rms values given in the table, and they confirm that no conformational transitions occurred during the MD run. This is in contrast to the simulation of Stouch et al., where the ranges of fluctuations were much larger and suggest that conformational transitions have occurred (see Table IX of ref 16). In Table 7 we compare some intramolecular distances and orientations of vectors in the headgroup with respect to the bialyer (ab) plane Po2 is the distance between the center of gravity of the two charge-sharing unesterified phosphate oxygens and the nitrogen. The distances between the phosphate moiety and the nitrogen atoms are a little large in the simulation, by 0.07-0.15 8, for molecule A and 0.06-0.24 8, for molecule B. This could be due to too strong an 0 - N repulsive interaction or too weak a P-N attractive force in the potential. The inclination of the P-N dipole and the vectors, drawn from the distance of the center of gravity of the two unesterified phosphate oxygens to the nitrogen atom, indicates that the charges are distributed almost parallel to the bilayer surface and, hence, have strong electrostatic interactions with the parallel charged group in the adjacent bilayer. In Table 7 we can see that the orientations of the charged headgroups in the simulation are very close to those observed experimentally. Table 8 compares the simulation results for the average glycerol conformation with the experimental structure. The a3 -59 65 -51( 10) 53(9) 52(8)
a4 -138 140 - 135(9) 142(9) 152(11)
Units are degrees. The root-mean-squared fluctuations from the MD simulations are given in parentheses.
a5 73 -75 68@) -69(9) -73(8)
a6 174 -176 170(7) - 173(8) -176(7)
10040 J. Phys. Chem., Vol. 99, No. 24, I995
Tu et al.
TABLE 7: GPC HeaderouD Structure and Orientationa inclination toward layer plane (deg) C1-C2 bond
P-N dipole
,P02-N dipole
TABLE 10: Root-Mean-Squared Deviation from the Crystal Structure of the Internal Geometric Parameters of DLG
P(:il
bond
MD (this work) MDi6 expt 3a intramolecular distances (A) O( 12)- -N
P- -N
MD A
3.20(1)
4'26(1)
expt A
3.13
4.17
MD B
3.20(1)
4'36(1)
expt B
3.14
4.22
Ejl:$N
P&- -N
5'27(1)
4.48(1)
3.95(1) 5.12 3.88 5'33(1)
4.15(1) 5.09 4.00
4.35 4.60(1) 4.39
a The root-mean-squared fluctuations from the MD simulations are given as the last significant digits in parenthess.
TABLE 8: GPC Glycerol Moiety Torsion Angle9 e1 e2 e3 expt A 172 -7 1 -63 expt B -63 61 -69 MDA (this work) 173(12) -68(12) -61(12) MD B (this work) 61(8) -73(11) -63(8) MDi6 62(7) -71(8) -62(7) a
e4 180
169 178(12) 164(11) 167(8)
Units are degrees. The root-mean-squared fluctuations from the
MD simulations are given in parentheses. TABLE 9: Lattice Parameters and Density of the DLG CrvstaP
expt 5.459
7.592 34.23 90.0 93.1 90.0 1.06 MD 5.548(2) 7.534(2) 34.08(1) 90.0(5) 92.6(3) 90.0(4) 1.05 " T h e MD values are averages over 40 ps. The MD root-meansquared fluctuations are given as the last significant digits in parentheses.
simulation and experimental results agree very well, within 3" for most and within 5" for all of the dihedral angles. The results from the constant NVE simulation of Stouch et al. are very similar to ours. We conclude that the glycerol region of GPC (and, hence, lecithin molecules) is modeled well by the Stouch potential. 3.2. 2,3-Dilauroyl-~-glycerol(DLG). The simulation of DLG was initiated from the crystal structure solved by Pascher and SundelLZs This molecule includes the glycerol moiety, the acyl ester linkages, and the hydrocarbon 'tails'. The crystal is monoclinic with two molecules per unit cell and one molecule in the asymmetric unit in the space group P21. The crystal structure can be described as bilayers with hydrocarbon chains packed parallel to each other at a 26.5" tilt angle with respect to the bilayer normal. The glycerol moieties are laterally linked within the bilayer plane by hydrogen bonds from the free hydroxyl group 011 to the carbonyl oxygen 0 3 2 in the ester linkage farthest from the hydroxyl in the b-translated molecules. There are no interlayer hydrogen bonds in the DLG crystal. Our analysis of the DLG crystal simulation is similar to that used for GPC. We begin with a comparison of the average lattice parameters and the density from the constant NPT simulation with the experimental data in Table 9. The a axis is 3% too long while the b and c axes are both less than 1% in error. The monoclinic symmetry was perfectly maintained during the simulation, and the p angle is less than 1" too large.
(A)
0.0217 0.0326 0.09
angle (deg)
dihedral (deg)
1.05 1.96 4.8
4.20 8.51
TABLE 11: DLG Chain Tilt Angle, Spacing, and Orientationa expt
MD
C2-chain
C3-chain
dCpc3
eb
26.2 27.1(2)
26.8 27.5(2)
4.7 4.7(0)
63.2 62.5(2)
@ 68.6 67.4(4)
a Units are 8, for distances and degrees for angles. The root-meansquared fluctuations from the MD simulations are given as the last significant digits in parentheses. e is the angle betweeen the interchain vector and the z axis. 4 is the angle between the projection of the interchain vector on the xy plane and the x axis.
The density from the simulation is 3% lower than the experimental density. Overall, the new hydrocarbon chain potential, combined with the model of Stouch et al. for the glycerollacy1 ester linkage, describes very accurately the lattice of the DLG crystal. In the X-ray crystal structure, the intralayer hydrogen bond length (0-0 distance) is 2.83 8,, whereas the average value from the simulation is 2.73 A. Using a distance-based criterion, we found no evidence for any interlayer hydrogen bonds during the simulation, in agreement with the crystal structure. In Table 10 we show the rms deviations from the experimental crystal structure of the DLG internal geometric parameters involving heavy atoms. For comparison we also include the deviations in the NVE simulation of Stouch et a1.I6 The deviations in the bond lengths and bond angles in both simulations are well within the experimental error, although the deviations in our simulation are significantly smaller than those of Stouch et al. The deviations in the torsion angles in our simulation are quite small and are about half as large as those from the simualtion of Stouch et al. The rms deviation of the average heavy atom positions with respect to the crystal structure was 0.38 8, in our constant NPT simulation. The hydrocarbon chain packing in the crystal is a feature we want to be able to reproduce well because errors in the chainpacking density will affect the molecular area and thus the headgroup network and because this, in tum, will lead to errors in the properties of simulated bilayers. Below the temperature of the chain melting (disordering) transition, if the area occupied by the headgroup is larger than the sum of the cross-sections of the two hydrocarbon chains, the chains tilt in order to establish close-packing contact. To analyze the chain packing, we have calculated the tilt angle (with respect to the bilayer normal), the distance between the centers of mass of the two chains, and the orientation of the vector connecting the chain centers of mass defined in terms of the polar angle 6 (with respect to the z axis) and the azimuthal angle q5 (between the projection onto the xy plane and the x axis). The comparison of the simulation and experimental data in Table 11 shows that the simulation quantitatively reproduces the chain packing in the DLG crystal: the average tilt and orientational angles from the simulation are generally within 1" of the experimental results, and the average distance between the chains is the same in the simulation and the crystal structure. The conformations of the glyceroUdiacy1ester group play an important role in determining the hydrocarbon chain packing. The torsion angles related to the glycerol C2-C3 bond, which constitutes the central link between the two hydrocarbon chains,
Glycerylphosphorylcholine and Dilauroylglycerol
J. Phys. Chem., Vol. 99, No. 24, 1995 10041
TABLE 12: DLG Glycerol Moiety Torsion AnglesO expt M D (this work) MDI6
81
82
83
-55 -65(8) -67(8)
61 56(8) 53(8)
173 175(7) 176(6)
4. Conclusions 84
54 54(7) 54(6)
a Units are degrees. The root-mean-squared fluctuations from the M D simulations are given in parentheses.
TABLE 13: DLG Acyl Ester Linkage Torsion Angle@ expt MD (this work) MD16
Pl 110 92(12) 89(10)
P2 170 168(8) 168(8)
P3 -159 -137(15) -137(13)
Y2
Y3
expt M D (this work) MDI6
Y’ 91 84( 1 1) 81(9)
P5 180
P4 176 179(8) 178(8)
177(10) 177(9) Y5
Y4 ~~
180 189(7) 190(7)
65 70( 11) 70(9)
180 183(7) 183(7)
-178 - 179(9) -179(7)
Units are degrees. The root-mean-squred fluctuations from the M D simulations are given in parentheses.
are antiperiplanarl-tsynclinal (in the nomenclature used by Hauser et al.,44 syn and anti denote torsion angles 90”, respectively, periplanar and clinal mean “approximately planar” and “inclined,” respectively, and +I- refers to the sign of the torsion angle). The antiperiplanar torsion angle 83, defined by Cl-C2-C3-031, guarantees the glycerol oxygen 031 continues in the zigzag plane of the glycerol carbon. Consequently, the natural D configuration at carbon atom C2 forces the glycerol oxygen atoms, 0 2 1 and 031, to which the fatty acids are attached, to lie in a +synclinal position. The orientation of the plane of the carboxyl groups toward the glycerol moiety is defined by the torsion angles pl and y l . The approximately anticlinal conformation of fi 1 orients the plane of the p carbonyl group perpendicular to the axis of the glycerol chain. The conformation of y l guarantees the perpendicular orientation of the two hydrocarbon chain zigzag planes. In order to achieve a parallel alignment of the chains, a bend in one or both of the chains is necessary. In the DLG crystal the p chain continues almost unbent (in the torsion p3) and the y chain forms a roughly perpendicular bend (in the torsions y 1 and y3). We list in Table 12 the dihedral angles of the glycerol moiety from the crystal structure and the average values from our simulation and that of Stouch et al. With the exception of 81, which is 10” off, our simulation results agree quite well with the experimental data. The agreement of the simulation of Stouch et al. is slightly worse. Note that the glycerol moiety torsion angles in GPC and DLG are quite different from each other (compare Tables 8 and 12), and the torsions in both GPC and DLG are significantly different from those in crystals of phosphatidylethanolamine and pho~phatidylcholine.~~ It is reassuring to see that the simulations can reproduce reasonably well the torsion angles in the two cases examined where the glycerol conformations differ substantially. Table 13 compares simulation and experimental results for the torsion angles in the acyl ester linkages. The largest deviations (about 20”) are in pl and p3. The smaller pl angle at the point of P-chain attachment to glycerol is compensated by the almost equally smaller p3. This leads to the correct orientation of the zigzag planes of the P-chain, Le., parallel to the y-chain. Likewise, the smaller, yet significant (about lo”) deviations in y l and y2 are offsetting and do not change the chain packing. Thus, in spite of some slight discrepancies between the simulated and experimental torsion angles in the glycerol and acyl ester region, as we have seen above, the hydrocarbon chain packing is correct in our DLG simulation.
We have performed constant N€T MD simulations of crystals of the lecithin fragments GPC, a model for the headgroup and glycerol regions, and DLG, a model for the glyceroVacy1 ester and hydrocarbon regions. By comparing the simulation results to experimental X-ray crystal structures, we have evaluted the quality of all-atom phospholipid force fields. For GPC we examined two existing force fields: the Stouch et a1.I6 and CHARMM 2229potentials. For DLG we employed a force field that is a combination of the Stouch et al. parameters for the g1yceroUacyl ester linkage and a compilation of new and existing potentials for the hydrocarbon chains. The average unit cell parameters and densities from the simulations of DLG and the Stouch et al. model for GPC were in good agreement, but those from the simulation of the CHARMM 22 model for GPC were in relatively poor agreement, with the experimental values. The molecular structures and hydrogen-bonding networks were also well-reproducedby the DLG simulation and the GPC simulation employing the Stouch et al. potential. The most significant deviations were in the torsion angles in the acyl ester linkages of DLG, but the deviations in adjacent torsions offset each other, leaving the hydrocarbon chain packing unchanged. In general, the results of our constant NPT simulations, in which the Ewald method was used to sum the electrostatic interactions, agree better with the experimental data than the constant energy and volume simulations of Stouch et al., in which the electrostatic interactions were truncated. We conclude that a combination of the Stouch et al. model for the headgroup and glycerollacy1 ester region and our compilation of potentials for the hydrocarbon chains appears to be a reasonable starting point for MD simulations of model biological membranes. Note Added in Proof, We have explored the consequences of using the second set of atomic charge parameters of ref 16 and find that there is little effect on the crystal properties of GPC and DLG. Acknowledgment. This work was supported by the National Institutes of Health under Grants R01 GM40712 and F32 GM14463. We thank Dr. A. MacKrell for providing us with the CHARMM 22 parameters before publication. References and Notes (1) Smith, G. S.; Sirota, E. B.; Safinya, C. R.; Clark, N. A. Phys. Rev. Lett. 1988, 60, 813. (2) Tristam-Nagle, S.; Zhang, R.; Suter, R. M.; Worthington, C. R.; Sun, W.-J.; Nagle, J. F. Biophys. J . 1993, 64, 1097. (3) Sun, W.-J.; Suter, R. M.; Knewtson, M. A,; Worthington, C. R.; Tristam-Nagle, S.; Zhang, R.; Nagle, J. F. Phys. Rev. E 1994, 49, 4665. (4) Alper, H. E.; Bassolino, D.; Stouch, T. R. J . Chem. Phys. 1993, 98, 9798. ( 5 ) Ahlstrom, P.; Berendsen, H. J. C. J . Phys. Chem. 1993,97, 13691. (6) Egberts, E.; Berendsen, H. 3. C. J . Chem. Phys. 1988, 89, 3718. (7) Damodaran, K. V.; Merz, K. M. Langmuir 1993, 9, 1179. (8) Damodaran, K. V.; Merz, K. M.; Gaber, B. P. Biochemistry 1992, 31, 7656. (9) Raghavan, K.; Reddy, M. R.; Berkowitz, M. L. Langmuir 1992,8, 233. (10) Venable, R. M.; Zhang, Y.: Hardy, B. J.; Pastor, R. W. Science 1993, 262, 223. (11) Heller, H.; Schaefer, M.; Schulten, K. J . Phys. Chem. 1993, 97, 8343. (12) Stouch, T. R. Mol. Simul. 1993, I O , 335. (13) Wilson, M. A.; Pohorille, A. J . Am. Chem. SOC.1994, 116, 1490. (14) Huang, P.; Perez, J. J.; Loew, G. H. J. Biomol. Struct. Dyn. 1994, I I , 921. (15) Egberts, E.; Marrink, S.-J.; Berendsen, H. J. C. Eur. Biophys. J. 1994, 22. (16) Stouch, T. R.; Ward, K. B.; Altieri, A,; Hagler, A. T. J. Comput. Chem. 1991, 12, 1033. (17) Williams, D. E.; Stouch, T. R. J . Comput. Chem. 1993, 14, 1066.
10042 J. Phys. Chem., Vol. 99, No. 24, 1995 (18) Nagle, J. Biophys. J . 1993, 64, 1476. (19) De Loof, H.; Harvey, S. C.; Segrest, J. P.; Pastor, R. W. Biochemistry 1991, 30, 2099. (20) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177. (21) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (22) Tobias, D. J.; Tu, K.; Klein, M. L. J . Phys. Chem., in press. (23) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (24) McCammon, J. A,; Harvey, S. C. Dyanmics of Proteins and Nucleic Acids; Cambridge University Press: New York, 1987. ( 2 5 ) Bareman, J. P.; Klein, M. L. J . Phys. Chem. 1990, 94, 5202. (26) Pearson, R. H.; Pascher, I. Nature 1979, 281, 499. (27) Abrahamsson, S.; Pascher, I. Acta Crystallogr. 1966, 21, 79. (28) Pascher, I.; Sundell, S. J . Mol. Biol. 1981, 153, 791. (29) MacKerell, A,; Schlenkrick, M.; Brickmann, J.; Karplus, M. Unpublished. (30) Discover molecular modeling system. BIOSYM Technologies, San Diego, CA. (31) Weiner, S. J.; Kollman, P. A,; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230.
Tu et al. (32) Smith, J. C.; Karplus, M. J . Am. Chem. Soc. 1992, 114, 801. (33) Ryckaert, J.-P.; McDonald, I. R.; Klein, M. L. Mol. Phys. 1989, 67, 957. (34) Williams, D. E. J . Chem. Phys. 1967, 47, 4680. (35) Ryckaert, J.-P.; Bellemans, A. J . Chem. Soc., Faraday Trans. 1978, 66, 95. (36) Tuckerman, M. E.; Martyna. G. J.; Berne, B. J. J . Chem. Phys. 1992, 97, 1990. (37) Rothlisberger, U.; Klein, M. L. Chem. Phys. Lett. 1994, 227, 390. (38) Andersen, H. C. J . Chem. Phys. 1980, 72, 2384. (39) Nos& S. J . Chem. Phys. 1984, 81, 511. (40) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. J . Chem. Phyr. 1992, 97, 2635. (41) Wentzcovitch, R. Phys. Rev. B 1991, 44, 2358. (42) Ciccotti, G.; Ryckaert, J.-P. Comput. Phys. Rep. 1986, 4, 345. (43) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J . Comput. Phys. 1977, 23, 327. (44) Hauser, H.; Pascher. I.; Pearson, R. H.; Sundell, S. Biochim. Biophys. Acta 1981, 650, 21. JP9422349
