Comparison between Molecular Dynamics Simulations and Wide

wide-angle X-ray scattering experiments. Though we find a generally good agreement between the simulated and experimental spectra, there are some ...
2 downloads 0 Views 81KB Size
Download PDF
2484

J. Phys. Chem. B 2007, 111, 2484-2489

Microscopic Structure of Phospholipid Bilayers: Comparison between Molecular Dynamics Simulations and Wide-Angle X-ray Spectra Marcello Sega* and Giovanni Garberoglio CNR-INFM and Department of Physics, UniVersity of Trento, Via SommariVe 14, I-38050 PoVo (Trento), Italy

Paola Brocca and Laura Cantu` Department of Chemistry, Biochemistry and Biotechnologies for Medicine, UniVersity of Milan, Via F.lli CerVi 93, I-20090 Segrate (MI), Italy ReceiVed: August 23, 2006; In Final Form: January 15, 2007

We present results of molecular dynamics simulations of fully hydrated dipalmitoylphosphatidylcholine and dimyristoylphosphatidylcholine bilayers in the disordered liquid crystalline phase (LR) and compare them to wide-angle X-ray scattering experiments. Though we find a generally good agreement between the simulated and experimental spectra, there are some deviations whose origin has been investigated by a reparametrization of the aliphatic chains’ force field. A detailed analysis of the various contribution to the X-ray spectra shows that a non-negligible contribution to the total scattered intensity comes from the headgroups and the headtail cross correlation.

1. Introduction In recent years, phospholipid membranes have been the subject of numerous investigations both from the experimental1-5 and computer simulation6-12 points of view. The interest in these artificial membranes resides in the fact that they are supposed to be good models for the study of many properties of the cell membrane. Many of the efforts in this field have been devoted to the study of dipalmitoylphosphatidylcholine (DPPC) and dimyristoylphosphatidylcholine (DMPC), mainly in the Lβ (gel) and LR (disordered liquid crystalline) phases, the latter being the most biologically relevant. A realistic description of the membrane properties via molecular dynamics (MD) simulation can be achieved by adopting a well tested force fieldssuch as CHARMM,13 AMBER,14 OPLS,15 or GROMOS16sa big enough system to minimize finite-size effects, and a long enough trajectory to reach the equilibration of the slowest degrees of freedom of interest. Up to now, none of the suggested parametrizations has emerged as the generally accepted description of these systems, although they have been used with some success. On the simulation side, this can be traced back to the large number of equally important parameters describing the molecule-molecule interactions that renders their optimization a highly nontrivial task as well as to the different methods that can be employed to perform the simulations, like, e.g., different thermodynamic ensembles or electrostatics treatments. Moreover, the experimental quantities to which computer simulations are compared may be characterized by large uncertainties as is the case, for example, of the area per lipid.17 In the investigation of lipid bilayers, the observables which are mostly used for comparison between computer simulations and experimental data are thermodynamic quantities such as the volume per lipid and “microscopic” quantities like the * Electronic address: [email protected].

deuterium order parameters of the lipid chains and the smallangle X-ray scattering (SAXS) spectra.18,19 Besides SAXS, wideangle X-ray scattering (WAXS) is also widely employed to investigate the microscopic structure of lipid bilayers. In fact, WAXS is an invaluable tool to obtain detailed information on the local organization of scattering centers, since the transferred momentum range accessible via this technique allows one to probe the system at distances of the order of the angstrom. Moreover, WAXS spectra are directly comparable to computer simulation results without strong assumptions intervening in the raw data treatment. However, to the best of our knowledge, a thorough comparison with WAXS data has never been performed, with the exception of the work by Klein and co-workers,20 who calculated the tail structure factor of the tails of DPPC in a bilayer in the Lβ phase and showed its similarity to WAXS data. No such analysis seems to be available for phospholipids in the LR phase, which characterizes the membranes of the living cell, and no attention has been devoted to the headgroups, which are commonly thought to be completely disordered20 and therefore not contributing significantly to the X-ray spectrum intensity. In this work, we compare the results of novel WAXS experiments with simulation results obtained by us using one of the mainstream potential models, namely that proposed by Berger and co-workers.9 Moreover, we investigate the contribution of the tail and headgroup moieties to the total scattering intensity and their role in determining the features of the main WAXS peak. 2. Methods 2.1. Experimental Setup. DMPC and DPPC were purchased from Avanti Polar Lipids Inc. (Alabaster, AL) and used without further purification. DMPC and DPPC samples for X-ray scattering measurements were prepared in pure water at a concentration of 30% bw. Samples were put in plastic capillaries (KI-BEAM, ENKI srl)

10.1021/jp065450d CCC: $37.00 © 2007 American Chemical Society Published on Web 02/22/2007

Microscopic Structure of Phospholipid Bilayers with a 2 mm internal diameter, 0.05 mm wall thickness, and very high X-ray transmission (98%), closed with polyethylene caps. The sample holder was temperature controlled by means of Peltier elements. Contemporary SAXS and WAXS measurements were performed on the ID02 high brilliance line at ESRF (Grenoble, France) with a beam cross section of 0.3 mm × 0.8 mm and a wavelength of 0.1 nm. The wavevector range obtained for WAXS spectra was 6-49 nm-1, and irradiation time was 0.1 s. Empty plastic capillaries and water samples were used to obtain the experimental subtracted spectra. Analogous DMPC and DPPC samples were prepared at 2% bw concentration in degassed water to perform accurate densitometry. Density measurements were performed with an Anton Paar DMA 5000 instrument. Measurements were taken after thermal equilibration and signal stability. 2.2. Simulation Details. The simulations have been performed using the open source GROMACS simulation package,21,22 describing the molecules at atomistic detail, except from the apolar CH, CH2, and CH3 groups, for which the united atom representation has been employed. The rectangular simulation cells, with typical volumes around 250 nm3, contained 128 DPPC or DMPC molecules and 3655 water molecules, with periodic boundary conditions applied in all directions. The water molecules have been simulated using the SPC23 model, that has been shown to be a valid choice for reproducing the properties of water at interfaces.24 As a starting point, we have used the potential parameters proposed by Berger and co-workers9 with the charges on the headgroup parametrized by Chiu and co-workers.25 We used the Ryckaert-Bellemans26 model for the dihedral interactions on the aliphatic chains. Two Berendsen thermostats with a time constant of 0.1 ps have been employed to maintain the temperature fixed for water and DPPC/DMPC molecules separately. The DPPC simulation has been performed at T ) 323 K, whereas the DMPC simulation has been performed at T ) 333 K, that is, the same temperatures as the experimental data with which we perform our comparison. The three diagonal components of the pressure tensor have been kept fixed at 1 atm by coupling separate Berendsen barostats to each of the x, y, and z directions, employing a time constant of 1.0 ps. As initial structures, we used the final structure of run E discussed in ref 27 for DPPC, available at http://moose.bio.ucalgary.ca/, and the final structure described in ref 28 for DMPC, available at http://www.apmaths.uwo.ca/∼mkarttu/. All bond lengths have been constrained using the SETTLE29 and LINCS30 algorithms for water and lipid molecules, respectively, thus allowing us to integrate the equations of motion with a time step of 2 fs. Lennard-Jones interactions have been calculated using a simple cutoff at 1.0 nm, as in the paper by Berger, and long-range corrections have been applied to calculate energies and pressures.31 The electrostatic interactions have been calculated using the particle-mesh Ewald method (PME),32 with a mesh grid of 0.12 nm and a spline order of 4. This setup will be called “simulation I/DPPC” and “simulation I/DMPC” for the two systems that we have considered. All the phospholipidic systems have been further equilibrated for 2 ns, checking the relaxation of various observables, such as volume, pressure, and potential energy. The spectra and other quantities have been calculated using production runs of 3 ns under the same thermodynamic conditions, averaging over configurations sampled every 20 ps. We have used the g_rdf program of the GROMACS package to calculate the spectra, which we have modified in order to compute also the partial contributions to be discussed later.

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2485 TABLE 1: Comparison between the Experimental and Simulated Volume and Area Per Lipidf simulation

volume per lipid (nm3)

I/DPPC I/DMPC II/DMPC DPPC (exp) DPPC (exp) DMPC (exp) DMPC (exp)

1.176 ( 0.004 1.075 ( 0.004 1.076 ( 0.004 1.232 ( 0.001a 1.229 ( 0.001c 1.127c 1.160 ( 0.020d

area per lipid (nm2) 0.647 ( 0.009 0.656 ( 0.009 0.652 ( 0.013 0.64 ( 0.01b 0.6633 ( 0.0083e

a Nagle et al.34 b Nagle et al.17 c Extrapolated from Nagle and Wilkinson.35 d This work. e Balagavy´ et al.36 f DPPC data refer to T ) 323 K, and DMPC data refer to T ) 333 K.

We have also performed another simulation of the DMPC membrane, called “simulation II/DMPC”, reparametrizing the force field in the same manner as done by Berger in the case of DPPC, trying to investigate in detail the origin of the small discrepancies observed between the experimental data and the results of simulation I/DMPC. In particular, we investigated the role of the tails’ parameters in determining the DMPC WAXS spectrum. In the original paper by Berger and co-workers,9 in order to discover the optimal parameters for the 15 carbon groups in each of the DPPC tails, simulations of a model system for the tails, namely pentadecane (CH3(CH2)13CH3), were performed, varying the Lennard-Jones interactions until the heat of vaporization and density of actual pentadecane were satisfactorily reproduced. Analogously, we followed the same procedure to find another possible set of parameters for DMPC’s tails, using tridecane (CH3(CH2)11CH3), as our model system. To this aim, we ran a set of 50 simulations of liquid- and gas-phase tridecane, varying the Lennard-Jones parameters to find which one would reproduce the experimental heat of vaporization and density of tridecane at T ) 298 K and p ) 1 atm.33 Following Berger, we assumed the Lennard-Jones diameters (σCH2, σCH3) of the CH2 and CH3 groups to be the same and the respective interaction energies, CH2, CH3, to be related by the formula CH3 ) (3/2)CH2. The simulations of the liquid phase were performed using a system of 128 tridecane molecules. For each of the new Lennard-Jones parameters, a 100 ps equilibration run was followed by a 200 ps production run. The gas-phase simulations consisted in 5 ns long runs of a single tridecane molecule, in the absence of periodic boundary conditions. Both the liquid and gas-phase simulations were performed at the constant temperature of T ) 298 K, whereas the liquid phase has been simulated at the constant pressure of 1 atm. Pressure and temperature have been kept fixed using the Berendsen weak coupling scheme, and the long-range corrections to energy and pressure have been included explicitly. 3. Results and Discussion 3.1. Volume and Area per Lipid. In order to compare our results with the available data and other simulations, we have computed the volume and area per lipid for DPPC and DMPC bilayers, which we report in Table 1, together with the experimental results. The volume per lipid VL has been computed as

VL )

Vbox - NWVW Nlip

(1)

where Vbox is the average volume of the simulation box, NW ) 3655 is the number of water molecules, and Nlip ) 128 is the

2486 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Sega et al.

number of lipids. The volume of a water molecule, VW, has been calculated with independent bulk water simulations at the same state points, obtaining the values VW ) 0.03122 and 0.03152 nm3 at T ) 323 and 333 K, respectively. Similarly, the area per lipid has been calculated as

AL )

Axy Nlip/2

(2)

where Axy is the average value of the simulation box surface parallel to the bilayer. The computer simulation results are sensibly below the experimental values for DPPC but not significantly different from the results obtained by other computer simulations on the same systems, such as the value given by Berger,9 VL ) 1.189 ( 0.009 nm3. 3.2. DPPC WAXS Spectra. The X-ray scattering intensity in the Born approximation can be written as follows:

I(q) ∝ 〈|

fj(q)eiq‚r ∑ j,m

|2〉

j,m

(3)

where rj,m is the position of the jth atom on the mth molecule and the symbol 〈‚‚‚〉 indicates an ensemble average. The overbar indicates an average over all the possible relative orientations of the wavevector q with respect to the simulation cell, to take into account the fact that the experimental samples consist in vesicles, that we approximate as unoriented bilayers. The quantity fj(q) is the form factor of atom j and is computed using the formula 4

f(q) ) c +

[ ( )] q

ai exp -bi ∑ 4π i)1

2

(4)

where c, ai, and bi are the Cromer-Mann coefficients for the given atomic species.37 To compute the X-ray scattering intensity (3), the following procedure has been employed: Since the equilibrium run is carried on at constant pressure, fluctuations in the box edges require a binning of the reciprocal space. The bin width ∆q is chosen to be ∆q ) 2π/L, being L the largest box vector of the first frame being analyzed. For every frame chosen for the sampling, then the following procedure is used: (a) The discrete set S of the allowed values of q is computed, according to the varying length of box edges and up to a predefined cutoff value. The real and imaginary parts of the sum Σ(q) ) ∑j,mfj(q)eiq‚rj,m appearing in eq 3 are hence computed for each of the wavevectors in S. (b) The angular average indicated with an overbar in eq 3 is then computed in the following way. For each of the vectors in S, its modulus q is computed and its corresponding bin assigned. The value of the histogram I(q) in that bin is then incremented by the modulus of Σ(q), normalized by the number of vectors that give contribution to the bin. The code for the computation of the scattering intensity as described above is part of the GROMACS simulation package, although we used a slightly modified version to allow the computation of the intra- and intermolecular contributions, as defined below, which is made available upon request. In presenting the WAXS experimental data, one subtracts the spectra of both the vessel and solvent, to obtain a pattern which is representative of the scattering intensity of the solute alone. In order to compare the simulated spectra with the experimental

Figure 1. Comparison between simulated and experimental spectrum of DPPC at 323 K: (solid line) experimental data; (dashed line) simulation result using Berger interaction potential (interpolated); (inset) the various contribution to the simulated WAXS spectrum for DPPC, according to eq 5; (squares) head contribution; (dashed line) tail contribution; (solid line) head-tail cross contribution.

ones, we have removed in the same way the contribution of the solvent from the computer generated spectra by subtracting the scattering intensity computed from a separate run of bulk water at the same state point. The simulation boxes of bulk water and lipid solution differ in size, and hence, also the accessible points in the reciprocal space are different. Therefore, before subtracting the two spectra, an interpolation using cubic splines has been performed. Finally, since the scattering intensity in eq 3 is defined up to a factor, in order to compare it with the experimental one, the computed spectra have been scaled so that the experimental and simulated intensities coincide at the position of the peak. We note in passing that the spectrum obtained by this subtraction is not the same as the one that would be obtained by restricting the sum to the phospholipids’ atoms only in eq 3, due to the phospholipid-water cross-contributions, which cannot be accessed in an actual experiment. We present in Figure 1 the simulated and experimental WAXS spectra of DPPC at T ) 323 K. The position of the well-known peak centered around q0 ) 14.3 nm-1 is well reproduced by the MD simulation as well as the general shape of the WAXS spectrum. The simulation data result in a slightly larger peak width, indicating a somehow shorter range of the order than in the actual sample. The disagreement at low wavevectors that can be noticed in the figure can be ascribed to the fact that the experimental data are taken on vescicles whereas our simulation box describes a flat system, although the isotropy of the real system is taken approximatively into account by using the spherical average on the wavevectors described above. The analysis of the WAXS spectrum alone does not give direct information on the spatial correlations between the lipid molecules, which is better addressed by the analysis of the moiety-moiety correlation function NtNt′ NmNm′

Itt′(q) ) 〈

ftj(q)ft′j′(q)eiq‚(r ∑ ∑ mm′ jj′

tmj-rt′m′j′)

〉

(5)

where the indices t and t′ indicate a particular class of molecules or moieties present in a sample (for example, t may represent the DPPC molecules, and t′, the water molecules), the indices m and m′ run on the different elements of the classes indicated by t and t′ (in our example Nt ) 128 and Nt′ ) 3655), and,

Microscopic Structure of Phospholipid Bilayers

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2487

finally, j and j′ run on the atoms composing the groups m and m′, respectively. The overline is the usual average over the orientation of the vector q, and 〈‚‚‚〉 is the average over all the configurations. It is straightforward to show that the total scattered intensity defined in eq 3 can be written as

I(q) )

Itt′(q) ∑ tt′

(6)

if the choice of the classes t is such that each of the atoms belongs to one and only one of the classes t (as in the case of the example given above). This property can be used to analyze the contribution to the total scattering intensity I(q) from the different moieties present in a sample. In the inset of Figure 1, we show the patterns corresponding to the cases where t ) t′ identifies either the atoms in the headgroup or the ones in the tail of DPPC. We also show the cross contribution where t represents the headgroup’s atoms and t′ represents the tail’s atoms, namely Ihead head(q), Itail tail(q), and 2Ihead tail(q). This partial decomposition shows that the WAXS peak is mostly, but not totally, due to the correlations between atoms in the hydrocarbon chains. The contribution due to the heads is dominant in the q < 11 nm-1 region and shows a well defined shoulder at the peak position, where it contributes with half of the intensity of the tails. Even though the headgroups are less structured than the tails, as already pointed out in ref 1, their contribution to the total scattering intensity is non-negligible. Moreover, it is shown that there is another important contribution, due to the head-tail cross contribution. As in the case of the heads, the intensity is greater in the low q region, although in this case a clear peak developes at q0. 3.2.1. Order Parameter. The spectral contribution from a given class t of atoms, Itt(q) can be further expanded as the sum of the contributions coming from the intra- and intermolecular (or moiety) correlations45

Itt(q) ) Itintra(q) + Itinter(q)

(7)

where the intramolecular part is defined as Nt

|

|2〉

tmj

Nt

tmj-rtm′j′)

〉

(9)

In order to allow a comparison of the intermolecular contributions from groups with different electronic densities, it is useful to normalize Itinter(q) to the scattering power of the single group at the given wavevector. This leads to the definition of the structure factor St(q) as St(q) ) Itinter(q)/Itintra(q) + 1. The quantity St(q) - 1 acts as an order parameter for the group t since it becomes zero when no spatial correlation between the given group on different molecules is found. We have calculated the order parameters St(q) - 1 for the whole DPPC, the tail, and the head moieties, which are reported in Figure 2. It is apparent from the figure that the DPPC order parameter takes the greatest contribution from the component coming from the tail-tail correlation, which dominates the peak of the

2π q0

(10)

obtaining for DPPC the value d ) 0.54 nm. Equation 10 can be derived from the Zernike-Prins formula38-41

S(q) ) 1 + F

Nm

∑ ∑ ∑ ftj(q)ftj′(q)eiq‚(r m m′*m jj′

Itinter(q) ) 〈

d ) 1.2285

(8)

and the contribution coming from the intermolecular correlations is Nt

structure factor observed at around 15 nm-1 and develops welldefined oscillations for larger wavevectors. This behavior is similar to that observed in the case of simple liquids31 and thus enforces the picture of a highly structured and close packed intermolecular organization of the tail groups. Turning to the case of the headgroups, one can observe that the shoulder found in the WAXS spectrum (see the inset of Figure 1) corresponds to a well defined peak in the order parameter. This peak is located at slightly larger wavevectors (i.e., shorter distances) with respect to the tails’ peak and is also approximatively three times smaller. Another signature of the little degree of order in the headgroups can be observed in the absence of defined oscillations in the order parameter after the first peak. Nevertheless, this finding shows that the headgroups are indeed structured, although less than the tails, and that the absence of a clear peak in the WAXS pattern does not necessarily imply the absence of local order. 3.2.2. Real Space Correlations. Since the major contribution to the scattering intensity is due to the tails, the position of the peak can be used to obtain an estimate of the average distance d between the scattering centers of the tails using the relation

Nm

ftj(q)eiq‚r ∑ ∑ m j

Itintra ) 〈

Figure 2. Order parameter for simulated DPPC at T ) 323 K using Berger’s force-field: (dotted line) total correlation; (solid line) tailtail contribution; (dashed line) head-head contribution

∫0∞[g(r) - 1]

sin(qr) 2 4πr dr qr

(11)

in the case of an ideal liquid composed of spherical particles in three-dimensional close packing. This is different from the actual practice in the case of gel-phase phospholipids.46 The quantity g(r) appearing in the previous equation is the pair distribution function, defined as

g(r) )

1 4πNFr2

N

〈

δ(r - |ri - rj|)〉 ∑ i,j

(12)

where the indices i and j run over all the N carbon groups of the tails. For a homogeneous bulk liquid, the constant F appearing in eqs 11 and 12 is just the density, and as a consequence, g(r) tends to the value of one for large values of the distance r.

2488 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Sega et al.

Figure 3. Radial distribution function of the CH2 and CH3 groups of DPPC tails: (dashed line) complete radial distribution function; (solid line) intermolecular contribution to the radial distribution function. We subtracted the delta function like contribution to the intramolecular term coming from the covalently bonded carbon groups.

Figure 4. WAXS spectrum for DMPC at 333 K: (full line) experimental data, (dashed line) results for simulation I/DMPC (interpolated), (dotted-dashed line) results for simulation II/DMPC (interpolated), (inset) the peak region.

However, in the case of DPPC, the tails’ carbon groups are inhomogeneously distributed within the simulation box and there is freedom in the choice of the value of F. We decided to set it to the density of the tail groups within the volume occupied by the lipid bilayer alone, i.e., F ) (N/NlipVL), where N denotes in this case the number of carbon groups in our system. With this choice, the meaning of g(r) is the ratio between the average local density at a distance r from a group in the lipids’ tails and the reference density of the tails’ groups within the bilayer, F. For distances larger than the bilayer thickness, it is then expected for g(r) to become less than one. The radial distribution function can be naturally separated into the contributions coming from carbon groups belonging to the same (intra) or to different (inter) molecules. The intermolecular contribution can be computed by subtracting from g(r) the intramolecular radial distribution function, calculated from eq 12 by restricting the sum over the groups of a single tail and averaging over all the molecules in the sample. The total and intermolecular contribution to the radial distribution function are reported in Figure 3. We notice that the intermolecular contribution is the leading one and has a significant peak at 0.53 nm, in perfect agreement with the result of eq 10, therefore confirming the picture of a liquidlike close packed distribution of the carbon groups in the tails as outlined above. 3.3. DMPC WAXS Spectra. We have performed the calculation of the WAXS spectrum also in the case of DMPC in order to test to what extent the force field developed for DPPC is transferable to other lipid membranes. The results, shown in Figure 4, indicate that one still has a generally good agreement between the simulated and the measured WAXS spectra. However, in this case, we notice that the peak at q0 ) 14.4 nm-1 is not reproduced as well as in the case of DPPC, and the whole pattern appears to be shifted toward lower wavevectors by an amount of ∼0.25 nm-1. We tried to see if a reparametrization of the force field would account for the observed discrepancy between the WAXS spectra. We decided to focus on the tails’ parameters because this part of the lipid is responsible for most of the measured intensity, as found above. Since the lipid tails of DMPC have 13 carbon groups (as opposed to the 15 in the case of DPPC), we have chosen tridecane as a model system for them. We then set up, following

Berger’s procedure, 50 bulk liquid and gas-phase tridecane simulations varying the Lennard-Jones parameters as described in section 2.2 and compared the results of these simulations with the experimental values of heat of vaporization Hvap and molecular volume Vmol taken from ref 33, in order to find the optimal values of the Lennard-Jones parameters. The molar heat of vaporization has been computed as ∆Hvap ) Ugas - Uliq + RT, where Ugas and Uliq are the molar potential energies of the gas and liquid phase simulations, respectively, and R is the gas constant. This procedure resulted in the following values: CH2 ) 0.4 kJ/mol and σCH2 ) 0.42844 nm, i.e., about 12% and 1% more than the pentadecane parameters found by Berger,9 respectively, corresponding to a heat of vaporization of 66.1 kJ/mol and an average volume per molecule of 0.406 nm3. These values have to be compared with the experimental ones reported in ref 33, namely Hvap ) 66.23 kJ/mol and Vmol ) 0.40674 nm3, and to those obtained using Berger’s parameters to model tridecane, namely Hvap ) 54.8 kJ/mol and Vmol ) 0.4133 nm3. It is known42,43 that keeping a constant σCH3/σCH2 ratio for alkanes of different chain length results in chain-length dependent Lennard-Jones parameters when the optimization is made to reproduce thermodynamic quantities, such as the average density of the enthalpy of evaporation. We would like to point out that the main reason for this reparametrization is to assess whether and how much the spectrum depends on the LennardJones values assigned to the tails’ group in the particular case of the DMPC bilayer only and it is hence not necessary to use transferable parameters, such as the ones developed in ref 43. We present in Table 2 a comparison between different set of Lennard-Jones parameters for aliphatic carbon chains. We report the results for the volume and area per lipid obtained with this new set of parameters in Table 1. We notice that the reparametrization did not alter the observed values within the indetermination. The WAXS spectrum of a simulation of DMPC with this new force-field is also reported in Figure 4. We notice that despite this refinement the position of the main peak is almost unaltered. This is again an indication of the fact that the structure of the lipid membrane is not determined by the properties of the tails alone, even in the spectral region dominated by the contribution of the tails themselves, and therefore, a complete reparametrization of the force field, which is outside the scope

Microscopic Structure of Phospholipid Bilayers

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2489

TABLE 2: Comparison between the Force-Field Parameters for the CH2 and CH3 Groups of DPPC and DMPCe force field Bergera 45A3b Chiu et al.c II/DMPCd

group

σ



C6 × 103

C12 × 105

CH2 CH3 CH2 CH3 CH2 CH3 CH2 CH3

0.39600 0.39600 0.40700 0.37480 0.405 0.375 0.40000 0.40000

0.38000 0.57000 0.41050 0.86720 0.42 0.92 0.42844 0.64007

5.8740 8.7770 7.4684 9.6138 7.413727 9.121396 7.0148 10.4816

2.2650 3.3850 3.3966 2.6646 3.271726 2.536583 2.8713 4.2911

a

Reference 9. b Reference 33. c Reference 43. d This work. e The units for σ, , C6, and C12 are nm, kJ mol-1, kJ nm6 mol-1, and kJ nm12 mol-1, respectively.

of this paper, would be necessary in order to reproduce all the details of the WAXS spectrum. We would like to stress, however, that the reparametrization has been performed only to look for a possible explanation of the little discrepancy between experimental and simulated data and is meant neither to support the idea of chain-length dependent parameters nor to suggest them as the best possible choice. 4. Conclusions In this paper, we have presented a thorough comparison between the results of wide-angle X-ray scattering experiments and computer simulations of DPPC and DMPC fully hydrated bilayers. The force field that has been developed by Berger and coworkers to describe DPPC has been shown to agree reasonably well with the experimental data. The observed position of the WAXS peak is reproduced within 1% and 2% in the case of DPPC and DMPC, respectively. The overall shape of the computer simulation peak resulted to be slightly larger with respect to experiment in the case of DPPC, and shifted toward lower wavevectors by an amount of ∼0.25 nm-1 in the case of DMPC. In order to understand the difference between the two cases we reparametrized the DMPC force field to take into account the shorter length of DMPC’s hydrocarbon chains with respect to that of DPPC. This reparametrization procedure provided a set of Lennard-Jones parameters that correctly reproduce both the heat of vaporization and the density of liquid tridecane, the system chosen as reference for the DMPC tails. We observed no sensible change in the peak position usign these new parameters, and we conclude that the peak position of the spectrum does not depend on the details of the interaction parameters of the tails’ groups alone, despite the fact that the scattered intensity is mainly due to the tails’ centers in that region. A detailed analysis of the various contributions to the WAXS peak of both systems showed that the intensity of the spectrum is mostly due to the interchain correlation of the carbon atoms on the lipid tail. However, the scattering intensity does not measure the tail structure only but takes also non-negligible contributions from the head-head and head-tail correlations. We then conclude that a precise reparametrization of the force field would have most likely to take into account all the different molecular groups of the lipid. Acknowledgment. M.S. and G.G. would like to thank Prof. R. Vallauri for useful discussions on the subject of this paper and a careful reading of the manuscript. The calculations have been performed on the HPC facility at the Physics Department of the University of Trento.

References and Notes (1) Sun, W. J.; Suter, R.; Knewtson, M.; Worthington, C.; TristramNagle, S.; Zhang, R.; Nagle, J. Phys. ReV. E 1994, 49, 4665. (2) Nagle, J. Biophys. J. 1993, 64, 1476. (3) Kucˇerka, N.; Liu, Y.; Chu, N.; Petrache, H.; Tristram-Nagle, S.; Nagle, J. Biophys. J. 2005, 88, 2626. (4) Mason, P.; Gaulin, B.; Epand, R.; Wignall, G.; Lin, J. Phys. ReV. E. 1999, 59, 3361. (5) Lemmich, J.; Mortensen, K.; Ipsen, J.; Hønger, T.; Bauer, R.; Mouritsen, O. Phys. ReV. E 1996, 53, 5169. (6) Egberts, E.; Marrink, S. J.; Berendsen, H. Eur. Biophys. J. 1993, 22, 423. (7) Feller, S.; Venable, R. M.; Pastor, R. Langmuir 1997, 13, 6555. (8) Moore, P.; Lopez, C.; Klein, M. Biophys. J. 2001, 81, 2484. (9) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 2002. (10) Jedlovszky, P.; Mezei, M. J. Chem. Phys. 1999, 111, 10770. (11) Marrink, S.; Lindahl, E.; Edholm, O.; Mark, A. J. Am. Chem. Soc. 2001, 123, 8638. (12) Petrache, H.; Feller, S.; Nagle, J. Biophys. J. 1997, 70, 2237. (13) Brooks, B.; Bruccoleri, R.; Olafson, B.; States, D.; Swaminathan, S.; Karplus, M. J. Comp. Chem. 1983, 4, 187. (14) Weiner, S.; Kollman, P.; Case, D.; Singh, U.; Ghio, C.; Alagona, G. J. Am. Chem. Soc. 1984, 106, 765. (15) Jorgensen, W.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (16) van Gunsteren, W.; Berendsen, H. Groningen Molecular Simulation (GROMOS) Library Manual; Biomos: Groningen, 1987. (17) Nagle, J.; Tristram-Nagle, S. Biochim. Biophys. Acta 2000, 1469, 159. (18) Pabst, G.; Koschuch, R.; Pozo-Navas, B.; Rappolt, M.; Lohner, K.; Laggner, P. J. Appl. Crystallogr. 2003, 36, 1378. (19) Klauda, J.; Kucˇerka, N.; Brooks, B.; Pastor, R.; Nagle, J. Biophys. J. 2006, 90, 2796. (20) Tu, K.; Tobias, D.; Blasie, J.; Klein, M. Biophys. J. 1996, 70, 595. (21) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Mod. 2001, 7, 306. We have actually used version 3.2.1 of the package. (22) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comp. Phys. Commun. 1995, 91, 43. (23) Berendsen, H.; Postma, J.; van Gunsteren, W.; Hermans, J. Intermolecular Forces; Reidel: Dordrecht, 1981, p. 331. (24) Tieleman, D.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871. (25) Chiu, S.; Clark, M.; Balaji, V.; Subramaniam, S.; Scott, H.; Jakobsson, E. Biophys. J. 1995, 69, 1230. (26) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1997, 18, 1463. (27) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871. (28) Gurtovenko, A.; Patra, M.; Karttunen, M.; Vattulainen, I. Biophys. J. 2004, 86, 3461. (29) Miyamoto, S.; Kollman, P. A. J. Comp. Chem. 1992, 13, 952. (30) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comp. Chem. 1997, 18, 1463. (31) Allen, M.; Tildesley, D. Computer Simulations of Liquids; Clarendon Press: Oxford, 1987. (32) Essman, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (33) Schuler, L.; Daura, X.; van Gunsteren, W. J. Comp. Chem. 2001, 22, 1205. (34) Nagle, J.; Wiener, M. Biochim. Biophys. Acta 1988, 942, 1. (35) Nagle, J.; Wilkinson, D. Biophys. J. 1978, 23, 159. (36) Balagavy´, P.; Dubnicˇkova´, M.; Kucˇerka, N.; Kiselev, M.; Yaradaikin, S.; Uhrı´kova´, D. Biochim. Biophys. Acta 2001, 1512, 40. (37) Hahn, T., Ed. International tables for crystallography; Reidel: Dordrecht, 1991; Vol. C, pp 500-502. (38) Keesom, W.; de Smedt, J. Proc. Acad. Sci. Amsterdam 1922, 25, 118. (39) Keesom, W.; de Smedt, J. Proc. Acad. Sci. Amsterdam 1923, 26, 112. (40) Zernike, F.; Prins, J. A. Z. Phys. 1927, 41, 184. (41) Warren, B. E. Phys. ReV. 1933, 44, 969. (42) Chiu, S.; Clark, M.; Jakobsson, E.; Subramaniam, S.; Scott, H. J. Phys. Chem. B 1999, 103, 6323. (43) Chiu, S.; Vasudevan, S.; Jakobsson, E.; Mashl, R.; Scott, H. Biophys. J. 2003, 85, 3624. (44) Tardieu, A.; Luzzatti, V.; Reman, F. J. Mol. Biol. 1973, 75, 711. (45) For the sake of simplicity, we will use from now on the term “intramolecular” only. (46) In the lipid gel phase, it is customary to consider the main peak as originated from the (x, y) ordering of the lipid chains extending in the z direction across the bilayer. Thus, on each plane cutting the bilayer parallel to the surface, the chains leave an hexagonal trace with a lattice parameter d associated to the peak position q0 by the relation d ) (2π/q0)2/x3 (ref 44).