Simulation-Based Methods for Interpreting X-Ray Data from Lipid Bilayers

2796 Biophysical Journal Volume 90 April 2006 2796–2807 Simulation-Based Methods for Interpreting X-Ray Data from Lipid Bilayers Jeffery B. Klauda,* Norbert Kucˇerka,y Bernard R. Brooks,* Richard W. Pastor,z and John F. Nagley *Laboratory of Computational Biology, National Institutes of Health, Bethesda, Maryland 20892; yPhysics and Biological Sciences Departments, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213; and zLaboratory of Biophysics, Center for Biologics Evaluation and Research, U.S. Food and Drug Administration, Rockville, Maryland 20852-1448 ABSTRACT The fully hydrated liquid crystalline phase of the dimyristoylphosphatidycholine lipid bilayer at 30
C was simulated using molecular dynamics with the CHARMM potential for five surface areas per lipid (A) in the range 55–65 A˚2 that brackets the previously determined experimental area 60.6 A˚2. The results of these simulations are used to develop a new hybrid zero-baseline structural model, denoted H2, for the electron density profile, r(z), for the purpose of interpreting x-ray diffraction data. H2 and also the older hybrid baseline model were tested by fitting to partial information from the simulation and various constraints, both of which correspond to those available experimentally. The A, r(z), and F(q) obtained from the models agree with those calculated directly from simulation at each of the five areas, thereby validating this use of the models. The new H2 was then applied to experimental dimyristoylphosphatidycholine data; it yields A ¼ 60.6 6 0.5 A˚2, in agreement with the earlier estimate obtained using the hybrid baseline model. The electron density profiles also compare well, despite considerable differences in the functional forms of the two models. Overall, the simulated r(z) at A ¼ 60.7 A˚2 agrees well with experiment, demonstrating the accuracy of the CHARMM lipid force field; small discrepancies indicate targets for improvements. Lastly, a simulation-based model-free approach for obtaining A is proposed. It is based on interpolating the area that minimizes the difference between the experimental F(q) and simulated F(q) evaluated for a range of surface areas. This approach is independent of structural models and could be used to determine structural properties of bilayers with different lipids, cholesterol, and peptides. INTRODUCTION Numerous studies over many years have focused on refining parameters restricts applications to systems at low hydration. the structure of lipid bilayers (1–5). However, quantitatively The structural model developed by Nagle and co-workers accurate structures of even pure bilayers have been difficult (8), here denoted the hybrid baseline model (HB), falls be- to obtain, especially in the most biologically relevant liquid tween the previous two. Specifically it consists of two func- crystalline (La) phase consisting of disordered and fully hy- tional types: Gaussians representing the lipid headgroups and drated lipids. Such bilayers are not crystalline with atomic the terminal methyls; and a baseline function consisting of positions determined at the A ˚ ngstro¨m level, but have atomic strips representing water and the methylene plateau joined by distribution functions with widths spread over 5 A ˚ . This pre- a smooth bridging function. With additional assumptions cludes an atomic-level structural description and substan- and data, HB also yields the surface area per lipid, A. Given tially limits the quality and quantity of structural data that can that most molecular simulation or modeling studies require A be obtained. Consequently, structural models are required to for at least the initial condition, the importance of this feature elucidate structural quantities in real space (e.g., electron in a structural model is clear. density profiles, surface areas/lipid, component densities) from An awkward aspect of HB involves the baseline function: experimental observations in reciprocal space, i.e., the scat- the electron density in the superposition region is due not tering form factors, F(q) (6,7). Broadly stated, a structural only to the water and the hydrocarbon chain methylenes, but model specifies the form of the electron density profile, r(z), also to the headgroup components. A more transparent model and the specific values of the parameters are determined by has no baseline function and more simply represents both the fitting to experiment. methylenes and the water by separate functions; such a model A variety of structural models have been applied to mem- has been advocated for reflectometry studies of monolayers branes. Wilkins et al. (5) employed constant electron densi- (9). The first part of this article develops a hybrid zero- ties for different regions to obtain the electron density of the baseline model, denoted H2, for analyzing and interpreting La phase, but the physically unrealistic discontinuities at the diffraction data from bilayers. edges of regions lead to spurious large amplitude high q The approach employed here is simulation based. Results oscillations in F(q) (5,7,8). The structural model of Wiener from molecular dynamics (MD) simulations have been fruit- and White (4) consisting exclusively of Gaussians is not fully compared with those from diffraction experiments. For confounded by discontinuities, but the large number of free example, Feller et al. (10) demonstrated that the distributions of certain lipid component groups were not Gaussian. More recently, Sachs et al. (11) and Benz et al. (12) compared the Submitted October 8, 2005, and accepted for publication December 29, 2005. simulations with experiment for various molecular properties Address reprint requests to Jeffery B. Klauda, E-mail: [email protected]. in real and Fourier space. Here the application uses simulations  2006 by the Biophysical Society 0006-3495/06/04/2796/12 $2.00 doi: 10.1529/biophysj.105.075697 Interpreting X-Ray Data from Bilayers 2797 to help motivate the functional form of H2, and to provide straints are required and to estimate the level of confidence in test cases for comparing HB and H2. The obvious advantage obtaining A); v), ‘‘Application to experimental x-ray data’’ of testing models with simulations is that much more detailed (comparison of density profiles obtained by H2 and HB from structural information is available from a simulation than the experimental F(q)); and vi), ‘‘Comparison of DMPC from experiments on real systems. Even if mistuned force simulations to experiment and a model-free method’’ (sim- fields or incomplete equilibration quantitatively distort the ulated results for F(q) and r(z) are compared to experiment, simulated structure, e.g., giving incorrect volume of water or targets are identified for CHARMM potential development, lipid molecules, the ensuing well-defined structure is still a and a simulation-based, model-free method for estimating A valid test case of the same generic type as typical bilayers. is proposed). Both force-field evaluation and the model-free The simulations used in this article were performed for method require consideration of the best statistical ensemble dimyristoylphosphatidycholine (DMPC) at five different for performing simulations. This is addressed in the Discus- fixed A, which bracket the previously determined value of sion and Conclusion section. It is argued that the constant 60.6 A ˚ 2 (6). Our primary test for models of r(z) is whether surface area rather than constant isotropic pressure ensem- they can determine A from the equivalent information avail- bles are more appropriate for the applications in this article able from x-ray experiments, which consists foremost of the because of the possibility of finite size effects and small defi- electronic scattering form factor F(q). Another important goal ciencies in the force field or methodology. of a structural model is to locate the component pieces of the lipid molecule within the bilayer and to determine the hydro- phobic thickness. The determination of A is a difficult test, METHODOLOGY one that neither HB nor H2 can pass, unless information in Molecular dynamics simulations were performed with the CHARMM addition to F(q) is provided to constrain the many parameters program (13) using the revised CHARMM27 (C27r) force field (14) and the required in any realistic electron density model. It should be modified TIP3P water model (15,16). The leapfrog Verlet algorithm was emphasized that this is not a criticism of the model method; used with tetragonal periodic boundary conditions and a time step of 1 fs. The Lennard-Jones interactions were smoothed by a switching function indeed, the advantage of the model method for r(z) is that ˚ (13). Constant particle number, pressure, surface area, and tem- over 8–10 A information from other experiments can be imposed on the perature ensemble (NPAT) simulations were run using the pressure-based model. This advantage is not offered by representing the elec- nonelectrostatic long-range correction (17) with a long-range cutoff of 30 A ˚. tron density profile by a Fourier series or Fourier transform. The particle mesh Ewald (18) method was used for the long-range (beyond 10 A˚ ) electrostatic contribution to the total energy with k ¼ 0.34 A ˚
1 and a In addition, structural models can also be extended to include fast-Fourier grid density of ;1 A ˚
1. All hydrogen atoms were constrained information from simulations, and this article develops guide- using the SHAKE algorithm (19). The extended system formalism was used lines regarding the kind of information that may be included. to maintain the temperature via the Hoover thermostat (20) with a thermostat The program that emerges from the preceding part of the coupling constant of 20,000 kcal mol
1 ps
2, and pressure was maintained introduction is to use simulations to produce and test a ge- with a barostat (21,22) with a piston mass of 2000 amu. neric model, which is then used to analyze the experimental The DMPC bilayer consisted of 36 lipids per monolayer (72 total) with 1848 water molecules with periodic boundary conditions in all directions with a data of Kucˇerka et al. (6). A second aspect of this article fixed A, i.e., the box lengths in the x and y direction are fixed. This system size involves the direct comparison of simulation and experi- has been shown to result in equivalent electron densities and other structural ment. By performing this comparison in q space, no model properties for systems larger than 72 lipids (23). Five trajectories with dif- is required, but the discrepancies are difficult to interpret. ferent cross-sectional areas (55, 59.7, 60.7, 61.7, and 65 A ˚ 2 per lipid) were Because the structural models represent the F(q) data very generated. The velocities were initialized at 203.15 K with a temperature increment of 10 K every 1 ps until the target temperature of 303.15 K was well, they can be used to carry out a comparison in real obtained, and the systems were then equilibrated for 3 ns. All averages were space. Nevertheless, the question arises, at which value of A evaluated for production runs of 10 ns with coordinates saved at 1 ps intervals. should one compare a simulation to experiment? Our answer The electron density profile r(z) along the bilayer normal z was obtained to this question leads to a simulation-based, model-free method as an average of the 10,000 snapshots following Feller et al. (24). To account for estimating the surface area. for temporal displacements of the entire bilayer along z, the center of the bilayer for each snapshot was taken to be zM, the mass weighted projection By way of outline, the following section describes the of the lipids along the z axis, and adjusted atomic positions zi were then methods used in the molecular dynamics simulations, and obtained from the raw atomic positions by subtracting zM. The small system discusses a common approximation related to the use of size suppressed undulations, so zM did not vary significantly with lateral atomic form factors to obtain F(q). The Results section is position in each snapshot. Based on the zi, the number of electrons in each divided into the following six subsections: i), ‘‘Component atom of lipid and water was then added to a histogram with a bin size of 0.1 A˚ in the z direction, Dz. Division by the bin volume and the number volumes from simulations’’ (an important first step for pro- of snapshots provided the electron density r(zj) for 660 values of zj, which viding constraints for the model); ii), ‘‘Simulated r(z) and includes water images to 633 A ˚. F(q)’’ (a comparison of these quantities from simulations of The continuous form factors, F(q), for symmetric bilayers (r(
z) ¼ r(z)), five areas, 55, 59.7, 60.7, 61.7, and 65 A ˚ 2); iii), ‘‘Structural are defined as models’’ (development of H2 guided by the simulation Z D=2 results and comparison with HB); iv), ‘‘Test of structural FðqÞ ¼ ½rðzÞ
rW cosðqzÞdz; (1) models’’ (fitting the simulated F(q) to determine what con-
D=2 Biophysical Journal 90(8) 2796–2807 2798 Klauda et al. where rW is the electron density of pure water. The discrete form factor from seven components listed above, as well as for a four-component model with simulation is determined at each value of qk, choline, phosphate, glycerol, and carbonyl combined into a single head- group distribution. In addition to these, we have obtained volumes for a six- Fðqk Þ ¼ +ðrðzj Þ
rW Þcosðqk zj ÞDz; (2) component model that combines the glycerol and carbonyl groups into one "j component, and a five-component model that additionally combines the where the electron density of water in simulations is 0.34 e A ˚
3 for TIP3P phosphate and choline into a single group. waters. The integrand in Eq. 1 at the upper and lower limits is zero because Deviations from a Gaussian distribution for the probability distributions r(z) is equal to rW , and similarly for Eq. 2 at r(zj) ffi rW. Values of F(qk) pm(z) of the component groups and combined distributions were quantified were obtained from the r(zj) for 800 values of qk evenly spaced from q ¼ 0 by kurtosis, g2 ¼ m4 =m2
3, and skew, g1 ¼ m3 =ðm2 Þ3=2 , where mi is the to qmax ¼ 0.8 A ˚
1; qmax is the upper experimental limit for DMPC (6). This ith sample moment about the mean. If a distribution is Gaussian, then g1 and procedure (11) assumes that the electrons are localized at the atomic nucleus, g2 are equal to zero. which is equivalent to assuming that the atomic form factors fi(q) are con- stants equal to fi(0). Benz et al. (12) have recently emphasized that the atomic form factors are RESULTS not constants so that one should calculate AFðqÞ ¼ +i2A fi ðqÞcosðqzi Þ which is only the same as the preceding procedure when fi(q) are constants. Component volumes from simulations However, Fig. 1 shows that the relevant fi(q)/fi(0) (25) deviate by only ;2% Spatial distributions of the component groups are shown by from 1.0 at the upper experimental range of q-values in reference (11). The deviation for our upper experimental range is only 5% because the dis- pm(z) in Fig. 2. The average deviations from unity of the sum tribution of electrons around nuclei is highly concentrated within a radial of all the component probabilities are of the order of 61%. distance of order sel ; 0.3 A ˚ . This distribution would require a spatial The region with the largest deviations occurred near the convolution of the electron density in the z direction, but only over the bilayer center where the average deviations were 62.4%. It distance sel which is typically five times smaller than the van der Waals may also be noted that, although the values of the volumes radii of atoms. This correction makes little difference to r(z) or F(q) in the experimental range of q, because the intrinsic disorder in the bilayer already Vm modulate the maximal values of the individual pm(zi) in broadens the distribution functions for the locations of the nuclei by sin . 2 Fig. 2, the locations of the maxima (which locate the mean A˚ and the total broadening s ¼ (s2in 1s2el Þ1=2 is negligibly different from the positions of the component groups along z) are independent broadening sin of the nuclei alone. Indeed, the A and r(z) obtained using the of the volumetric analysis. atomic form factor correction to F(q) were nearly identical (within 0.1%) to Table 1 lists the lipid component volumes for the five values obtained without the correction. It should be noted that the use of atomic form factors is only exact for atoms. Because lipids are molecules, surface areas simulated using a six-component volumetric their valence electrons are displaced from atomic orbitals. Consequently, the analysis. Four-, five-, and seven-component analyses were use of atomic form factors is not exact. One needs molecular orbitals and also performed. Standard deviations obtained by comparing orientation dependence of chemical bonds. However, this complication the four volumetric analyses were ;1.0 A ˚ 3 (60.09%) in the makes as little difference to F(q) as the use or nonuse of atomic form factors total volume VL. Consistent with reference (27), standard described above. Electron density weighted histograms rm ðzj Þ were obtained for each of deviations in the sum of the volumes of the phosphate and the m ¼ 1,. . .,7 groups: water, choline, phosphate, glycerol, carbonyl, choline were smaller than the deviations in the individual methylenes on the tails, and the terminal methyls on the tails. Following the components. The total headgroup volume VH is nearly inde- method of Petrache et al. (26), these rm(zj) were converted into probability pendent of simulated area A. The constancy of VH is expected distributions, i.e., pm ðzÞ ¼ rm ðzÞVm =nm . The sum of all probabilities, because the headgroup is largely immersed in water. This pT ðzj Þ ¼ +m pm ðzj Þ, should ideally be unity for each zj bin, and this method obtains the component volumes Vm by minimizing +j ðpT ðzj Þ
1Þ2 . The simulation result supports the assumption used in structural deviations from unity test the assumption that the component volumes are modeling that the value of VH determined experimentally for independent of z. Petrache et al. (26) obtained component volumes for the the gel phase can be used for determination of the fluid phase structure. The total chain volume VC shows a small sys- tematic increase as A is increased; this is consistent with more disordered chains requiring greater volume. The water volume VW is independent of A with the volume of water essentially equal to that in the bulk. Simulated r(z) and F(q) Fig. 3 shows total electron density profiles r(z) for three of the simulated areas. The simulation at A ¼ 65 A ˚ 2 is nearly symmetric and fairly smooth, which is consistent with this simulation having reached equilibrium. As the simulated area is reduced, the simulated electron densities become less smooth and more asymmetric, suggesting that equilibration FIGURE 1 Normalized atomic form factors fi(q)/fi(0) for carbon, oxygen, takes longer, possibly due to stronger excluded volume con- phosphorus, and nitrogen atoms within the experimental q-range (0 , q , straints in the headgroup region. However, the ‘‘high- ˚
1). 0.8 A frequency’’ roughness of these electron density profiles has a Biophysical Journal 90(8) 2796–2807 Interpreting X-Ray Data from Bilayers 2799 FIGURE 2 The bottom panel shows the component ˚ 2 simulation, pm(z) along the probabilities for the A ¼ 60.7 A bilayer normal z for water (w), choline (chol), phosphate (phos), glycerol (gly), and carbonyls (co) on the left, and on the right for combinations of some of these components, phosphate 1 choline (PC), carbonyl 1 glycerol (CG), and water 1 choline, with chain methylenes (CH2) and termi- nal methyls (CH3) and their sum in the middle. The Gibbs dividing surfaces are indicated by vertical dashed lines labeled DC for the hydrocarbon boundary and 0.5DB for the Luzzati water boundary. The top panel shows devia- tions of ptot(zi) from unity with the right half from the six- component analysis and the left half from the seven- component analysis. negligible effect on the calculated form factors F(q) within some of the lipid component groups in the bilayer (4,8). the experimental range 0 , q , 0.8 A ˚
1. Fig. 4 shows the Nevertheless, the distribution functions for any component corresponding F(q). The F(q) curves vary significantly, group need not be purely Gaussian and indeed, deviations which demonstrates that experimental measurements of F(q) were observed in earlier simulations (10). A comparison is should be important for determining A and bilayer structure. shown in Fig. 5 for various lipid components at 60.7 A ˚ 2. The values of kurtosis g2 in the distributions for choline, pho- sphate, glycerol, and carbonyl for the simulation at 60.7 A ˚2 Structural models are small
0.07,
0.19,
0.14, and
0.06, respectively, and A major issue in structural modeling is the number of adjust- similar small values are calculated for other A. In general, the able parameters. It is desirable that a model be able to rep- distributions are more Gaussian for the phosphate 1 choline resent all interesting features of lipid bilayers. On the other (PC) and carbonyl 1 glycerol (CG) combined components, hand, a model with too many parameters can fit the data by with g2 ¼ 10.06 and
0.07, respectively. Similarly, g1 in different combinations of the parametric values; i.e., the the individual group distributions is reduced from about parameters are underdetermined. In general, simple func-
0.2 to 10.1 when the headgroups are combined. In con- tional forms with few parameters that still provide a good trast, a substantially larger kurtosis, g2 ¼ 11.03, is obtained representation of the data and physical features are preferable for the distribution of methyls from both monolayers (Fig. 5, to more general forms with more parameters. Here a new bottom panel). The skew is zero to within statistical error by structural model is developed with a robust number of param- symmetry. The distribution of terminal methyls from only eters based upon our simulation results. one monolayer (not shown) is strongly skewed toward the The most realistic structural models currently use the headgroups of that monolayer (g1 ¼ 0.3). Gaussian functional form to represent the distributions of The new structural model, H2, consists of functional forms that provide excellent representations of the electron TABLE 1 Volumetric results from simulations using the densities for the five components shown in Fig. 5, six-parameter volumetric analysis for total lipid volume (VL) and component volumes for water (VW), chain methylene (VCH2), terminal methyl (VCH3), phosphate (Vphos), choline (Vchol), carbonyl 1 glycerol (VCG), total head (VH), total chains (VC), and r ¼ VCH3 / VCH2 Simulated Experiment ˚ ) 2 A (A 55 59.7 60.7 61.7 65 60.6 VL (A ˚ 3) 1061.3 1072.0 1072.3 1070.4 1074.6 1101 VW 29.5 29.5 29.5 29.5 29.5 30.0 Vchol 109.4 109.7 105.5 109.7 108.1 – Vphos 69.7 68.5 72.85 68.0 69.2 – VCG 142.6 145.3 145.6 145.6 147.3 – VCH2 26.3 26.7 26.8 26.7 26.9 27.7 VCH3 54.0 53.7 53.0 52.9 52.8 52.6 VH 321.6 323.5 323.9 323.4 324.5 331 VC 739.7 748.5 748.4 747.1 750.1 770 FIGURE 3 The electron density profiles, r(z), as a function of z along the r 2.05 2.01 1.98 1.98 1.97 1.9 bilayer normal for simulated areas 55 (solid gray), 60.7 (solid black), and 65 Experimental column from Kucˇerka et al. (6). ˚ 2 (dashed black). A Biophysical Journal 90(8) 2796–2807 2800 Klauda et al. One Gaussian represents the contribution of the phosphate group in the upper leaflet to the electron density profile, pffiffiffiffiffiffi
1=2  2 2  GP ðz; zP ; sP Þ ¼ CP ðsP 2pÞ exp
ðz
zP Þ =2sP ; (4a) and a similar GP Gaussian with parameter
zP represents the phosphate in the lower leaflet, so rP ðzÞ ¼ GP ðz; zP ; sP Þ 1 GP ðz;
zP ; sP Þ: (4b) A single Gaussian models the terminal methyls from both leaflets, rCH3 ðzÞ ¼ GCH3 ðz; 0; sM Þ FIGURE 4 Form factors, F(q), from three of the five simulated areas, 55 pffiffiffiffiffiffi
1=2  2 2  ˚ 2 (dashed black). (solid gray), 60.7 (solid black), and 65 A ¼ CCH3 ðsCH3 2pÞ exp
z =2sCH3 : (5) By symmetry, combining the methyl distribution from both H2 r ðzÞ ¼ rP ðzÞ 1 rCH3 ðzÞ 1 rCG ðzÞ 1 rCH2 ðzÞ 1 rBC ðzÞ; (3) leaflets results in a skew of zero, though there remains a substantial positive kurtosis (Fig. 5). A second Gaussian for where the notation for the densities is rP(z) for the phosphate the methyl density does not significantly improve the overall groups, rCH3(z) for the terminal methyls, rCG(z) for the fit of the model to F(q), and is not included in rCH3(z) to carbonyl 1 glycerol, rCH2 ðzÞ for the methylenes on the hydro- avoid additional adjustable parameters. carbon chains, and rBC(z) for the water 1 choline (BC). The H2 uses just one Gaussian for the carbonyl and glycerol functional forms are described next. groups in each leaflet, rCG ðzÞ ¼ GCG ðz; zCG ; sCG Þ 1 GCG ðz;
zCG ; sCG Þ: (6) Combining the carbonyl and glycerol groups in each mono- layer in a single Gaussian reduces the number of parameters. This simplification arises because the distributions of these two groups overlap considerably (Fig. 2). Each Gaussian has parameters for its width s, and integrated size, C. GP and GCG also each have a parameter for the position z along the bilayer normal; GCH3 is constrained by symmetry to zCH3 ¼ 0. There are a total of eight parameters for the first three terms on the right side of Eq. 3. However, the number of electrons, nei , is known for each component group and equals the molecular area A multiplied by the integral of the Gaussian over z. Therefore, CP 3 A ¼ 47 for each of the two phosphate group Gaussians, CCH3 3 A ¼ 36 for the single methyl Gaussian, and CCG 3 A ¼ 67 for each of the two carbonyl-glycerol Gaussians. These physical constraints re- duce the number of independent Gaussian parameters to five. Although Gaussians provide good approximations for small, localized groups, they are clearly inappropriate for representing the many methylene groups on the hydrocarbon tails of the lipids (Fig. 5). As illustrated by the pCH3 1 pCH2 curve in Fig. 2, these methylenes and the terminal methyls together comprise the entire hydrophobic core of the bilayer. The composite probability distribution is well represented by the sum of two classical error functions (also used in Schalke et al. (9) to model monolayers) FIGURE 5 Results of independently fitting Gaussians to the phosphate, terminal methyl, and CG distributions and the other functional forms in H2 pHC ðzÞ ¼ 0:5½erfðz;
DC ; sCH2 Þ
erfðz; 1 DC ; sCH2 Þ; to the water 1 choline and methylene distributions for the A ¼ 60.7 A ˚2 (7a) simulation. The solid lines are results from MD and the dashed H2. The bottom panel shows the terminal methyl distribution on an expanded scale. where the error function (erf) is defined by Biophysical Journal 90(8) 2796–2807 Interpreting X-Ray Data from Bilayers 2801 Z pz
m ffi2s The total electron density profile in H2 is obtained by 2 2 erfðz; m; sÞ ¼ pffiffiffiffi exp½
x dx; (7b) summing the components of Eq. 3, specified in Eqs. 4–9. p 0 The first column of Table 2 shows how the separate com- with location m and width s. One parameter is required in ponents yield a total of 16 parameters. Thus far, a total of five H2 for the average locations DC and
DC of the boundaries constraints have been noted for the number of electrons, re- of the hydrocarbon interfaces, otherwise identified as the ducing the number of independent parameters to 11. Gibbs dividing surfaces between the hydrocarbon region Now the next type of constraints that involve volumes is and the headgroup region. Another parameter sCH2 gives introduced. Experimentally, the volume VL of the lipid mol- the widths of these surfaces (68% of the change from total ecule is the most accurate datum. This allows the elimination hydrocarbon to no hydrocarbon occurs within DC 6 sCH2 ). of nW as a free parameter because the volume AD/2 of half However, to obtain the contribution of just the methylenes the experimental or simulation unit cell is just VL 1 nW VW. to the electron density, it is necessary to subtract the ter- As shown by Nagle and Wiener (28), this constraint is equi- minal methyl distribution from pHC, taking into account valent to the relation that the number of electrons (neM ¼ 9) and the volume VCH3 e of the terminal methyls are different from the methylenes AFð0Þ ¼ 2ðnL
VL rW Þ; (10) (28). For the contribution of the methylenes to the electron where neL is the number of electrons in the lipid molecule and density profile, F(0) is the integral of (r(z)
rW). Because H2 combines the water 1 choline into a single distribution, the total volume rCH2 ðzÞ ¼ CCH2 pHC ðz; DC ; sCH2 Þ
ð8r=9ÞGCH3 ðz; 0; sCH3 Þ; constraint used in H2 is (8) e nchol where the parameter CCH2 is the electron density of the ADBC ¼ VL
; (11) methylene region. CCH2 is proportional to 8/VCH2, and rW the parameter defined by r ¼ VCH3/VCH2 is employed in the which is derived under the assumption that the electron den- terminal methyl subtraction. Furthermore, the integral of sity of choline region equals that of water. rCH2 ðzÞ3A should be constrained to be the total number of The simplest and most powerful volumetric relation for chain methylene electrons (192 for DMPC). This constraint H2 is reduces the number of independent parameters required for the methylenes from four to three. A ¼ VC =DC ; (12) The final term in Eq. 3 is the water 1 choline distribution. The water distribution shown in Fig. 2 is not well described TABLE 2 Parameter count for H2 and HB models by a simple form. One error function does not provide a good fit and two error functions proliferate the number of param- H2 HB eters. However, Fig. 2 suggests that the composite distribu- P C I P C I tion function consisting of water 1 choline component of the e P or PC head 3 n 2 3 R 2 headgroup can be well represented by error functions CG head 3 ne 2 3 VL 2 CH3 2 ne 1 2 r 1 rBC ðzÞ ¼ rW ½1
0:5ðerfðz;
DBC ; sBC Þ CH2 4 ne ,sCH2,r 1 1 – 1 BC or water 3 rW, VL,sBC 0 1 rW 0
erfðz; 1 DBC ; sBC ÞÞ: (9) Baseline function – – – 2 wb, zb 0 Area 1 VC 0 1 VC 0 H2 exploits this by using an electron density contribution, DH1 – DH1
1 – DH1
1 rBC(z), consisting of two parameters (DBC and sBC), mul- Totals 16 11 5 13 8 5 tiplied by the known electron density of pure water, rW. If the feature represented by rBC(z) corresponded only to Constraints ne No. of electrons water, then DBC would be the Luzzati thickness defined as r Ratio of methyl to methylene volume DB and shown in Fig. 2. This is not the case because the DBC rW Known water electron density in Eq. 9 includes the choline component. The integral of VL Total lipid volume (only lipid) rBC(z) 3 A is the total number of electrons of choline plus VC Chain volume (lipid volume minus headgroup) the number of electrons corresponding to nW water mole- sCH2 Width of methylene error function sBC Width of BC error function cules per lipid in the simulation cell. However, this rel- R Ratio of headgroup peak areas ationship does not immediately reduce the number of wb Width of bridge in baseline function independent parameters because nW is a parameter that zb Position of bridge in baseline function cannot be measured experimentally for fully hydrated sam- P is the number of parameters for each feature, column C abbreviates the ples (3). Therefore, nW should not be taken from the simu- names of the constraints, and I is the number of independent degrees of lations for the purpose of testing models. freedom in the fitting. Biophysical Journal 90(8) 2796–2807 2802 Klauda et al. which immediately yields A from the fitted DC in Eq. 8 and the electron density in these components in excess of from VC, which is obtained by subtracting the headgroup the baseline function. The constraints r, VC, and DH1 are the volume VH (29) from the total lipid volume VL. same as those applied to H2, as indicated in Table 2. The Experiment obtains estimates for the volumetric r ratio constraints in H2 on ne for the component groups have two (4,28), so this is constrained in H2. These three volumetric counterparts in HB. The first is the R constraint on the ratio constraints (VL, VC, and r) therefore reduce the number of of the integrated sizes of the two Gaussian headgroup peaks independent parameters in H2 from 11 to 8. Three additional and the second is a VL constraint in Eq. 10. Table 2 lists the constraints are needed to maintain robustness in the model total number of independent parameters as five when A is fits. The widths of the error functions of the BC and CH2 counted as a parameter and the bridge in the baseline func- distributions, si, were too flexible in the unconstrained fits. tion is constrained as described above (6). Therefore, these values were constrained to within 60.1A ˚ Although the baseline function reduces the number of from the simulated value by soft Bayesian constraints. A final parameters in HB, its primary description does not include a constraint for H2 refers to the distance DH1 obtained from gel most important feature, namely, the hydrophobic boundary phase studies; DH1 is the distance between the location DHH/ DC that is included explicitly in H2. Therefore, the A cannot 2 of the maximum in the electron density and the location DC be directly determined from Eq. 12 for HB. Instead, DC for of the hydrocarbon Gibbs dividing surface. The use of these HB is obtained from the headgroup peak location DHH/2, constraints reduces the number of independent parameters to using DC ¼ DHH/2
DH1, where DH1 is obtained from the five. gel phase (29); this is equivalent to the bootstrap method of HB has been amply described in previous applications McIntosh and Simon (30). (6,29), so the focus is on the differences with H2. HB consists of four terms, Test of structural models HB r ðzÞ ¼ rb ðzÞ 1 rCH3 ðzÞ 1 rPC ðzÞ 1 rCG ðzÞ; (13) As a first test, H2 was fit to the simulated F(q) without which are the electron densities for the baseline, rb(z), constraining DH1 or the widths of the error functions. Only methyl, rCH3(z), phosphate 1 choline, rPC(z), and carbonyl form factors at q , 0.8 were used in all fits because that is the 1 glycerol, rCG(z). The major difference is that HB reduces experimentally accessible range. The A obtained from this the number of model parameters with a baseline function eight-parameter fit deviated significantly from the actual rb(z) to represent both the methylenes and the water. This simulated surface area, e.g., for the simulation at 60.7 A˚ 2 the ˚ 2 baseline function employs a smoothly varying bridge between predicted area was 66.0 A . In addition, the two fitted Gibbs the known electron density of bulk water, rW, and a meth- dividing surfaces, DBC and DC, were .1 A ˚ too close to the ylene plateau, rCH2 . The bridge has two independent param- bilayer center and their widths were too large. Constraining eters, one for the location of the center zb of the bridge and the si-values, but not DH1, improved the value of A, but only one for its width wb. For gel phases the difference in the to 64.0 A ˚ 2. Despite the disagreement with A, both of these electron densities of the methylene plateau and water is small fits provided excellent agreement with F(q) and the total r(z). (;5%), so structure determination is rather insensitive to the This demonstrates that the unconstrained H2 with eight or bridge parameters. The difference is larger for fluid phases six fitted parameters is underdetermined and parameter flex- (;20%), but simulations have enabled the location of the ibility results in poor component determination. Consequently, bridge to be constrained relative to the headgroup peaks (6). all the constraints for H2 listed in Table 2 are required. The width of the bridge has also been constrained to be the Even with only five parameters, the fits of the model to the width of the region that simulated r(z) contains both hydro- simulated F(q) data have such small deviations that they carbon and water, ;8 A ˚ as seen in Fig. 2. These constraints cannot be distinguished graphically from the simulated data play a similar role as the sCH2 and sBC constraints used in shown in Fig. 4. This indicates that the model is more than H2, but are considerably different in detail. H2 requires four adequate to account for primary F(q) data from x-ray dif- parameters, DC, sCH2 , DBC, and sBC to describe the baseline fraction. The predicted A obtained from H2 with all con- features that are incorporated by only two parameters for the straints are compared to the simulated A in Table 3. There is bridge in the HB baseline function. The fit to the F(q) data is excellent agreement with the simulated A, where the H2 more sensitive to the H2 parameters DC and DBC, which vary determined A has negligible bias and an average root mean strongly depending upon the lipid. square deviation of 0.1 A ˚ 2. In HB the headgroup and terminal methyls are also rep- TABLE 3 Area A in A ˚ 2 for the H2 and HB structural models resented by Gaussians with the difference that they are when fit to F(q) obtained from simulations performed at the superimposed on the baseline function. Therefore the methyl exact A shown in the top row Gaussian is a negative trough (like the last term in Eq. 8 for Simulation 55.0 59.7 60.7 61.7 65 RMSD Bias H2), and represents the deficit in electron density compared to the more electron dense methylenes. Similarly, the head- H2 54.8 59.9 61.0 61.7 64.6 0.1 10.02 HB 54.7 60.0 60.6 61.8 65.3 0.1 10.06 group Gaussians, GPC and GCG, are scaled to represent only Biophysical Journal 90(8) 2796–2807 Interpreting X-Ray Data from Bilayers 2803 Fig. 6 shows that the H2 parameter fit yields good repre- TABLE 4 Values of the H2 parameters fit to the simulated sentations of the electron densities of the individual compo- and the experimental form factors nents, although there are small, but noticeable, deviations in Asim ¼ 60.7 A ˚2 Fit to x-ray the carbonyl 1 glycerol and phosphate peaks. The methyl rCH3 1.98* 1.9* trough tends to be slightly higher than simulations because CP 0.77 0.78 kurtosis is absent in GCH3. However, the model results agree zP 17.72 17.83 well overall with the simulated r(z). The parameters for the sP 2.37 2.18 ˚ 2 are CCG 1.10 1.11 constrained H2 fit to the simulated F(q) with A ¼ 60.7 A zCG 13.64 13.88 listed in Table 4. sCG 2.32 2.28 HB was fit with five independent parameters and the CBC 0.34 0.33 constraints listed in Table 2. It fits the simulated F(q) in the DBC 15.17 15.68 experimental range 0 , q , 0.8 so well that, like H2, one sBC 2.77y 2.96y CCH2 0.30 0.29 cannot discern any deviations from the simulated F(q) curves DC 12.27 12.70 on the scale of Fig. 4. The bottom panel of Fig. 6 shows that sCH2 2.31y 2.32y rHB(z) from the fitted model agrees well with the simulated CCH3 0.59 0.59 r(z). Fig. 6 also shows the individual terms of HB and allows sCH3 2.89 2.22 comparison with the simulated contributions from the molec- A (VC/DC) 61.0 60.6 DH1 5.28* 4.95* ular components. GPC is located very close to the phosphate ˚ 2/0.1A @A/@r (A ˚) 0.32 0.19 distribution and GCG is located near the carbonyl distribu- @A/@DH1 (A ˚ 2/0.1A ˚) 0.53 0.45 tion. The size of the GCG is considerably smaller than the RMSD 0.013 0.022 sum of the carbonyl and glycerol contributions because the *Hard constrained values. baseline function contains a fraction of the carbonyl and y Soft constrained values. glycerol electrons. GPC is larger than GCG by the constrained factor R ¼ 1.76 because the electron density of the phosphate Application of the models to experimental is much larger and a smaller proportion of its electrons are x-ray data included in the baseline function. The results for A are listed in Table 3. The HB model has previously been applied to DMPC exper- imental F(q) and volumetric data (6). This section applies H2 to the same data. The constrained parameters VL, r, and DH1 were set to values obtained from experiment rather than the simulated values shown in Tables 1 and 4. The experimental uncertainties for r are estimated to be of order 60.1 and for DH1 of the order of 60.1 A ˚ . Table 4 examines the sensitivity of the fitting results on these parameters, i.e., @A=@r and @A=@DH1 . Clearly, the value of A depends strongly on DH1 with a change in A of 0.45 A ˚ 2 for every 0.1 A ˚ change in DH1. The H2 fits are less sensitive to a change in r, where @A=@r ¼ 0.19 A ˚ 2 per 0.1 A ˚ change in r. The resulting A for these H2 fits is 60.6 6 0.5 A ˚ 2 with a confidence based on the uncer- tainties in r and DH1. The fits to the experimental F(q) data are very good as shown in Fig. 7 and Table 4 for H2 and by Kucˇerka et al. (6) for HB. Because the fits to the simulated F(q) have negligible RMSD, the H2 RMSD in Table 4 contains mostly exper- imental error in F(q). The fits to the F(q) are equally good when r and DH1 are varied within their estimated uncertainty. Therefore, the accuracy of the values used in these con- straints cannot be deduced from the F(q) data and their un- certainties propagate uncertainty in the determination of A. However, the model form factors for different values of the FIGURE 6 Results of fitting H2 (top panel in red) and HB (bottom panel constraints begin to differ for q-values that exceed the cur- in blue) for A ¼ 60.7 A ˚ 2 with the total r(z) and the component r(z) rent experimental range, as seen in Fig. 7; this emphasizes (simulation results in black). The CG, PC, and CH3 component contribu- tions for the HB model are shown as differences from the water level rW and the desirability for obtaining data to the highest possible the total electron density is the sum of the baseline and the component q-value. None of the preceding model fits change the sign of contributions. F(q) near q ¼ 0.7, and the locations of the maxima in the Biophysical Journal 90(8) 2796–2807 2804 Klauda et al. FIGURE 7 H2 form factors fit to the experimental F(q). The red H2 curve FIGURE 9 A comparison of the experimental form factors (6) with those shows the result for the parameters in Table 4 and the other two H2 plots are from simulations at two areas. The experimental F(q) was scaled to MD 60.7 A˚ 2 and MD 61.7 A ˚ 2 was artificially rescaled to the experimental F(q) to for the altered values of DH1 and r given in the legend. better view the residuals for that simulation. lobes and the crossing points where F(q) ¼ 0 are nearly the unilamellar samples and one for the oriented samples. If identical. it is assumed that these two relative scaling factors were The r(z) of the HB and H2 models are compared in Fig. 8. obtained correctly, then this permits only one scaling factor Overall agreement is satisfactory, although there are distinct to compare to simulations; this gives the total root mean differences in the electron densities for various positions square deviation (RMSD) listed in the ‘‘One factor’’ column within the bilayer. HB has a higher phosphate peak than H2 of Table 5. If two separate scaling factors for each sample and a lower carbonyl-glycerol shoulder. Kucˇerka et al. (6) type are employed, the RMSD in the last column of Table 5 reported A ¼ 60.6 6 0.5 A ˚ 2 using the HB model, in agree- is obtained. Only a small decrease in RMSD is obtained by ment with H2. employing both scaling factors, which is consistent with the relative scaling factor having been chosen correctly by Kucˇerka et al. (6). Comparison of DMPC simulations to experiment The results in Table 5 show that the simulations fit the and a model-free method experimental data best for A between 60.7 and 61.7 A ˚2 The simulated form factors are compared with experiment in within a standard error of the model-based value (60.6 6 0.5 Fig. 9 for two simulated surface areas, and the deviations A˚ 2). Assuming that the RMSD is parabolic with respect to A, from experiment are listed in Table 5. The absolute scale of the minimum RMSD occurs at 61.1 A ˚ 2. This simulation- the experimental F(q) is unknown and simulations can help based estimate for the area from the experimental data is to obtain it (11). There were two scaling factors embedded in independent of the structural models and is referred to here the experimental F(q) data from Kucˇerka et al. (6), one for as the model-free method. The simulated F(q) in Fig. 9 cross zero for a short range of q-values near q ¼ 0.7 where the experimental F(q) are very small. In contrast, neither H2 nor HB cross zero in that TABLE 5 Comparison of experimental (6) and simulated F (q) F(q) scaling ˚ ) 2 A (A One factor Two factors 55.0 0.22 (0.007) 0.19 (0.009) 59.7 0.072 (0.006) 0.066 (0.006) 60.7 0.044 (0.003) 0.042 (0.003) 61.7 0.047 (0.002) 0.047 (0.002) 65.0 0.12 (0.001) 0.12 (0.002) The RMSD was obtained from the difference of the F(q) from simulations at different areas and the experimental F(q). The standard error of the RMSD for the simulations, given in parentheses, were calculated from 2.5- ns blocks. The experimental F(q) were scaled to best fit the simulated F(q) FIGURE 8 Comparison of the r(z) obtained from HB and H2 (from Table 4) either with a single factor for both experimental samples (labeled as ‘‘One fit to experimental form factors. factor’’) or with ‘‘Two factors’’, one for each sample type. Biophysical Journal 90(8) 2796–2807 Interpreting X-Ray Data from Bilayers 2805 region. Although this is a clear difference, the experimental x-ray F(q), to provide better values of structural parameters data alone do not afford a clear indication that the simu- for lipid bilayers. The results of DMPC simulations reported lations are incorrect. The more important comparison is that here have guided the development of a new structural model, the root mean square residuals for the entire q range are larger H2, which includes additional structural features in a more for the simulations (0.042) than for the models (0.022). transparent way than the previously employed model, HB. The F(q) do not lend substantial insight into the origin of The tests with simulations were designed to mimic the way the differences between simulation and experiment and to experimental data are analyzed, with a nonlinear least squares where one might look to improve the simulation. For this, fitting to the F(q) data constrained by additional data, such as the electron densities of the simulation and experiment are the volume of the lipid, and outside information, such as the compared in Fig. 10. The small differences between HB number of electrons in component groups. and H2 due to different functional forms are averaged as a Because H2 includes the hydrocarbon thickness DC composite experimental result. The comparison of simula- explicitly, in principle it is not necessary to use the DH1 con- tion and experiment in Fig. 10 illustrates three regions with straint obtained from the gel phase. Such a feature would differences. The first is the water region, where the simulated provide a substantial advantage to H2 over HB. In practice, electron density for the bulk water region is higher than real however, without constraint DH1 H2 does not obtain satis- water due to the known inaccuracies of the TIP3P water factory values of area A as shown by fits to simulated data. model (31). This also is evident in the lower water volumes H2 obtains accurate values of A with the DH1 constraint, in Table 1. At 303 K and 1 bar it was found that the density of accurately fits the F(q) data in the experimental q range, and water with TIP3P is 1.6% higher than experiment and is the reproduces the simulated total and component r(z) (Fig. 6). cause for the higher electron density away from the bilayer HB was tested on the simulated data and it performed about center. The second region of discrepancy is the higher and as well as H2. Although both models have many parameters more prominent shoulder on the headgroup peak at z ¼ 14 A ˚ to provide realistic representations of bilayer structure, both near the location of the CG group. Third, the simulated have five independent degrees of freedom when the neces- methyl density at the bilayer center and the chain methylene sary number of constraints are applied. plateau density near z ¼ 8 A ˚ is consistently higher than the Having passed the simulation test, H2 was applied to experimental results. This discrepancy is consistent with the experimental data for DMPC (Fig. 7). The overall model under prediction of the chain volume VC in Table 1. Fig. 10 results were in excellent agreement with the earlier results indicates that these differences in the overall electron density obtained with the HB model. The predicted surface area per profile are fairly minor; in particular, it is encouraging that lipid for both models is 60.6 6 0.5A ˚ . Although Fig. 8 shows the locations of the headgroup component distributions are small differences in r(z) in the carbonyl 1 glycerol shoulder nearly identical for experiment and simulation. and the height, though not the position, of the maximum, there is near-perfect agreement for the other regions. It appears that neither structural model is clearly superior to the other. DISCUSSION AND CONCLUSIONS We believe that both models, with their rather different func- The primary goal of this article is to use simulations to tional forms, are valuable because their combined use pro- improve modeling of experimental structural data, especially vides an estimate of uncertainties in r(z). An equally important purpose of this article is to demon- strate how to use experimental data to improve simulations. As has been emphasized recently (6,11,12), the primary comparison of simulations to experimental diffraction data should be between the F(q) obtained from simulations and the experiment because this is a direct test that does not involve structural modeling. The first significant result of this test was that the simulations agree fairly well with the experimental F(q) when the simulated area was close to the value obtained by modeling (Fig. 9). Indeed, finding the A that best fits the experimental F(q) is a model-free simula- tion-based method for obtaining A. If the potentials used in simulation are accurate, then this model-free method will be applicable to other systems such as lipid mixtures and bilay- ers with incorporated peptides. This method would be superior to structural modeling, because it avoids the need for more FIGURE 10 The r(z) obtained from the structural models fit to the structural model parameters than can be successfully fitted experimental F(q); average of H2 and HB (blue) and components of H2 to the available experimental data. The model-free method ˚ 2. (red). The black curves show the simulations for A ¼ 60.7 A for DMPC results in A ¼ 61.1 A ˚ 2, which is within the Biophysical Journal 90(8) 2796–2807 2806 Klauda et al. confidence of the structural modeling value, A ¼ 60.6 6 priori, it is still prudent to consider that shortcomings in the 0.5A ˚ 2. This suggests that the current potentials are already simulation potentials could contribute to a nonzero value of reliable for many purposes. Nevertheless, the statistically sig- g. For example, the surface tension for pure liquids such as nificant differences between the simulated and experimen- water is highly sensitive to the potentials, their cutoffs, and tal F(q) (Table 5 and Fig. 9) do indicate deficiencies in the the lack of polarizability (40,41). Such shortcomings would CHARMM potential. also distort the surface tension of bilayers but would not Because discrepancies in reciprocal space are difficult to necessarily have a large effect on the structure, provided that interpret for improving real space potentials, the use of real it is simulated at the correct value of A. In contrast, sim- space modeling of the experimental data provides a more ulations in the NPT ensemble allow these small differences insightful comparison to the simulations. As shown in the pre- in surface tension to distort A (42) and thereby produce poor vious subsection, the comparison highlights the well-known agreement with F(q). While obtaining agreement of exper- deficiency in the density of TIP3P water (31), and confirms iment and large simulation systems constrained to g ¼ 0 is an the volumetric analysis of the simulations that the hydro- ultimate goal, simulating only at NPT appears unduly restric- carbon chain volume is smaller than experiment (Table 1). A tive and limiting because it magnifies small flaws in the new but small discrepancy is also observed in the carbonyl potentials. This point has also been convincingly made by region of the headgroup, and the methyl trough is insuffi- Ane´zo et al. (42), which emphasizes that obtaining the cor- ciently deep. These all provide clear targets for ongoing rect area per lipid is a poor measure of the force-field quality development of CHARMM potentials. or methodology. The preceding discussion pivots on what ensemble and In conclusion, we propose that simulations be performed values of thermodynamic parameters should be employed in at several areas in the NPAT ensemble (or at several surface simulations. There are two distinct approaches. The approach tensions in the NPgT ensemble) as part of a broad-based employed here is to carry out simulations in the NPAT ensem- analysis of a bilayer or biomembrane. The best fit to exper- ble at or near the experimentally derived surface area. A imental data provides a simulation-based model-free value parameter set is considered well tuned if simulated and ex- for A. When the model F(q) fits the experimental data as well perimental properties agree at this surface area. Equivalent as it does for DMPC in this article, and the electron density results would be expected from simulations carried out in the profiles from different models agree, then comparison of the NPgT ensemble, where g is the surface tension evaluated at simulated electron density with the models provides insight the experimental surface area (32–34). However, there is one both into deficiencies of the simulation and into the structure property that is poorly obtained in this approach. The bilayer of the bilayer. surface tension, g, which is identically zero experimentally We acknowledge Horia I. Petrache, Jonathan N. Sachs, and Douglas J. for flaccid bilayers (35), is 19.8 6 2.9 dyn/cm/monolayer for ˚ 2. Finite Tobias for their helpful comments on this manuscript. the present DMPC system at our best A ¼ 60.7 A size effects have been proposed (36) as the reason why the This research was supported in part by the Intramural Research Program of the National Institutes of Health (National Heart, Lung and Blood Institute) surface tension should differ from zero in simulations, even and by the Extramural Program of the Institute of General Medicine, grant if the potentials were perfectly tuned. Subsequent theoretical GM44976 (J.F.N. and N.K.). work (37) supports this notion, though it leads to somewhat smaller surface tensions than presently obtained in CHARMM- based simulations. System size dependence of the area has REFERENCES been observed in some recent simulation studies (38) but not 1. Janiak, M. J., D. M. Small, and G. G. Shipley. 1979. Temperature and in others (39). Another approach is to carry out simulations compositional dependence of the structure of hydrated dimyristoyl in the constant isotropic pressure ensemble (NPT). This is lecithin. J. Biol. Chem. 254:6068–6078. equivalent to the NPgT ensemble with imposition of the 2. McIntosh, T. J., and S. A. Simon. 1986. Area per molecule and distribution of water in fully hydrated dilauroylphosphatidylethanol- requirement that g ¼ 0. Under most conditions the area of amine bilayers. Biochemistry. 25:4948–4952. bilayers contracts and the bilayer becomes correspondingly 3. Nagle, J. F., and S. Tristram-Nagle. 2000. Structure of lipid bilayers. thicker when simulated at NPT with the present CHARMM Biochim. Biophys. Acta. 1469:159–195. potentials (12,33). This thickening has a strong effect on 4. Wiener, M. C., and S. H. White. 1991. Fluid bilayer structure F(q). The agreement of the simulated and the experimental determination by the combined use of x-ray and neutron-diffraction II. Composition-space refinement method. Biophys. J. 59:174–185. F(q) will therefore become poor, primarily because the simu- 5. Wilkins, M. H. F., A. E. Blaurock, and D. M. Engelman. 1971. Bilayer lated area A is less than the experimental area. For this reason structure in membranes. Nat. New Biol. 230:72–76. in part, the conclusion of Benz et al. (12) using NPT is con- 6. Kucˇerka, N., Y. Liu, N. Chu, H. I. Petrache, S. Tristram-Nagle, and J. F. siderably more critical of the CHARMM potentials than the Nagle. 2005. Structure of fully hydrated fluid phase DMPC and DLPC conclusions we draw in this article, which uses the NPAT ap- lipid bilayers using x-ray scattering from oriented multilamellar arrays and from unilameller vesicles. Biophys. J. 88:2626–2637. proach and locates the best value of the simulated surface area. 7. Liu, Y. F., and J. F. Nagle. 2004. Diffuse scattering provides material Although the potential for finite size effects is a good parameters and electron density profiles of biomembranes. Phys. Rev. reason not to impose the g ¼ 0 constraint on simulations a E. 69:040901. Biophysical Journal 90(8) 2796–2807 Interpreting X-Ray Data from Bilayers 2807 8. Wiener, M. C., R. M. Suter, and J. F. Nagle. 1989. Structure of the 24. Feller, S. E., R. M. Venable, and R. W. Pastor. 1997. Computer fully hydrated gel phase of dipalmitoylphosphatidylcholine. Biophys. J. simulation of a DPPC phospholipid bilayer: structural changes as a 55:315–325. function of molecular surface area. Langmuir. 13:6555–6561. 9. Schalke, M., P. Kruger, M. Weygand, and M. Losche. 2000. 25. Cromer, D., and J. Mann. 1968. X-ray scattering factors computed from Submolecular organization of DMPA in surface monolayers: beyond numerical Hartee-Fock wave functions. Acta Crystallogr. A. 24:321–324. the two-layer model. Biochim. Biophys. Acta. 1464:113–126. 26. Petrache, H. I., S. E. Feller, and J. F. Nagle. 1997. Determination of com- 10. Feller, S. E., D. X. Yin, R. W. Pastor, and A. D. MacKerell, Jr. 1997. ponent volumes of lipid bilayers from simulations. Biophys. J. 72:2237–2242. Molecular dynamics simulation of unsaturated lipid bilayers at low 27. Armen, R. S., O. D. Uitto, and S. E. Feller. 1998. Phospholipid com- hydration: parameterization and comparison with diffraction studies. ponent volumes: determination and application to bilayer structure calcu- Biophys. J. 73:2269–2279. lations. Biophys. J. 75:734–744. 11. Sachs, J. N., H. I. Petrache, and T. B. Woolf. 2003. Interpretation of 28. Nagle, J. F., and M. C. Wiener. 1988. Structure of fully hydrated small angle x-ray measurements guided by molecular dynamics simu- bilayer dispersions. Biochim. Biophys. Acta. 942:1–10. lations of lipid bilayers. Chem. Phys. Lipids. 126:211–223. 29. Tristram-Nagle, S., Y. F. Liu, J. Legleiter, and J. F. Nagle. 2002. 12. Benz, R. W., F. Castro-Roman, D. J. Tobias, and S. H. White. 2005. Structure of gel phase DMPC determined by x-ray diffraction. Biophys. Experimental validation of molecular dynamics simulations of lipid J. 83:3324–3335. bilayers: a new approach. Biophys. J. 88:805–817. 30. McIntosh, T. J., A. D. Magid, and S. A. Simon. 1987. Steric repulsion 13. Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. between phosphatidylcholine bilayers. Biochemistry. 26:7325–7332. Swaminathan, and M. Karplus. 1983. CHARMM: a program for macro- molecular energy, minimization, and dynamics calculations. J. Comput. 31. Jorgensen, W. L., and C. Jenson. 1998. Temperature dependence of TIP3P, Chem. 4:187–217. SPC, and TIP4P water from NPT Monte Carlo simulations: seeking tem- peratures of maximum density. J. Comput. Chem. 19:1179–1186. 14. Klauda, J. B., B. R. Brooks, A. D. MacKerell Jr., R. M. Venable, and R. W. Pastor. 2005. An ab initio study on the torsional surface of 32. Feller, S. E., Y. H. Zhang, and R. W. Pastor. 1995. Computer-simu- alkanes and its effect on molecular simulations of alkanes and a DPPC lation of liquid/liquid interfaces. II. surface- tension area dependence of bilayer. J. Phys. Chem. B. 109:5300–5311. a bilayer and monolayer. J. Chem. Phys. 103:10267–10276. 15. Durell, S. R., B. R. Brooks, and A. Bennaim. 1994. Solvent-induced 33. Skibinsky, A., R. M. Venable, and R. W. Pastor. 2005. A molecular forces between two hydrophilic groups. J. Phys. Chem. 98:2198–2202. dynamics study of the response of lipid bilayers and monolayers to trehalose. Biophys. J. 89:4111–4121. 16. Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. 1983. Comparison of simple potential functions for sim- 34. Zhang, Y. H., S. E. Feller, B. R. Brooks, and R. W. Pastor. 1995. ulating liquid water. J. Chem. Phys. 79:926–935. Computer-simulation of liquid/liquid interfaces. I. Theory and appli- cation to octane/water. J. Chem. Phys. 103:10252–10266. 17. Lagu¨e, P., R. W. Pastor, and B. R. Brooks. 2004. Pressure-based long- range correction for Lennard-Jones interactions in molecular dynamics 35. Ja¨hnig, F. 1996. What is the surface pension of a lipid bilayer mem- simulations: application to alkanes and interfaces. J. Phys. Chem. B. brane? Biophys. J. 71:1348–1349. 108:363–368. 36. Feller, S. E., and R. W. Pastor. 1996. On simulating lipid bilayers with 18. Darden, T., D. York, and L. Pedersen. 1993. Particle mesh Ewald: an an applied surface tension: periodic boundary conditions and undula- NLog(N) method for Ewald sums in large systems. J. Chem. Phys. 98: tions. Biophys. J. 71:1350–1355. 10089–10092. 37. Marsh, D. 1997. Renormalization of the tension and area expansion 19. Ryckaert, J. P., G. Ciccotti, and H. J. C. Berendsen. 1977. Numerical modulus in fluid membranes. Biophys. J. 73:865–869. integration of the Cartesian equations of motion of a system with 38. Lindahl, E., and O. Edholm. 2000. Mesoscopic undulations and thick- constraints: molecular dynamics of n-alkanes. J. Comp. Phys. 23:327– ness fluctuations in lipid bilayers from molecular dynamics simulations. 341. Biophys. J. 79:426–433. 20. Hoover, W. G. 1985. Canonical dynamics: equilibrium phase-space 39. Marrink, S. J., and A. E. Mark. 2001. Effect of undulations on surface distributions. Phys. Rev. A. 31:1695–1697. tension in simulated bilayers. J. Phys. Chem. B. 105:6122–6127. 21. Nose, S., and M. L. Klein. 1983. A study of solid and liquid carbon 40. Feller, S. E., R. W. Pastor, A. Rojnuckarin, S. Bogusz, and B. R. tetrafluoride using the constant pressure molecular-dynamics tech- Brooks. 1996. Effect of electrostatic force truncation on interfacial and nique. J. Chem. Phys. 78:6928–6939. transport properties of water. J. Phys. Chem. 100:17011–17020. 22. Andersen, H. C. 1980. Molecular-dynamics simulations at constant 41. Lamoureux, G., A. D. MacKerell, and B. Roux. 2003. A simple pressure and/or temperature. J. Chem. Phys. 72:2384–2393. polarizable model of water based on classical Drude oscillators. 23. de Vries, A. H., I. Chandrasekhar, W. F. van Gunsteren, and P. H. J. Chem. Phys. 119:5185–5197. Hunenberger. 2005. Molecular dynamics simulations of phospholipid 42. Ane´zo, C., A. H. de Vries, H. D. Ho¨ltje, D. P. Tieleman, and S. J. bilayers: influence of artificial periodicity, system size, and simulation Marrink. 2003. Methodological issues in lipid bilayer simulations. time. J. Phys. Chem. B. 109:11643–11652. J. Phys. Chem. B. 107:9424–9433. Biophysical Journal 90(8) 2796–2807