help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Originally published as Biophys J. BioFAST on March 30, 2007.
doi:10.1529/biophysj.106.087130
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
biophysj.106.087130v1
92/12/4157    most recent
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Sonne, J.
Right arrow Articles by Peters, G. H.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Sonne, J.
Right arrow Articles by Peters, G. H.
Biophysical Journal 92:4157-4167 (2007)
© 2007 The Biophysical Society

Reparameterization of All-Atom Dipalmitoylphosphatidylcholine Lipid Parameters Enables Simulation of Fluid Bilayers at Zero Tension

Jacob Sonne * {dagger}, Morten Ø. Jensen {ddagger}, Flemming Y. Hansen *, Lars Hemmingsen § and Günther H. Peters * {dagger}

* Department of Chemistry, Technical University of Denmark, Lyngby, Denmark; {dagger} MEMPHYS-Center for Biomembrane Physics, University of Southern Denmark, Odense, Denmark; {ddagger} Department of Life Sciences and Chemistry, Roskilde University, Roskilde, Denmark; and § Department of Natural Sciences, The Faculty of Life Sciences, University of Copenhagen, Copenhagen, Denmark

Correspondence: Address reprint requests to Günther H. Peters, E-mail: ghp{at}kemi.dtu.dk.


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
Molecular dynamics simulations of dipalmitoylphosphatidylcholine (DPPC) lipid bilayers using the CHARMM27 force field in the tensionless isothermal-isobaric (NPT) ensemble give highly ordered, gel-like bilayers with an area per lipid of ~48 Å2. To obtain fluid (L{alpha}) phase properties of DPPC bilayers represented by the CHARMM energy function in this ensemble, we reparameterized the atomic partial charges in the lipid headgroup and upper parts of the acyl chains. The new charges were determined from the electron structure using both the Mulliken method and the restricted electrostatic potential fitting method. We tested the derived charges in molecular dynamics simulations of a fully hydrated DPPC bilayer. Only the simulation with the new restricted electrostatic potential charges shows significant improvements compared with simulations using the original CHARMM27 force field resulting in an area per lipid of 60.4 ± 0.1 Å2. Compared to the 48 Å2, the new value of 60.4 Å2 is in fair agreement with the experimental value of 64 Å2. In addition, the simulated order parameter profile and electron density profile are in satisfactory agreement with experimental data. Thus, the biologically more interesting fluid phase of DPPC bilayers can now be simulated in all-atom simulations in the NPT ensemble by employing our modified CHARMM27 force field.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
METHODS
 RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
In molecular dynamics (MD) simulations of biological membranes, lipid bilayers are the most common model system, and with the increase in computer power, simulations of larger and more complex bilayer systems are becoming accessible (1Go–8Go). To represent bilayers accurately in MD simulations, it is important to have a well-parameterized force field in combination with appropriate macroscopic boundary conditions (ensembles). Although these topics may seem somewhat technical in a biological context, the functional importance of cell membranes should not be overlooked. For example, it has been shown that phospholipase activity changes the mechanical properties of bilayers (9Go), which could have therapeutic application in delivery of cytotoxic anti-cancer drugs (10Go). Further, bilayer properties affect mechanosensitive gating in McsL (11Go–13Go) and phospholipase C activity (14Go), and have been suggested to be related to anesthesia (15Go,16Go). See Jensen and Mouritsen (17Go) for a recent review of bilayer influence on protein function.

Commonly used united-atom force fields produce fluid (L{alpha}) phase bilayers in the isothermal-isobaric NPT ensemble, provided that the correct electrostatic cutoff strategy is used (18Go–21Go). In this study, the NPT-notation indicates that the simulation box length in the direction normal to the bilayer z is coupled to the barostat independently from the lateral directions x and y and the pressure tensor holds, Px = Py = Pz. Therefore the surface tension {gamma} is zero (1Go). As explained by Berger et al. (22Go), flaccid bilayers in experiments are able to adjust their area per lipid A. Minimizing the contact between the acyl chains and water and maximizing the acyl-chain entropy creates a free energy minimum, which together dictate the equilibrium area per lipid Aeq, i.e., Formula Therefore, it has been argued that the (tensionless) NPT ensemble is appropriate for simulating lipid bilayers (22Go–24Go).

In contrast, bilayers simulated using the all-atom CHARMM22 force field (25Go,26Go) in the NPT ensemble show a dramatic lateral contraction and overly ordered lipid acyl chains (1Go). Despite optimizations in both the headgroup and the acyl chains in the subsequent CHARMM27 parameter set (27Go,28Go), the area per DPPC lipid molecule is underestimated by at least 15 Å2 (2Go), compared to experimental data (29Go). In fact, the area per DPPC lipid in NPT simulations (2Go) is close to that of the gel phase (4Go). Recently, the CHARMM27 force field was further optimized (30Go), but the revised parameters for the acyl chains do not produce fluidlike bilayers in the NPT ensemble (31Go). Thus, to mimic the biologically relevant fluid phase, some bilayer simulations using CHARMM parameters apply a positive surface tension (1Go,2Go,5Go,13Go,28Go,30Go,32Go). It has been suggested that the gel-like properties in NPT simulations can be attributed to finite size effects, and therefore, applying a positive surface tension, which stretches the bilayer to the experimentally determined area per lipid, is appropriate (1Go,28Go,33Go). However, finite size effects seemingly account for an area contraction of less than 1 Å2/lipid (19Go), yet, the area per lipid is underestimated by at least 15 Å2 compared to experimental data for a DPPC bilayer (72 lipids) simulated in the NPT ensemble with CHARMM27 parameters (2Go). Thus, even though MD simulations of bilayers are subject to finite size effects, these can by no means account for the gel-like bilayer properties in NPT simulations using CHARMM27 parameters. It appears that the gel-like properties are mostly a result of the force field not being optimized for lipid bilayer simulations in the NPT ensemble, but for simulations with an applied surface tension.

Although there is no agreement on whether application of a surface tension is appropriate when simulating lipid bilayers, tensionless NPT simulations are appealing from a practical point of view since the area per lipid need not be known before simulation, if the parameters are optimized for this ensemble. The advantage of having parameters that are optimized for NPT simulations becomes more prominent when, for example, studying lipid mixtures and bilayers with embedded proteins where experimental data are scarce. Therefore, realistic NPT bilayer simulations have been one ultimate goal for MD force fields (30Go,34Go).

Along these lines, the purpose of this study is a reparameterization of the CHARMM phosphatidylcholine (PC) lipid parameters, which permits simulations of fluid phase DPPC lipid bilayers in the NPT ensemble. The lipid headgroup charge is the target for our reparameterization. Our article is structured as follows: In Methods, we describe the parameterization strategy and in Simulation Details, we give the technical details of our calculations and simulations. We present and discuss our results in Results and Discussion, respectively, and summarize our findings in Summary and Conclusions. The partial charges resulting from our study are provided in an Appendix.


   METHODS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
We start this section by outlining relevant aspects of the CHARMM optimization strategy before we describe the methods used in this study.

CHARMM27 parameterization strategy
In the CHARMM force field, electrostatic interactions are accounted for by partial point charges located at the atomic centers (26Go–28Go,35Go,36Go). Since electrons are delocalized, this representation cannot be exact and there is no rigorous way of determining such atomic charges. Nevertheless, finding suitable partial charges that, with reasonable precision, account for molecular properties is one task in parameterizations of this force field. Complying with the CHARMM strategy, Foloppe and Mackerell, Jr. (27Go) optimized the CHARMM partial charges of dimethylphosphate to reproduce quantum mechanical interaction energies with a strategically placed water molecule using Mulliken population analysis to propose the initial partial charges. The optimized charges are now used for the phosphate moiety in the headgroup of phospholipids (28Go).

In the optimization of the CHARMM27 parameters, Foloppe and MacKerell, Jr. (27Go) distinguish between macromolecular target data and small molecule target data. The ability of the CHARMM27 force field to reproduce the former has the highest priority. In the case of lipid bilayers, macromolecular target data would be bilayer properties such as lipid densities, electron densities, and order parameters, and small molecule target data would be, e.g., torsional energy surfaces in the lipid headgroup and water-lipid interaction energies.

Our parameterization strategy
Pressure profile calculations for lipid bilayers suggest that the equilibrium area per lipid is determined by a delicate balance between large and opposing forces originating from bonded and nonbonded interactions (37Go,38Go). Therefore, reparameterization of any of the terms in the energy function could, in principle, affect the area per lipid. NPT and NPz AT simulations of a crystalline (all-trans) C36 alkane showed good agreement with experiments (39Go,40Go), which indicates that the alkane-alkane interactions are well parameterized. With the recent refinement of the torsional potential for CHARMM27 alkanes, the parameters for the acyl chains of phospholipids also seem highly optimized (30Go). Our attention therefore turns to the lipid headgroup region. The inter-headgroup interactions are determined by the Lennard-Jones parameters and the partial charges. We believe the latter to be a more promising optimization target, since numerous united-atom force-field bilayer simulations have shown that the bilayer properties are quite sensitive to details in the treatment of the electrostatic headgroup interactions (18Go–20Go). Moreover, pressure profiles derived from atomistic simulations indicate that electrostatic attractions significantly contribute to the positive surface tension in the headgroup region. This indicates that the electrostatic forces in this region are, on average, contractive (37Go,38Go). Therefore, the lipid headgroup charges will be the target of our reparameterization.

We determine initial partial charges of the whole lipid headgroup and upper acyl chain from ab initio data, using (I) Mulliken population analysis and (II) a restricted electrostatic potential (RESP) fitting approach (3Go), i.e., we adjust the partial atomic charges to fit the quantum mechanical electrostatic potential (see Simulation Details for further details). The RESP charges generate a realistic electrostatic potential around the molecule of interest and the method is optimal for reproducing intermolecular interactions (3Go,41Go,42Go). However, determination of the Mulliken charges does not add any significant overhead to the calculation of the quantum mechanical electrostatic potential and we therefore include both methods in this study.

In the RESP method, the derived charges are known to be strongly dependent on the conformation of the molecule. Reynolds et al. (43Go) addressed this problem by deriving RESP charges from different conformers and then estimated the final set of partial charges as the Boltzmann-weighted average of the charge sets found for the different conformers. We adopt a similar approach by extracting 69 dipalmitoylphosphatidylcholine (DPPC) configurations from a 15 ns MD simulation of a DPPC lipid bilayer, where the area was fixed at the experimental value of 62.9 Å2 per lipid (29Go). This experimental estimate was later adjusted to 64 Å2 (4Go,44Go), but for extracting lipid conformations, this adjustment should be of minor importance. These 69 DPPC molecules were capped to form dipentanoatephosphatidylcholine (DPePC), shown schematically in Fig. 2. Since numerous studies have shown that the CHARMM alkane parameters are highly optimized (28Go,30Go,39Go,40Go), we did not include the complete acyl chains in our ab initio calculations. From the 69 DPePC configurations, we determined the Mulliken and RESP charges and calculated the final charges as a simple average over the conformers, which are already Boltzmann-weighted from the MD simulation.


Figure 2
View larger version (44K):
[in this window]
[in a new window]

  		FIGURE 2  Schematic structure of dipentanoatephosphatidylcholine (DPePC) with indication of the average atomic charges determined by (A) Mulliken population analysis and (B) the RESP procedure. Atoms in the lipid backbone are black. The atoms in the three methyl groups on nitrogen (N3) have identical charges, but for clarity only one methyl group is shown explicitly. A bar indicates that the charge is negative. The atom labels also indicate the symmetry constraints that were used. For example, the two atoms labeled H5 have the same charge by definition. The CHARMM27 charge group definition is shown with shaded boxes and the total charge of the groups are given by the qXs (X is the group number I–VII). The CHARMM27 group charges are shown in parentheses. The two dashed lines indicate transfer schemes A and B (see Our Parameterization Strategies for details). Note that the precision of the charges leaves DPePC with an apparent nonzero net charge due to round-offs. This is not the case when all six decimals are included (see Table 2 and Appendix).

 
When changing the atomic partial charges, the remainder of the force-field parameters could be readjusted iteratively to fix, e.g., ab initio potential energy surfaces and vibrational data, i.e., small molecule target data (27Go). However, our main motivation for assigning new partial charges to DPPC lipids is to obtain fluidlike properties of DPPC bilayers, i.e., to optimize the macromolecular (bilayer) properties only. Therefore, testing whether our new charges reproduce small molecule target data is beyond the scope of this work. Our new DPePC partial charges were therefore transferred to DPPC without modifying the force field further. We used two schemes for this transfer:
Scheme A. All charges from DPePC were used in DPPC except for the terminal C26 methyl groups (see Fig. 2).
Scheme B. All charges from DPePC were used in DPPC except for the terminal C24-C26 ethyl groups.

Excluding the methyl and ethyl groups in Scheme A and B, respectively, leaves DPePC, and therefore also DPPC, with a nonzero net charge. To obtain a neutral DPPC molecule we compensated this net charge with small and equal counter charges on each nonacyl atom in DPPC. We choose not to cluster our charges in charge groups, with integral values for the group charges, since most membrane simulations using the CHARMM force field evaluate electrostatic interactions by PME summation, rendering charge clustering superfluous. The lack of charge groups abrogates transferability of our charges to other phospholipids such as phosphatidylethanolamines.

Testing the new parameters
Combining the Mulliken population analysis (I) and RESP (II) with transfer Scheme A and B yields four sets of partial charges, which we subsequently tested in four MD simulations of a DPPC lipid bilayer. The simulation results, referred to as I.A, I.B, II.A, and II.B, were compared with experimental data and with NPT and NPz{gamma}T simulations using the CHARMM27 parameters.

As benchmark properties for the comparison, we resorted to the volume and area per lipid, the order parameter profile, and the electron density profile. The volume per lipid is calculated by subtracting the water volume nwVw from the box volume (45Go), where nw is the number of water molecules and Vw is the average volume of one bulk water molecule in the simulation. The area per lipid is obtained as the area of the simulation box xy-plane divided by the number of lipids in one leaflet. The order parameter for the ith acyl methyl(ene) group |SCD,i| is calculated as |SCD,i| = |<3/2 cos2 {theta}i – 1/2>|, where {theta}i is the angle between a C-H bond vector in the ith methyl(ene) group and the bilayer normal, i.e., the z-axis. The brackets denote averaging over the C-H bonds in the ith methyl(ene) group, lipids, and time. The electron density profile was calculated by binning the difference Z – q for all atoms along the z-axis, where Z is the atom number and q is the atomic partial charge. Thus, we neglect bilayer undulations and assume that the electrons are located at the atomic centers, thereby ignoring variations in the atomic form factors. Both assumptions are reasonable when simulating a small fluid bilayer patch (32Go). To elucidate the changes in the lipid structure and in the hydration caused by the reparameterization, we have calculated pair distribution functions for selected atoms pairs.

Atomic partial charges
Without geometry optimization, the quantum mechanical electrostatic potential around ~69 DPePC conformers was evaluated at the RHF/6-31G(d) level in ~100,000 grid points using Gaussian98 (46Go). In addition to the electrostatic potential, we also calculated the atomic charges from a Mulliken population analysis. From the quantum mechanical electrostatic potential, atomic partial charges were determined by two successive restricted electrostatic potential (RESP) fits. In the first fit, we used the CHARMM27 charges as the initial guess and used no symmetry constraints. The charges were restrained by a hyperbolic penalty function with weight 0.0005 AU to avoid large charge separation. In the successive fit, the output charges from the first fit were used as input, and the hyperbolic penalty weight was increased to 0.001 AU. In the second fit we enforced symmetry constraints for equivalent atoms (see Fig. 2). The same definition of equivalent atoms was used to symmetrize the charges from the Mulliken population analyses. Lastly, we averaged partial charges over the 69 DPePC configurations to give the average partial charge distribution in the DPePC molecule from the Mulliken and RESP procedures, respectively. The RESP fitting was performed using AMBER 4.1 with the Restrained ESP Fit package 2.3 (47Go).

MD simulations
The DPPC lipid bilayer, from which we extracted the lipid conformations for the atomic charge determination, consists of 72 lipids solvated with ~2000 water molecules resulting in a total of ~16,000 atoms. This corresponds to ~29 water molecules per lipid (29Go). The water molecules were placed around the bilayer using Solvate (48Go) and subsequently the water layer surrounding the bilayer was cropped to a rectangular, periodic simulation box. The water molecules were represented by the TIP3 water model (49Go). The bilayer system was equilibrated in the CHARMM27 force field for 15 ns with the area fixed at the experimental value, 62.9 Å2 (29Go). After this equilibration, we changed the partial charges to sets I.A, I.B, II.A, and II.B in four simulations, which were then continued. We also carried out two reference simulations, referred to as III and IV, where the original CHARMM charges were used. The simulation times and pressure coupling schemes are summarized in Table 1. In all simulations, we used a time step of 1.0 fs and the target temperature of the Langevin thermostat was 325 K with a damping coefficient of 5 ps–1. The pressure was controlled by the Nosé-Hoover Langevin barostat (50Go) with a piston oscillation time of 100 fs and a damping time of 50 fs. Electrostatic interactions were evaluated using the PME method (51Go,52Go) with a grid spacing below 1 Å. MD simulations were carried out using NAMD (53Go).


View this table:
[in this window]
[in a new window]

  		TABLE 1  Overview of MD simulations

 

   RESULTS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
Atomic partial charges
In Fig. 1, we present average charges and standard deviations for Mulliken and RESP procedures. The charges are also compiled in Table 2. To get a better overview of the charge distribution, the obtained charges are mapped onto DPePC in Fig. 2 with only two decimals, which leaves DPePC with a nonzero net charge due to round-off errors. With six decimals the new charges fulfill electroneutrality (see Table 2). Qualitatively, there is a good agreement between the CHARMM27 atomic charges and the charges obtained from the Mulliken method. This agreement is expected since the initial CHARMM partial charges were taken from a Mulliken population analysis (25Go). In general, the same is true for the obtained RESP charges; however, the RESP charges for N3, C6, and C24 have the opposite sign compared to the corresponding CHARMM27 charges. A similar sign inversion relative to CHARMM27 is not observed for our Mulliken charges and therefore this prediction appears to be specific for the RESP method. Comparing the error bars on the charges from the two methods, we see that the RESP charges exhibit a larger configurational dependence than the Mulliken charges. This is a commonly known difference between the two methods (3Go) and was part of our rationale for using several DPePC configurations.


Figure 1
View larger version (20K):
[in this window]
[in a new window]

  		FIGURE 1  Overview of the obtained partial charges for dipentanoatephosphatidylcholine (DPePC). Average values and standard deviations for the new RESP charges and new Mulliken charges are determined from 69 Boltzmann-weighted DPePC conformers. For comparison, we have included the charges from the CHARMM27 force field. The atom IDs used on the abscissa axis are defined in Fig. 2.

 

View this table:
[in this window]
[in a new window]

  		TABLE 2  Partial charges for DPePC and DPPC in different parameter sets

 
Even though we do not make use of charge groups in our simulations, we calculated the net charge of the CHARMM27 charge groups for the new charge sets. These are shown in Fig. 2 by the shaded boxes and the corresponding group charges are denoted qX (X is the group number I–VII). The group charges for CHARMM27 are shown in parentheses. Even though the RESP method predicts that N3 has a positive charge, the direction and magnitude of the NP-dipole moment is similar to that predicted by the Mulliken method. Comparison of qI and qII indicates that the new charges, both Mulliken and RESP, give a smaller NP-dipole moment compared to CHARMM27 charges. Interestingly, the methylene group charges (qV and qVI) are no longer neutral and therefore do not have alkane properties as in the CHARMM27 parameter set.

Since there is no rigorous definition of atomic partial charges, any of the three charge distributions in Fig. 1 can be valid and we need to test their quality in lipid bilayer simulations. Table 2 and Fig. 7 provide an overview of our charges after transfer to DPPC.


Figure 7
View larger version (25K):
[in this window]
[in a new window]

  		FIGURE 7  Overview of the new partial charges for the DPPC lipids as obtained from the RESP method. (A) Transfer scheme A where the atomic charges of C1 to C4 in the acyl chain were modified compared to CHARMM27. (B) Transfer scheme B where the atomic charges of C1 to C3 in the acyl chains were modified compared to CHARMM27. In both panels the atoms in the lipid backbone are shaded. The atoms in the three methyl groups on nitrogen have identical charges and for clarity, only one methyl group is shown explicitly. The solid zigzag line in the sn2 chain indicates that these atoms have been assigned new charges, which are equal to the corresponding charges in the sn1 chain. For clarity, the atoms are only shown explicitly in the sn1 chain. The dashed lines in the sn1 and sn2 chains symbolize the rest of the acyl chain where the CHARMM27 charges were used. A bar indicates that the charge is negative. Note that rounding of the charges in the figure leave DPPC with a nonzero net charge, which is not the case when all six decimals are included (see Table 2).

 
Testing new partial charges
We tested parameter sets I.A, I.B, II.A, and II.B in MD simulations of a DPPC lipid bilayer. For reference, we have also included results from two simulations using the CHARMM27 force field, one with no applied surface tension and one with an applied surface tension of 61 mN/m (2Go). These reference simulations are referred to as III and IV, respectively. All six simulations are summarized in Table 1. In the following we compare the area per lipid, the volume per lipid, the order parameter profile, and the electron density profile for these six simulations mutually and with experiments.

Area and volume per lipid
From Fig. 3 we see that the areas per lipid in simulations I.A, I.B, and II.B monotonically decrease to ~55 Å2 within the first 2 ns of the simulations, resembling the behavior of the first 2 ns of simulation III. In simulation III, the area per lipid reaches the experimental value for the gel phase (4Go) after 15 ns and is still decreasing moderately. Based on the immediate decrease in the area per lipid observed in simulations I.A, I.B, and II.B we decided to end these three simulations and to discard them from further analysis. In contrast, the area per lipid found in simulation II.A is stable and the average value of 60.4 ± 0.1 Å2 compares favorably to the commonly accepted experimental value of 64 Å2 (4Go). As previously reported, it is necessary to apply a surface tension of 61 mN/m to obtain an area per lipid of 64.5 ± 0.3 Å2 using the CHARMM27 force field (simulation IV) (2Go). The error estimates for the area per lipid are obtained as the standard error of the mean area calculated in 250 ps data blocks.


Figure 3
View larger version (13K):
[in this window]
[in a new window]

  		FIGURE 3  Area per lipid as function of time for six simulations using different parameter sets and ensembles (see Table 1). The curves are 200 point-running averages of the area per lipid calculated every 0.5 ps. The numbers on the right ordinate axis are experimental values for the area per lipid (4Go).

 
Simulations II.A, III, and IV give lipid volumes of (12.0 ± 0.3) x 102, (11.7 ± 0.6) x 102, and (12 ± 1) x 102 Å3/lipid. The experimental estimate is ~12.30 x 102 Å3/lipid and varies a few Å3/lipid depending on the technique used (4Go). Thus, all the simulation results agree with the experimental data within the statistical uncertainty in the simulations. In the determination of the lipid volume, the water volume was calculated from the water electron density in Fig. 5 A. Simulation III is not completely equilibrated (see Fig. 3), and therefore we used only the last 3 ns to estimate the lipid volume in this simulation. In the following sections we also use only the last 3 ns of this simulation for the data analysis.


Figure 5
View larger version (16K):
[in this window]
[in a new window]

  		FIGURE 5  Electron density data. (A) Electron density profiles from simulations II.A, III, and IV. (B) Absolute form factors obtained by Fourier transform of panel a. We have included estimates of the absolute form factors obtained from experiments (29Go). These were scaled such that Fexp(Formula) = F sim(Formula) (open circles). In panel a, the standard error in the mean, estimated from the variations in the average electron density profiles calculated in 250 ps time blocks, is below 0.01 e/Å3 in all three simulations.

 
Order parameters
In Fig. 4 we show the order parameter profile for the sn2-chain from simulations II.A, III, and IV as well as experimentally determined deuterium order parameters at 41°C and 50°C (54Go,55Go). The experimental profile measured at 50°C is close to the simulation conditions. The profile obtained at 41°C is the upper limit for |SCD| in a fluid phase DPPC bilayer, since the main phase transition temperature is 41°C (55Go). In the region from C7 to C16, the profile from simulation II.A resembles most closely the experimental profile for the fluid phase found at 41°C. In the region nearer to the glycerol backbone (C3-C6), the order parameters in simulation II.A are lower than the experimental values at 50°C. Such deviations are not found in simulation IV (CHARMM27 parameters with {gamma} = 61 mN/m), which resembles the experimental profiles at 50°C quite closely. The order parameter profile calculated from simulation III (CHARMM27, {gamma} = 0 mN/m) clearly indicates highly ordered acyl chains characteristic for the gel phase. Overall, the changes in the sn1 chain (not shown) are similar to the changes presented for the sn2 chain. The order parameter for the C2 carbon in the sn1 chain is 0.167 ± 0.002 in simulation II.A and 0.138 ± 0.003 in simulation IV.


Figure 4
View larger version (11K):
[in this window]
[in a new window]

  		FIGURE 4  The order parameter profiles from experiments (Exp.) and three simulations (II.A, III, and IV) for C2 to C16. Experimental data from Douliez et al. (54Go). For C2 the data refer to the pro-R position. For simulation III, the order parameters are calculated from the last 3 ns. The standard error in the mean, estimated from the variations in the average order parameter profiles calculated in 250 ps time blocks, is <0.005 for all methylene groups in all three simulations.

 
Electron density profile
Fig. 5 A shows the average electron density profile from simulation II.A, III, and IV. The overall shapes of the three profiles are similar and the characteristic methyl troughs in the bilayer center and the headgroup peaks are apparent in all three simulations. The profile from simulation III has more sharp features than the two other profiles, but the differences between simulations II.A and IV are more subtle. To obtain the bilayer form factors Fsim(q), we Fourier-transform the real-space profiles in Fig. 5 A, which yields the Fsim(q) curves in Fig. 5 B (32Go). We have also included estimates of the absolute bilayer form factors from x-ray scattering experiments (29Go) and fitted these to our simulated F(q). For the fitting, we scaled the experimental first order (h = 1) form factors Fexp(Formula) to the simulated form factors Fsim(Formula) and scaled the higher order form factors of the given experimental sample by the same factor. For h ≥ 2, the experimental F(q) is larger than observed in both simulations. The root mean-square deviations of the fitted experimental data points from the corresponding simulated F(q) are 0.4, 1.1, and 0.2 e/Å2 for simulations II.A, III, and IV, respectively. Judged by the form factors, simulation IV therefore exhibits the best agreement with the experimental data.

Pair distribution functions
In general, the pair distribution function g(r) from the fourth acyl carbon (C24) to nonhydrogen headgroup atoms within the same lipid exhibits no qualitative differences when comparing simulations II.A and IV (data not shown). The only exception is g(r) for the atom pair C24sn1–C20sn2 shown in Fig. 6 A. In simulation II.A, this g(r) has a bimodal structure with two peaks at 5.5 Å and 7.5 Å, respectively. In simulation IV there is only one peak in at 5.4 Å, i.e., close to the position of one of the peaks in the II.A simulation. We have also calculated the function gxy(r), which is a pair-distribution function based on the xy-projection of the distance vector between the atoms. In Fig. 6 A, gxy(r) is shown for simulations II.A and IV. We find essentially no change in gz(r), a pair distribution function based on the z-projection of the distance vector between the atoms, when comparing simulation II.A and IV (data not shown).


Figure 6
View larger version (10K):
[in this window]
[in a new window]

  		FIGURE 6  Pair distribution functions g(r). (A) Between the fourth acyl carbon (C24sn1) and the second acyl carbon (C20sn2) in the same lipid. The value g(r) is the ordinary pair distribution function with r equal to the distance between the atoms and gxy(r) is the pair distribution function where r is the projection of the distance between atoms on the xy-plane. (B) Between C24 (fourth acyl carbon) and oxygen in water. The value g1(r) is the pair distribution function averaged over the two chains. The value g2(r) is the pair distribution function between the center of mass of the two C24s in one lipid and oxygen in water. All g(r)s are averaged over lipids and time.

 
The pair distribution function g1(r) from the fourth acyl carbon (C24 in Fig. 2 and C4 in Fig. 7) to the water oxygen is shown in Fig. 6 B. The value g1(r) for simulations II.A and IV indicate that the hydration of the glycerol backbone region is higher in simulation II.A. Fig. 6 B also shows the pair distribution function g2(r) between the center of mass of the two C24s (one in each acyl chain) and water. Consistent with g1(r), g2(r) is higher in simulation II.A compared to simulation IV. The standard error for g2(r), for distances <2 Å, is of the same magnitude as g2(r) in simulation IV.

We find only minor differences in the orientation of water molecules within 6 Å of C24 when comparing simulations II.A and IV (data not shown).


   DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
MD simulations of DPPC lipid bilayers using the all-atom CHARMM27 parameter set yield highly ordered bilayers with gel-like properties in the NPT ensemble (1Go,2Go). Since finite size effects only account for a part of the bilayer ordering, it appears that the CHARMM27 lipid parameters need optimization provided that one resorts to NPT simulations.

In this study we determined four new sets of partial charges for the DPPC lipid headgroup from Mulliken population analysis and from the RESP fitting procedure. These procedures do not comply with the traditional strategy used for optimizing partial charges in the CHARMM force fields (25Go–28Go,36Go). To this end we note that Foloppe and Mackerell, Jr. (27Go) explain that the CHARMM27 parameters are "primarily optimized to reproduce macromolecular target data while maximizing agreement with small molecule target data". In our opinion, the macroscopic target data are the key properties of the bilayer and consequently a correct representation of these must be the ultimate goal of any simulation deploying CHARMM lipid parameters. Thus, even though the methods we used to obtain the new atomic partial charges deviate from the traditional CHARMM parameterization strategy, this should be of minor importance if the new parameters represent DPPC bilayers more accurately than the conventional CHARMM27 parameters.

We tested the four new sets of partial charges in MD simulations using the charges presented in Table 2, i.e., with six decimals. Using atomic partial charges with six decimals in the MD simulations should not be crucial for the properties of the bilayer. However, keeping this precision was convenient since the round-off errors introduced to fulfill electroneutrality were found smaller than the configurational dependence of the atomic charges.

Which parameter set?
Fig. 3 shows that within 15 ns, the area per lipid in simulation III decreases to ~48 Å2, which is close to the experimental value for the DPPC gel phase (4Go). Based on the similarity of the first few nanoseconds, the same behavior is expected for simulations I.A, I.B, and II.B. Thus, further equilibration of these simulations would most likely give bilayers with gel-like properties. Since the goal of this study is to develop a new parameter set that gives more fluidlike bilayers in the NPT ensemble, we stopped simulations I.A, I.B, and II.B. The area per lipid in simulation II.A oscillates at an approximate average value of 60.4 Å2, and considering the variations in the area per lipid as determined from different experiments (4Go), this value is satisfactory. The area per lipid from simulation II.A is also in excellent agreement with the results of MD simulations using other force fields (18Go,19Go,22Go). To obtain an area per lipid that is close to the experimental value of 64 Å2 using the CHARMM27 parameters, it is necessary to apply a surface tension of 61 mN/m (simulation IV) (2Go). The volume per lipid is comparable to the experimental value of 1230 Å3 (4Go) for all three simulations.

The acyl chains in simulation III were found to be highly ordered when compared to the experimental data for a fluid (L{alpha}) phase bilayer, which is consistent with the underestimation of the area per lipid in this simulation. The chain order data are similar to results from bilayer simulations using the CHARMM22 parameter set with zero applied surface tension (1Go). The order parameter profiles for simulations II.A and IV indicate that these simulations are in the fluid phase. In simulation II.A, the lower part of the acyl chains is slightly more ordered than expected from experiments, but is still fluidlike when compared to simulation III. In the region C3-C5 the chain order is somewhat underestimated. In the region C4 to C6, where experiments cannot resolve the order parameter profile (54Go), the shape of the II.A profile deviates from that resulting from simulation IV and deviates from what seems to be the generic shape of simulated order parameter profiles (2Go,20Go,22Go,28Go). The profile computed from simulation IV resembles the experimental profiles at 50°C and reproduces the plateau from C4 to C6, which indicates that the atypical shape of the II.A order parameter profile from C4 to C6 is an effect of the II.A parameters rather than a consequence of melting of the acyl chains.

Compared to simulations II.A and IV, the electron density profile from simulation III shows sharp headgroup peaks, a deep narrow methyl trough in the bilayer center and plateau regions in between. These well-defined features indicate that the positions of the headgroup phosphate moiety and the terminal methyl groups are relatively well defined, which is consistent with the overall picture of a very ordered bilayer structure in this simulation. The sharp features in the electron density profile of simulation III are also reflected in the form factor FIII(q), which is nonzero for q > 1 Å–1 and resembles the experimental gel phase F(q) (4Go). However, the bilayer in simulation III is not in a fully developed gel (Lß') phase since the acyl chains are not tilted relative to the bilayer normal. The difference between the electron density profiles from simulations II.A and IV are subtle and therefore the electron density profiles are compared with experimental data in Fourier space. Fitting the experimental estimate for the absolute form factors to the simulated factors, we find that the experimental data points consistently lie slightly higher than predicted by both simulations II.A and IV which, except for the third-order data, is consistent with the findings of Sachs et al. (56Go). Overall, the comparison of absolute form factors indicates that simulation IV conforms more closely to x-ray scattering experiments (29Go) than the results from simulation II.A. For this comparison, we used the currently available experimental data. However, new and improved x-ray results for DPPC show much better agreement with both simulations II.A and IV than the data used for comparison here (44Go). For example, it appears that the experimental form factors we have used for the fitting in Fig. 5 B are overestimated by ~0.3 e/Å2 and ~0.5 e/Å2 for h = 3 and h = 4, respectively. Until the new DPPC x-ray data is released for general use, we cannot carry out a more thorough comparison to determine which one of simulations II.A or IV reproduces the new experimental data the best.

The II.A parameters were previously used in a pressure profile study (38Go). The pressure profiles reported in that study qualitatively resemble pressure profiles calculated from other atomistic force fields (13Go,37Go). Since pressure profiles cannot be measured experimentally, we will settle with this qualitative agreement.

The origin of the fluid phase
For future optimizations studies of the CHARMM lipid parameters, it is useful to pin down why the II.A parameters improve the bilayer properties. It is likely that simulations I.A, I.B, and II.B would eventually attain gel-like properties as does simulation III, which indicates that area per lipid is relatively insensitive to most of the changes that we have made in the headgroup, such as the inversion of the sign on the N3 and C6 charges in parameter set II.B. Apparently, the reduction of the group charge of the choline and phosphate moieties (i.e., reduction of the NP-dipole moment) in simulations I.A, I.B, and II.B does not affect the area per lipid either. As seen in Fig. 7, the only difference between simulations II.A and II.B is essentially the partial charge of the C24 methylene group (see Fig. 2). Since only simulation II.A is able to maintain fluid phase properties of the bilayer, the fourth methylene group seems to be responsible for the fluid phase properties of the DPPC bilayer in that simulation.

Changes in the lipid configurations were analyzed through pair distribution functions g(r) calculated between the seemingly important fourth acyl carbon (C24) and other nonhydrogen headgroup and upper acyl-chain atoms in the same lipid. We only compare simulations II.A and IV since the lateral density (area per lipid) is approximately the same and therefore any differences in g(r) can be attributed to the new partial charges rather than differences in g(r) between fluid and gel-like bilayers, which are less important here. The pair distribution functions in general reveal no qualitative structural changes in the lipid configurations when comparing simulations II.A and IV. The only exception is the C24sn1–C20sn2 g(r), where it appears that one effect of introducing the II.A parameters is to stabilize a new lipid configuration at r ~ 7.5 Å in addition to the state that dominates simulation IV at r ~ 5.5 Å. The lateral pair distribution functions indicate that in simulation II.A, the lateral distance between the upper parts of the sn1 and sn2 chains has increased, which supposedly contributes to the increased fluidity of the resulting state of the bilayer. It is also interesting to note that the hydration of the glycerol backbone is higher in simulation II.A compared to simulation IV even though the area per lipid is slightly higher in the latter. Further, since the center of mass of the two C24s in one lipid should, on average, lie between the two acyl chains, the nonzero value of g2(r) for distances below 3 Å in simulation II.A indicates that water is occasionally found between the two acyl chains. This is essentially never the case in simulation IV. This result conforms well to the increased separation between the upper acyl chains as indicated by the two intralipid g(r)s. The increased hydration of the glycerol backbone region in simulation IIA may also be important for the fluidity of the bilayer.

Thus, we have identified several structural differences between simulations II.A and IV relating both to the lipid structure and the level of hydration. However, the analysis does not allow us to decide whether these structural changes alone are responsible for the enhanced fluidlike properties of the bilayer in simulation II.A.

Future optimization
In summary, our investigations of the area and volume per lipid, the order parameter profile and the electron density profile for DPPC clearly indicate that, in the NPT ensemble, the II.A parameters reproduce fluid phase bilayer properties better than the CHARMM27 parameters. The volume per lipid compares favorably with experiments for simulations II.A, III, and IV. Further, we find that for a pure DPPC lipid bilayer simulated with CHARMM27 charges, fluidlike properties can also be obtained by applying an appropriate positive surface tension which, however, requires that the area per lipid has been predetermined experimentally.

Still, in simulation II.A the area per lipid is slightly underestimated and the order parameter is slightly overestimated in the lower acyl chain region compared to experimental data. This could be due to finite size effects in the simulation (19Go), but it could also indicate that further fine-tuning of the lipid parameters in the glycerol backbone region is needed to approach the experimental data even further. Indeed, the extreme and unexpected sensitivity to the charges in the fourth methylene group stresses the importance of obtaining good lipid parameters in the glycerol backbone region. In, e.g., the GROMOS96 45A3 force field, the area per lipid is also sensitive to the charges in this region (57Go). Since the glycerol backbone is a common motif in all glycerophospholipids, fine-tuning the CHARMM parameters in this region could be useful also for other lipid systems. Adjusting, for example, the acyl-chain Lennard-Jones parameters could also prove fruitful (22Go), but then the acyl-chain parameters would depart from the highly optimized CHARMM27 alkane parameters.


   SUMMARY AND CONCLUSION
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
As previously described in the literature, we confirm that dipalmitoylphosphatidylcholine (DPPC) lipid bilayer simulated in the NPT ensemble using the CHARMM27 force field have gel-like properties while fluid (L{alpha}) phase properties can be obtained by applying an appropriate positive surface tension (1Go,2Go,28Go,30Go). Since the applied surface tension has to be adjusted based on experimental data (1Go,2Go), we believe that this approach makes it difficult to take full advantage of the predictive power of MD simulations. In an attempt to obtain fluid phase properties of DPPC bilayers simulated in the NPT ensemble using the CHARMM energy function, we assigned new partial charges to the headgroup and upper acyl chains of DPPC using Mulliken population analysis and a RESP fitting procedure. The new Mulliken partial charges do not have an immediate effect on the bilayer properties, but the RESP charges do. Using the new RESP charges, we find a dramatic improvement of the bilayer properties compared to simulations using the CHARMM27 charges with zero applied surface tension. Thus, the new RESP partial charges presented in this study allow for simulating a DPPC lipid bilayer in the fluid phase at constant pressure and zero applied surface tension using the CHARMM energy function.


   APPENDIX: NEW PARTIAL CHARGES
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
ACKNOWLEDGEMENTS
 REFERENCES
 
Table 2 shows the average Mulliken and RESP charges obtained for dipentanoatephosphatidylcholine (DPePC). The Mulliken and RESP methods are denoted I and II, respectively. Also shown are the dipalmitoylphosphatidylcholine (DPPC) charges transferred from DPePC using transfer schemes A and B described in Our Parameterization Strategy. The CHARMM27 partial charges are shown for reference. Fig. 7 gives an overview of the new RESP charges listed in Table 2 after the transfer to DPPC. To increase clarity in the figure, the charges are shown with two decimals only. Due to round-offs, the molecules have a nonzero net charge, which is not the case when all six decimals are used (see Table 2).


   ACKNOWLEDGEMENTS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
REFERENCES
 
The authors thank John F. Nagle for valuable suggestions and for sharing unpublished results with us. We also gratefully acknowledge the constructive ideas of our reviewers as well as the helpful suggestions of Richard Pastor.

This work was supported by the Danish National Research Foundation via a grant to the MEMPHYS-Center for Biomembrane Physics.


   FOOTNOTES
 
Morten Ø. Jensen's present address is D. E. Shaw Research, New York, NY 10036.

Submitted on April 17, 2006; accepted for publication January 5, 2007.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 SUMMARY AND CONCLUSION
 APPENDIX: NEW PARTIAL CHARGES
 ACKNOWLEDGEMENTS
 REFERENCES
 
1. Feller, S. E., and R. W. Pastor. 1999. Constant surface tension simulations of lipid bilayers: the sensitivity of surface areas and compressibilities. J. Chem. Phys. 111:1281–1287.[CrossRef]

2. Jensen, M. Ø., O. G. Mouritzen, and G. H. Peters. 2004. Simulations of a membrane-anchored peptide: structure, dynamics, and influence on bilayer properties. Biophys. J. 86:3556–3575.[Abstract/Free Full Text]

3. Bayly, C. I., P. Cieplak, W. D. Cornell, and P. A. Kollman. 1993. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 97:10269–10280.[CrossRef]

4. Nagle, J. F., and S. Tristram-Nagle. 2000. Structure of lipid bilayers. Biochim. Biophys. Acta. 1469:159–195.[Medline]

5. Carrillo-Tripp, M., and S. E. Feller. 2005. Evidence for a mechanism by which {omega}-3 polyunsaturated lipids may affect membrane protein function. Biochemistry. 44:10164–10169.[CrossRef][Medline]

6. Bandyopadhyay, S., J. C. Shelley, and M. L. Klein. 2001. Molecular dynamics study of the effect of surfactant on a biomembrane. J. Phys. Chem. B. 105:5979–5986.

7. Ash, W. L., M. R. Zlomislic, E. O. Oloo, and D. P. Tieleman. 2004. Computer simulations of membrane proteins. Biochim. Biophys. Acta Biomembr. 1666:158–189.[Medline]

8. MacCallum, J. L., and D. P. Tieleman. 2006. Computer simulation of the distribution of hexane in a lipid bilayer: spatially resolved free energy, entropy, and enthalpy profiles. J. Am. Chem. Soc. 128:125–130.[CrossRef][Medline]

9. Callisen, T. H., and Y. Talmon. 1998. Direct imaging by Cryo-TEM shows membrane break-up by phospholipase A2 enzymatic activity. Biochemistry. 37:10987–10993.[CrossRef][Medline]

10. Andresen, T. L., S. S. Jensen, and K. Jørgensen. 2005. Advanced strategies in liposomal cancer therapy: problems and prospects of active and tumor specific drug release. Prog. Lipid Res. 44:68–97.[CrossRef][Medline]

11. Sukharev, S. I., P. Blount, B. Martinac, F. R. Blattner, and C. Kung. 1994. A large-conductance mechanosensitive channel in E. coli encoded by MscL alone. Nature. 368:265–268.[CrossRef][Medline]

12. Colombo, G., S. J. Marrink, and A. E. Mark. 2003. Simulation of MscL gating in a bilayer under stress. Biophys. J. 84:2331–2337.[Abstract/Free Full Text]

13. Gullingsrud, J., and K. Schulten. 2004. Lipid bilayer pressure profiles and mechanosensitive channel gating. Biophys. J. 86:3496–3509.[Abstract/Free Full Text]

14. James, S. R., R. A. Demel, and C. P. Downes. 1994. Interfacial hydrolysis of phosphatidylinositol 4-phosphate and phosphatidylinositol 4,5-bisphosphate by turkey erythrocyte phospholipase C. Biochem. J. 298:499–506.[Medline]

15. Cantor, R. S. 1997. The lateral pressure profile in membranes: a physical mechanism of general anesthesia. Biochemistry. 36:2339–2344.[CrossRef][Medline]

16. Cantor, R. S. 2003. Receptor desensitization by neurotransmitters in membranes: are neurotransmitters the endogenous anesthetics? Biochemistry. 42:11891–11897.[CrossRef][Medline]

17. Jensen, M. Ø., and O. G. Mouritsen. 2004. Lipids do influence protein function—the hydrophobic matching hypothesis revisited. Biochim. Biophys. Acta Biomembr. 1666:205–226.[Medline]

18. Anézo, C., A. H. de Vries, H. Holtje, D. P. Tieleman, and S. Marrink. 2003. Methodological issues in lipid bilayer simulations. J. Phys. Chem. B. 107:9424–9433.

19. Wohlert, J., and O. Edholm. 2004. The range and shielding of dipole-dipole interactions in phospholipid bilayers. Biophys. J. 87:2433–2445.[Abstract/Free Full Text]

20. Patra, M., M. Karttunen, M. T. Hyvonen, E. Falck, P. Lindqvist, and I. Vattulainen. 2003. Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interactions. Biophys. J. 84:3636–3645.[Abstract/Free Full Text]

21. Sum, A., R. Faller, and J. de Pablo. 2003. Molecular simulation study of phospholipid bilayers and insights of the interactions with disaccharides. Biophys. J. 85:2830–2844.[Abstract/Free Full Text]

22. Berger, O., O. Edholm, and F. Jahnig. 1997. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 72:2002–2013.[Abstract/Free Full Text]

23. Jähnig, F. 1996. What is the surface tension of a lipid bilayer membrane? Biophys. J. 71:1348–1349.[Free Full Text]

24. Marrink, S. J., and A. E. Mark. 2001. Effect of undulations on surface tension in simulated bilayers. J. Phys. Chem. B. 105:6122–6127.

25. MacKerell, A. D., Jr., D. Bashford, M. Bellott, R. L. Dunbrack, Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, and M. Karplus. 1998. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B. 102:3586–3616.

26. Schlenkrich, M., J. Brickmann, A. D. MacKerell, Jr., and M. Karplus. 1996. Empirical potential energy function for phospholipids: criteria for parameter optimization and applications. In Biological Membranes: A Molecular Perspective from Computation and Experiment. K. Merz, Jr. and B. Roux, editors. Birkhäuser, Boston, MA.

27. Foloppe, N., and A. D. MacKerell, Jr. 2000. All-atom empirical force field for nucleic acids. I. Parameter optimization based on small molecule and condensed phase macromolecular target data. J. Comput. Chem. 21:86–104.[CrossRef]

28. Feller, S. E., and A. D. MacKerell, Jr. 2000. An improved empirical potential energy function for molecular simulations of phospholipids. J. Phys. Chem. B. 104:7510–7515.

29. Nagle, J. F., R. Zhang, S. Tristram-Nagle, W. Sun, H. I. Petrache, and R. M. Suter. 1996. X-ray structure determination of fully hydrated L{alpha} phase dipalmitoylphosphatidylcholine bilayers. Biophys. J. 70:1419–1431.[Abstract/Free Full Text]

30. 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 alkanes and its effect on molecular simulations of alkanes and a DPPC bilayer. J. Phys. Chem. B. 109:5300–5311.[Medline]

31. Skibinsky, A., R. M. Venable, and R. W. Pastor. 2005. Muscle and contractility—a molecular dynamics study of the response of lipid bilayers and monolayers to trehalose. Biophys. J. 89:4111–4121.[Abstract/Free Full Text]

32. Klauda, J. B., N. Kuçerka, B. R. Brooks, R. W. Pastor, and J. F. Nagle. 2006. Simulation-based methods for interpreting x-ray data from lipid bilayers. Biophys. J. 90:2796–2807.[Abstract/Free Full Text]

33. Feller, S. E., and R. W. Pastor. 1996. On simulating lipid bilayers with an applied surface tension: periodic boundary conditions and undulations. Biophys. J. 71:1350–1355.[Abstract/Free Full Text]

34. Benz, R. W., F. Castro-Roman, D. J. Tobias, and S. H. White. 2005. Experimental validation of molecular dynamics simulations of lipid bilayers: a new approach. Biophys. J. 88:805–817.[Abstract/Free Full Text]

35. Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. 1983. CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4:187–217.[CrossRef]

36. MacKerell, A. D., Jr., B. Brooks, C. L. Brooks, III, L. Nilsson, B. Roux, Y. Won, and M. Karplus. 1998. CHARMM: the energy function and its parameterization with an overview of the program. In Encyclopedia of Computational Chemistry, Vol. 1. P. v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F. Schaefer, III, and P. R. S. Schreiner, editors. John Wiley & Sons, Chichester, UK.

37. Lindahl, E., and O. Edholm. 2000. Spatial and energetic-entropic decomposition of surface tension in lipid bilayers from molecular dynamics simulations. J. Chem. Phys. 113:3882–3893.[CrossRef]

38. Sonne, J., F. Y. Hansen, and G. H. Peters. 2005. Methodological problems in pressure profile calculations for lipid bilayers. J. Chem. Phys. 122:124903.[CrossRef][Medline]

39. Jensen, T. R., M. Ø. Jensen, N. Reitzel, K. Balashev, G. H. Peters, K. Kjaer, and T. Bjørnholm. 2003. Water in contact with extended hydrophobic surfaces: direct evidence of weak dewetting. Phys. Rev. Lett. 90:086101/1–086101/4.

40. Jensen, M. Ø., O. G. Mouritsen, and G. H. Peters. 2004. The hydrophobic effect: molecular dynamics simulations of water confined between extended hydrophobic and hydrophilic surfaces. J. Chem. Phys. 120:9729–9744.[CrossRef][Medline]

41. Singh, U. C., and P. A. Kollman. 1984. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 5:129–145.[CrossRef]

42. Woods, R. J., and R. Chappelle. 2000. Restrained electrostatic potential atomic partial charges for condensed-phase simulations of carbohydrates. J. Mol. Struct. THEOCHEM. 527:149–156.[CrossRef]

43. Reynolds, C. A., J. W. Essex, and W. G. Richards. 1992. Atomic charges for variable molecular conformations. J. Am. Chem. Soc. 114:9075–9079.[CrossRef]

44. Kucerka, N., S. Tristram-Nagle, and J. Nagle. 2006. Closer look at structure of fully hydrated fluid phase DPPC bilayers. Biophys. J. 90:L83.[Abstract/Free Full Text]

45. Feller, S. E., R. M. Venable, and R. W. Pastor. 1997. Computer simulation of a DPPC phospholipid bilayer: Structural changes as a function of molecular surface area.