A Molecular Dynamics Investigation of the Influence of Hydration and

The interaction of local anesthetics with lipid membranes. Patricio A. Zapata-Morin , F.J. Sierra-Valdez , J.C. Ruiz-Suárez. Journal of Molecular Gra...
0 downloads 0 Views 441KB Size
Download PDF
14326

J. Phys. Chem. B 2006, 110, 14326-14336

A Molecular Dynamics Investigation of the Influence of Hydration and Temperature on Structural and Dynamical Properties of a Dimyristoylphosphatidylcholine Bilayer Carl-Johan Ho1 gberg and Alexander P. Lyubartsev* DiVision of Physical Chemistry, Arrhenius Laboratory, Stockholm UniVersity, S-106 91 Stockholm, Sweden ReceiVed: March 10, 2006; In Final Form: May 30, 2006

A series of molecular dynamics simulations of dimyristoylphosphatidylcholine bilayers, with different levels of hydration and temperature, were performed to examine the influence of hydration on properties of lipid membranes. Structural and dynamical properties such as area per lipid, electron densities, order parameters for all CH bonds and water, diffusion, and reorientation autocorrelation functions were determined and were all found to be affected by changes in the hydration level. The simulations give an overall picture of the bilayer going to a more ordered state when the hydration level is reduced. Lipid headgroups were found to adopt an orientation more parallel to the membrane plane when the water content was decreased. Dynamical properties such as lipid diffusion and relaxation of reorientation time correlation functions were found to become slower with the removal of water. Our simulation results generally agree with experimental data in cases where such data are available. One important conclusion drawn is that while structural properties are affected only moderately dynamical properties are affected very strongly by a decrease of water content.

Introduction Understanding of structural and dynamical properties of lipid membranes on the molecular level is a matter of fundamental interest. Lipid membranes provide the basic structural unit of all living cells, and their molecular properties have a decisive influence on the cell functioning. Not surprisingly, as soon as the molecular dynamics technique had been established, computer simulations of lipid bilayers became an important source of information on the behavior of lipid molecules in membranes and their interaction with the environment1-7 (for a recent review, see ref 8). During the last few years, new versions of the CHARMM9 and GROMOS10 force fields have appeared that are specially tuned for phospholipid simulations. Also, the size of simulated systems and simulation time have increased from tens of lipids and hundreds picoseconds in earlier works to hundreds of lipids and tens or even hundreds of nanoseconds in works of the last few years.11-16 In most of the existing molecular dynamics simulations of lipid membranes the amount of water has been about 20-30 molecules per lipid, which corresponds to a fully hydrated bilayer. There are however experimental data showing that structural, dynamical, and thermodynamical properties of membranes depend on the degree of hydration.17-19 A change of the degree of hydration influences also mechanical properties of a membrane and its phase behavior. A change of membrane properties caused by a change of hydration affects in turn membrane proteins and their biological functioning.20,21 Simulation data for lipid membranes at a low degree of hydration are rather fragmented.22 Armen and co-workers23 studied a number of phospholipids (dipalmitoylphosphatidylcholine (DPPC), dioleoylphosphatidylcholine (DOPC), and palmitoyloleoylphosphatidylcholine (POPC)), each type of lipid at its own level of hydration. In the work by Mashl et al.,24 different properties of DOPC bilayers were studied at different * Author to whom correspondence should be addressed. E-mail: sasha@ physc.su.se.

levels of hydration from 5 to 30 water molecules per lipid. The lengths of molecular dynamics runs were, however, limited to 1.5 ns, which may be considered too short in comparison with the several nanosecond slow fluctuation modes of conformational changes in bilayers observed in simulations of Lindahl and Edholm.11 It is known from the large amount of experimental data that a fully hydrated dimyristoylphosphatidylcholine (DMPC) bilayer transforms from the liquid-crystalline phase to the gel phase at about 25 °C. At a low degree of hydration, the transition temperature increases up to about 40 °C.25,26 We therefore run two series of MD simulations, one at 30 °C, for which transition to the gel phase is expected at hydration levels below 15% water (6-7 water molecules per lipid), and one at 50 °C, where no transition to the gel phase takes place. The hydration level was changed in these simulations from 5 to 29 waters per lipid. At lower temperature (30 °C), the lengths of the simulation runs were extended to 200 ns. Such long simulations are especially important in the vicinity of the phase transition region, that is at lower temperatures and low hydration, since in this case the diffusion and overall conformational changes are much slower. Our aim is to give a detailed overview of the computed structural and dynamical properties and to analyze how the most important molecular properties of the lipid bilayer, accessible from molecular dynamics simulations, are affected by the amount of water and temperature of the system. We would also like to say something about how the force field behaves at low hydration levels and examine if the force field reproduces properties of the DMPC lipid not yet examined at different levels of hydration. Simulated Systems and Computational Details Three different lipid bilayer systems, each consisting of 128 (64 × 2) DMPC lipids and 578, 1250, or 3655 water molecules, corresponding to approximately 5, 10, or 29 water molecules per lipid, were simulated at two temperatures, 30 and 50 °C.

10.1021/jp0614796 CCC: $33.50 © 2006 American Chemical Society Published on Web 07/06/2006

MD Investigation of a DMPC Bilayer

Figure 1. Dimyristoylphosphatedylcholene molecule. Names of atoms used in the text are shown.

They are referred to as low, medium, and full hydration correspondingly. A simulation with 2400 water molecules (19 H2O/lipid) was also carried out at 30 °C. Its results turned out to be almost identical with the corresponding case of 29 H2O/ lipid (with the exception of some dynamical properties), and therefore it was omitted in the discussion of the results. An additional simulation with 5 water molecules per lipid at 30 °C starting from a gel phase was carried out to control convergence. For DMPC lipids, we used the united-atom model for hydrocarbon groups. Previous studies of phospholipid bilayers (see references in the Introduction) have shown that both unitedatom and all-atom force fields yield generally similar results, while the united-atom model provides a more than 3-fold savings of computer time due to the lesser number of atoms. The lipid force field parameters for bonded and nonbonded interactions as well as atomic partial charges are based on the GROMOS force field,4 with modification by Berger et al.6 The DMPC molecular structure with the names of atoms referenced in the text is shown in Figure 1. For water the simple point charge (SPC) model was used.27 The starting configuration for the fully hydrated (3655 water molecules) lipid bilayer simulated at 30 °C was downloaded from the Web site of the Biocomputing group at the University of Calgary (http://moose.bio.ucalgary.ca). In this starting configuration the lipids were in the liquid-crystalline phase. For the bilayers with less water, all water molecules were first removed and then added to the desired amount with an equal number of water molecules on both sides of the bilayer. The starting water configuration for these simulations was in a crystal structure. Then energy minimization and a 1 ns preequilibration run with isotropic cell fluctuations have been performed. All simulations were preformed with the Gromacs simulation package, version 3.2.28 All bond lengths were kept constant at their equilibrium values using the LINCS algorithm.29 The time step was set to 5 fs. The use of a 5 fs time step for the unitedatom lipid model and rigid water has been examined and validated in previous publications.30 The temperature in the simulated systems was regulated separately for lipids and water with the Nose´-Hoover thermostat.31 The relaxation constant for the thermostat was set to 0.5 ps. The pressure was set to 1 atm and was regulated by the Parrinello-Rahman barostat32 with a relaxation constant of 2.0 ps anisotropically, allowing the x-, y-, and z-dimensions of the rectangular simulation box to change independently of each other. The long-range electrostatic forces were treated by the particle mesh Ewald (PME) summation algorithm.33 It is widely recognized that an accurate account of long-range electrostatic interactions is important in the simulation of lipid membrane systems where the observed properties result from a very delicate balance of different forces, in which the dipole-dipole interactions of the headgroups play a major role. The Ewald sum computes electrostatic interactions as coming from all periodic images of the simulation cell, so such a simulation setup corresponds to lipid bilayers arranged in a multilamellar stack, which is the case for many experimental studies. The question to what extent the periodicity introduced by the Ewald approach

J. Phys. Chem. B, Vol. 110, No. 29, 2006 14327 affects simulation results is a matter of ongoing discussions, and approaches based on a reaction field treatment of longrange electrostatics are also considered. Still, in a number of works14,30,34-36 it was demonstrated that the Ewald technique provides better stability of different membrane properties relative to the system size and better agreement with experimental data in comparison with simulations employing the simple electrostatic cutoff. In a methodological study30 it was argued that when the cell dimension exceeds the mean headgroup-headgroup distance many times, the long-range ordering artifacts are counteracted by stronger short-range fluctuations and do not make a noticeable influence on the results. In another methodological study34 it was shown that, while using PME, the convergence of different structural simulation results relative to the system size occurs already at 36 lipids per leaflet. Though the latter conclusion was made for constant volume simulations, it is reasonable to suggest that it holds for constant pressure simulations that sample similar configurations in a relatively limited range of volumes, with the caveat that constant pressure sampling requires a longer time to reach convergence. The results of ref 14 confirm a stable area per lipid in the range of 32-512 lipids per leaflet in constant pressure simulations. To justify using PME in simulations with different hydration levels, we note that in ref 30 possible sources of artifacts in the lattice sum techniques were determined as (i) the buildup of artificial lateral dipole moments interacting favorably with each other across the periodic lattice and (ii) the absence of interaction between charges separated by half of the box size along one of the lattice directions. As it was noted above, possible effects of these factors for a large enough simulation cell are small, and it seems plausible to suggest that these factors, both acting on a longer range in the lateral direction, remain essentially the same at lower hydration when the bilayer surfaces approach each other and the distance between them is small relative to the simulation cell size. The dynamical properties can be also affected by the treatment of the long-range electrostatics,30 and recent work36 argues in favor of using the PME approach for the computation of diffusion and time correlation functions. For Lennard-Jones interactions, a cutoff radius was set at 1.0 nm. The dispersion correction from the interactions outside the cutoff was included into energy and pressure. A neighbor list scheme was used with the list update frequency at every 10th time step. The eventual center of mass motion of the system was removed every 500 fs. Simulations at 30 °C were run 200 ns from the starting conditions described above. The last 100 ns were used for the trajectory analysis. Such a long equilibration time was chosen to ensure that all important structural parameters have reached equilibrium values, taking in mind that the system may be close to the phase transition region. Simulations at 50 °C were run during 50 ns. Their starting configurations were taken from the final configurations of the corresponding simulations performed at 30 °C. As a consequence of this and due to higher temperature, the equilibrium time was set to 5 ns. Results and Discussion Area per Lipid. One of the most fundamental properties of a lipid bilayer and one of the most common ways to determine whether the bilayer system has reached equilibrium is the area per lipid. When the area per lipid reaches a stable value, most other structural properties of the lipid bilayer do not change either. Simulated area per lipid can also be compared with experimental values available from volumetric data and X-ray spacing.

14328 J. Phys. Chem. B, Vol. 110, No. 29, 2006

Ho¨gberg and Lyubartsev TABLE 2: Summary of Simulation Resultsa number of H2O/lipid 5

10

30 °C 50 °C

Area per Lipid (Å2) 57.9(0.2) 61.5(0.3) 60.7(0.2) 64.1(0.2)

63.1(0.3) 65.4(0.3)

30 °C 50 °C

d-Spacing (Å) 42.3(0.2) 45.7(0.2) 41.6(0.2) 45.0(0.2)

62.9(0.3) 61.1(0.3)

30 °C 50 °C

Membrane Thickness (Å) from Total Electron Density 35.3 (0.25) 33.7 (0.3) 32.5 (0.3) 34.1 (0.25) 32.7 (0.3) 32.2 (0.3)

30 °C 50 °C

Membrane Thickness (Å) from Lipid Electron Density 29.0 (0.25) 27.3 (0.3) 27.1 (0.3) 28.4 (0.25) 26.9 (0.3) 26.3 (0.3)

30 °C 50 °C

Headgroup Tilt Angle (deg) 85.1 (0.4) 81.9 (0.4) 83.8 (0.4) 81.1 (0.4)

DMPC Diffusion Coefficient (10-8cm2/s) 30 °C corrected 0.7(0.1) 3.2(0.3) 30 °C uncorrected 1.0 (0.2) 6.0(0.3) 50 °C corrected 1.5(0.3) 7.3(0.8) 50 °C uncorrected 2.6 (0.3) 11.8(1) H2O Diffusion Coefficient (10-5cm2/s) 30 °C bulk ) 4.98 (0.04) 0.073(0.01) 0.594(0.01) 50 °C bulk ) 6.96 (0.03) 0.140(0.003) 0.85(0.01)

Figure 2. Area per lipid as a function of simulation time for the cases (a) T ) 303 K and (b) T ) 323 K: 29 H2O/lipid, black; 10 H2O/lipid, red; 5 H2O/lipid, blue.

TABLE 1: Subaverages of the Area per Lipid and Standard Deviations (in Å2) Computed from Different Parts of MD Trajectories for Simulations at 30 °C 29 H2O/lipid

10 H2O/lipid

5 H2O/lipid

time (ns)

〈A〉

〈δA〉

〈A〉

〈δA〉

〈A〉

〈δA〉

0-50 50-100 100-150 150-200

63.18 63.55 63.17 63.18

1.35 1.22 1.14 1.33

61.01 61.23 61.22 61.74

1.21 0.96 0.88 1.18

58.72 58.11 58.0 57.91

1.01 1.08 0.75 0.63

The time evolution of the area per lipid for the six simulations is presented in Figures 2a and 2b. The area per lipid seems to be stable in all the simulations. A more detailed statistical analysis is presented in Table 1 with average values of the area and its standard deviations computed over 50 ns windows of the trajectories. One can see that deviations of the average values are well within the standard deviations. One can note that for the lowest hydration standard deviation becomes noticeably smaller after about 100 ns of the simulation. Generally, the system with the lowest hydration shows overall slower dynamics (more detailed analysis is given below). To be on the safe side, we have chosen to compute average properties over the last 100 ns of the trajectories. Simulations carried out at a higher temperature (50 °C), which were started from the last configurations of the corresponding lower-temperature simulations, do not show any noticeable trend in behavior of the area per lipid from the beginning. Averaging of these trajectories was started after 5 ns. Note also that in both cases there are clearly seen fluctuations of the area per lipid on the time scale of several nanoseconds, which confirms results of earlier work11 and which

29

79.3 (0.4) 80.9 (0.4) 5.8 (0.3) 22 (1) 18.0 (2) 32.4 (2) 2.25(0.02) 3.09(0.03)

30 °C 50 °C

Lipid Relaxation Time (ns) 128 30.1 47.7 15.5

28.8 14.5

30 °C 50 °C

Headgroup Relaxation Time (ns) 103 13.7 27.0 12.3

11.7 9.7

30 °C bulk ) 1.39 50 °C bulk ) 1.10 a

H2O Relaxation Time (ps) 71.8 (0.2) 14.8 (0.06) 29.9 (0.15) 7.98 (0.03)

5.06 (0.02) 3.03 (0.01)

See text for details. Statistical uncertainty is shown in parentheses.

shows the necessity of long, many tens of nanoseconds, simulations of membranes to obtain reliable equilibrium properties. The average area per lipid for the production part of each simulation is presented in Table 2. The statistical uncertainty was obtained by dividing the production part of the trajectories into n windows, computing the average value in each window, taking the standard deviation of this subaverage, and dividing it by xn-1. Values n ) 10 and 20 were used for computation of statistical uncertainty of the area per lipid yielding similar results. The same procedure was used for evaluation of the statistical error of other properties. For the fully hydrated systems the area per lipid obtained was 63.1 Å2 for 30 °C and 65.4 Å2 for 50 °C. The result for the higher temperature perfectly coincides with experimental results (65.4 Å2) at the same temperature.37 For the lower temperature our result is somewhat higher than the experimental one for the same work (60 Å2). There exist other experimental data for the average area per lipid ranging between 59 and 65 Å2 (refs 38 and 39 and references therein), with a recent more accurate value of 60.6 ( 0.5 Å2 (ref 40) which is 2.5 Å2 lower than our simulated result. From Table 2 it is clear that the area per lipid decreases when the water content decreases. This happens at the both simulated temperatures. The change is larger when the water content is decreasing from 10 to 5 H2O/lipid in comparison to the change between 29 and 10 H2O/lipid. The explanation for this is that for a fully hydrated bilayer the water molecules penetrate the bilayer down to approximately the lipid ester oxygen. The water

MD Investigation of a DMPC Bilayer

J. Phys. Chem. B, Vol. 110, No. 29, 2006 14329

molecules interact with different parts of the lipid headgroup and make the bilayer swell (see also density distributions of water and lipids below). On the initial stage of dehydration, water molecules are removed mostly from the water layer between the bilayers. When the water content falls further, even water molecules in the headgroup layer begin to disappear. If most of the water molecules are removed, then the remaining molecules have less capacity to keep the lipids apart, and the distance between the lipids decreases. Even screening of the charged fragments of the lipid headgroups by water molecules is decreasing upon dehydration which leads to stronger attraction between neighboring headgroups which in turn favors lower area per lipid. The effect of decreasing area per lipid at lower hydration has been observed experimentally in osmotic stress measurements41 of lecithin bilayers where the reduction of water activity led not only to a decrease of the bilayer separation but also to a decrease of the area per lipid. For a given degree of hydration, the system adopts a conformation corresponding to the minimum of free energy upon variation of these two structural parameters, the bilayer spacing and area. Computations of the chemical potential of water in a bilayer system at different hydrations and free or fixed areas could give additional insight into this effect. An interesting observation for the different simulations is the difference in fluctuations of the area over time. For the bilayer at full hydration and 30 °C, the standard deviation for the area per lipid obtained from the last 100 ns of the simulations is 1.23 Å2 compared to 1.03 and 0.69 Å2 for the medium and low hydrations, respectively. The fluctuations of the total area 〈δA〉 are related to the compressibility modulus KA of the system12

KA )

kTA 〈δA〉2

(1)

where A is the average total area of the simulated bilayer. Compressibilities computed according to eq 1 are 272, 390, and 860 mN/m for full, medium, and low hydrations, respectively. This may be compared with the experimental value KA ) 108 ( 35 mN/m44 for full hydration at T ) 30 °C. Note that the undulation contribution to the compressibility is practically absent in our simulation due to relatively small system size. Simulations of larger fragments of a glycerolmonoolein (GMO) bilayer12 show that undulations can decrease the compressibility modulus by a factor of 2-3, which brings our result close to the experimental one. Note also that even after averaging over 100 ns of trajectories, uncertainty in evaluation of fluctuations of the area per lipid is rather high, about 0.1-0.15 Å2, leading to 20-30% uncertainty in the compressibility. Still, our data allow us to conclude that dehydration leads to smaller fluctuations of the area and correspondingly to a higher compressibility modulus, which reflects the more ordered properties of the less hydrated bilayers. Given the smaller area per lipid for these systems, the lipids have less translational freedom resulting in a more stable area per lipid. To ensure convergence of the simulation performed at low hydration and 303 K, which may correspond to gel-phase conditions, a control simulation was made starting from the gel phase. The gel phase was prepared by a simulation with the CHARMM27 force field9 which resulted in a clear gel phase with 49 Å2 area per lipid and about 30° tilt of lipid tails. (Our preliminary simulations with the CHARMM27 force field showed its preference for the gel phase. These simulations are in progress now, and results will be presented elsewhere.) Then

Figure 3. Area per lipid as a function of simulation time for 5 H2O/ lipid starting from the gel phase at 303 K.

hydrogens on the lipids were removed, the configuration was preequilibrated with the GROMOS force field parameters during 5 ns at a fixed volume, and then constant pressure simulations were switched on. The evolution of the area per lipid of the latter simulation is shown in Figure 3. One can see that after about 5 ns it gives about the same area per lipid (and the same other properties) as the simulation started from the liquidcrystalline phase. This fact demonstrates that the results presented below are, from a simulation point of view, calculated in well-equilibrated systems. According to the phase diagram26 for the DMPC/water mixture, at this temperature the low hydration system should be in the gel phase or in some mixed phase combining properties of both. Recently a ripple phase consisting of liquid-crystalline and gel-phase domains was observed in simulations of a DPPC bilayer with the GROMOS force field upon cooling at full hydration.42 Our low hydration simulation seems to correspond to a mixed phase structure that may be different from that described in ref 42. Thus, an analysis of the configurations did not show any presence of separate domains, or any tilt of the lipid tails typical for the gel phase. Electron Densities. Electron density, defined as charge per nm3, gives information about the structure along the bilayer normal. In molecular dynamics simulations, electron density is formed from the atom densities multiplied by the corresponding atom number and corrected for the partial atom charges. The electron density is related by the Fourier transform to the structure factor measurable in X-ray experiments and may serve as one of the basic criteria for validation of molecular dynamics simulations.43 The electron densities for the simulations with the highest and lowest hydrations are presented in Figure 4. We display separately the lipid and water densities as well as the total density at T ) 30 °C. All the plots have a very characteristic appearance that agrees well with experiment40,44 and previous simulation studies,6 including both positioning and magnitudes of the characteristic features of the profiles. From the electron density profile the bilayer head-to-head distance, defined as the distance between the two maxima in the plot (DHH), can be estimated. It is usual practice to define DHH from the maxima of the total electron density. For the fully hydrated bilayer at 30 °C this yields the result DHH ) 32.5 Å, which is slightly less than experiment (35.3 Å, ref 40). As a consequence of this, the area per lipid in our simulations is slightly higher than the area reported in the same experimental work. An interesting fact is that if one considers electron density formed by the lipids only, then their maxima occur at a closer

14330 J. Phys. Chem. B, Vol. 110, No. 29, 2006

Figure 4. Electron density distribution of lipids and water along the z-axis. Total electron densities: green solid line for T ) 303 K, 29 H2O/lipid; green dashed line for T ) 303 K, 5 H2O/lipid. Lipid electron densities: solid black line for T ) 303 K, 29 H2O/lipid; black dashed line for T ) 323 K, 29 H2O/lipid; solid blue line for T ) 303 K, 5 H2O/lipid; blue dashed line for T ) 323 K, 5 H2O/lipid. Water electron densities: black dotted line for T ) 303 K, 29 H2O/lipid; black dashdotted line for T ) 323 K, 29 H2O/lipid; blue dotted line for T ) 303 K, 5 H2O/lipid; blue dash-dotted line for T ) 323 K, 5 H2O/lipid.

distance to the membrane mid-plane than the maxima coming from the total density; see Figure 4. The difference in maxima positions comes from the strongly varying water density in the headgroup region and leads to a 5-6 Å difference between the apparent head-to-head distance (defined by the total electron density) and the true average head-to-head distance (defined by the lipid electron density). For convenience, we give results obtained by both definitions in Table 2. Different definitions of the average head-to-head distances may have some implications for the determination of the average area per lipid from X-ray data since DHH is often used as a parameter for volumetric partitioning between “water” and “lipids”.44 The head-to-head distance defined from the maximum of the lipid electron density corresponds to the average distance between centers of the headgroups. To evaluate the volume occupied by lipids, one needs to add the effective van der Waals diameter to these DHH values. Assuming 5-6 Å as a reasonable value for the van der Waals diameter, we obtain an effective membrane thickness about equal to that determined from the maxima of the total electron density, which thus seems to be a fair evaluation of the membrane thickness. Upon decrease of the water content in the bilayer system, the bilayer becomes thicker. In the present study the effect is clearly visible only when the water content is 5 H2O/lipid. This can be understood from the fact that when the lipids are dehydrated they move to a more ordered conformation indicated by the smaller area per lipid. To be able to fit in a smaller area, the tails of the lipids have to adopt a more extended conformation, which causes the membrane to be thicker. For the more hydrated systems the lipid headgroups are still surrounded by enough water, and the effect is not seen. Another observation is that the minimum of electron density around the bilayer center becomes more pronounced when the bilayer is dehydrated. The minimum is called a lipid trough and can be prescribed to the lower electron density of methyl groups at the ends of the lipid tails. The stronger minimum of the electron density can be understood as a consequence of the more ordered state of lipids at a low degree of hydration. Also, it may mean that protrusions of lipids of one layer into another decrease at low hydration.

Ho¨gberg and Lyubartsev The electron density profile for water is also displayed in Figure 4. One can see that it is becoming broader at low hydration indicating a thicker bilayer when the water content is decreased. In Figure 4 it can also be seen that the electron density for the headgroup increases while the water content decreases. This can be interpreted in terms of the headgroup orienting itself to the bilayer x-y-plane. Without the water to interact with, this is a favorable position for the polar headgroup. This is confirmed also by observing separate density profiles for the headgroup nitrogen and phosphorus (data not shown) and by an increase of the headgroup tilt angle (see next section). The distance between the peaks for nitrogen and phosphorus is decreasing when removing water indicating that the two atoms are closer to each other with respect to the z-axis and the headgroup is more parallel to the bilayer x-y-plane. Distribution maxima from the phosphate and choline groups overlap more strongly with each other at low water content, resulting in the higher maxima of the electron density coming from the polar groups in Figure 4. Headgroup Tilt Angle. The results for the angle between the headgroup P-N vector, pointing from the phosphorus to the nitrogen atom, and the bilayer normal are presented in Table 2. The simulated results show that for the fully hydrated system the angles are 79.3° and 80.9° for 30 and 50 °C, respectively, indicating headgroups lying essentially flat down the bilayer plane. The headgroups are orienting themselves to the bilayer x-y-plane even more with decreasing hydration, the effect being more pronounced for the lower temperature. This behavior of the headgroup is in correspondence with previous experimental studies of other phosphocholine lipids.45 Previous simulation of the fully hydrated DPPC bilayer46 reported a result of 80° for the angle between the headgroup P-N vector and the bilayer normal, which is similar to our result for DMPC. Simulations carried out for some other lipids gave results of 60° for phosphatidylcholines47 and 68° for phosphatidylethanolamine.48 The change in the headgroup tilt angle may be related to the change in the area per lipid. The perpendicular to the membrane components of the headgroup dipoles repel each other, and though this repulsion is relatively small, it is of long-range character and creates a noticeable contribution to the balance of forces regulating the average area per lipid.14 When at low hydration the headgroups are oriented more in the x-y-plane; the perpendicular dipole component, and thus the repulsion force, are decreasing, favoring a smaller area per lipid. Note however that this relationship does not always hold upon the change of temperature: At full hydration, simulation at 50 °C shows a larger tilt angle and a larger area per lipid. Order Parameters. Deuterium order parameters for the lipid headgroup, glycerol moiety, and tails were calculated. The deuterium order parameter is related to the orientation of the principal frame of the electric field gradient around the deuterium nucleus substituting one of the hydrogens. It is directly measurable in NMR experiments. It is defined as

1 Si ) 〈3 cos2(θi) - 1〉 2

(2)

where θi is the angle between the ith C-H vector and the bilayer normal. The broken brackets denote an average over all molecules and time. In our simulations we do not have explicit hydrogens. It is however common practice to assume that C-H vectors for the ith carbon are perpendicular to the vector connecting the i - 1 and i + 1 carbons and that the i - 1, i + 1 carbons, and the two hydrogens adopt a tetrahedral conforma-

MD Investigation of a DMPC Bilayer

J. Phys. Chem. B, Vol. 110, No. 29, 2006 14331

TABLE 3: Order Parameters for the Lipid Headgroup and Glycerol Moietya experimental 5 H2O/lipid 10 H2O/lipid 29 H2O/lipid

a

b

R 30 °C 50 °C

0.062 0.068

0.058 0.057

30 °C 50 °C

0.119 0.102

0.106 0.098

30 °C 50 °C

-0.044 -0.044

-0.012 -0.022

30 °C 50 °C

0.215 0.247

0.239 0.234

0.250 0.233

30 °C 50 °C

-0.050 -0.033

-0.022 -0.023

G1 S -0.023 -0.024

-0.103/0.152*

-0.15

30 °C 50 °C

-0.214 -0.189

-0.1925 -0.187

G2 -0.191 -0.187

-0.143/0.212*

-0.20

30 °C 50 °C

-0.183 -0.188

-0.204 -0.189

G3 R -0.193 -0.191

0.156/0.212*

30 °C 50 °C

-0.243 -0.250

-0.271 -0.268

G3 S -0.273 -0.261

-0.167/0.224*

0.039 0.049 β 0.091 0.089 γ

-0.006 -0.008

0.035/0.050

0.04

-0.025/-0.042 -0.03

0.006/0.010

G1 R 0.005/-

0

-

-0.23

a Experimental order parameters are taken from ref 49 (given under column a and measured for bicelles/liposomes) and from ref 50 (column b) for the case of full hydration. An asterisk indicates that only the absolute value was measured. Statistical error for all values is within 0.005.

tion. The glycerol carbons g1 and g3 have two hydrogens that are distinguishable relative to the “CO” bond in the beginning of the second tail. They are also distinguishable in NMR order parameter measurements and conventionally labeled by the letters “r “and “s”.49 In Table 3 the order parameters for the R-, β-, and γ-headgroup carbons as well as glycerol carbons are presented (see Figure 1 for the definitions). The order parameters for the tails are given in Figures 5a and 5b. We have also presented the experimental NMR order parameters available in the literature for the case of full hydration. Our results are in good agreement with the experimental ones, except for the case of the g1 carbon, for which we obtained the order parameter for one of the hydrogens as 0.25 while experiments show values between -0.1 and -0.15. The positive value of the order parameter at g1 in the simulated system can be interpreted as the carbon chains are directed preferably along the z-axis (g2, g3, and the tail carbons) but make a step at the g1 carbon so that the initial piece of the longer lipid tail lies more in the membrane plane. Another noticeable deviation is for the β-carbon (0.09 and -0.03 for the simulated and experimental order parameters, respectively). When going from the full to low hydration, the order parameters for the headgroup and tails are increasing. This corresponds to a more ordered system. For some glycerol carbons, the order parameter is oppositely decreasing. The trends are however rather weak. Note also that the order parameter is not directly related to the degree of ordering in the system. It is rather a combination of “order” and “orientation”. For example, a zero order parameter may mean both a completely unordered (isotropic) system and a perfectly ordered one at the magic angle of 54.7° to the reference axis (eq 2). We have therefore

Figure 5. Deuterium order parameters given as absolute values for the DMPC tails, (a) first tail amd (b) second tail: 29 H2O/lipid, black; 10 H2O/lipid, red; 5 H2O/lipid, blue. Simulations at 303 K are marked by solid lines, and simulations at 323 K by dashed lines.

Figure 6. Distribution of angles between R, g1r, and tail1-4 CH vectors and the bilayer normal, for 303 K: 29 H2O/lipid and R carbon, black solid line; 5 H2O/lipid and R carbon, blue solid line; 29 H2O/ lipid and the g1 carbon, black dashed line; 5 H2O/lipid and the g1 carbon, blue dashed line; 29 H2O/lipid and the fourth carbon in tail one, black dotted line; 5 H2O/lipid and the fourth carbon in tail one, blue dotted line.

calculated the distribution of θ angles for some of the CH bonds and displayed them in Figure 6. One can see that distributions of all angles become slightly narrower at low hydration. The maxima of the distributions are

14332 J. Phys. Chem. B, Vol. 110, No. 29, 2006 also somewhat shifted. For the “g1r” CH bond, this leads to a decrease of the order parameter because the center of the distribution for cos θ is shifted from 0.8 to 0.77 (or from 36.5° to 39° for the angle itself). The order parameters for lipid tails are shown in Figures 5a and 5b. For the fully hydrated bilayer they are in a good agreement, with deviations of less than 0.02 from many experimental studies and earlier simulations.6,37,49,51 One can even note an “odd-even” effect with order parameters being slightly lower in the odd positions.52 When the water content is reduced, the order parameters for both tails are becoming more negative. It is seen from Figure 6 that this is caused by the narrowing of the distribution of the θ angle around 90°. The effect of hydration on the tail order parameters is also more pronounced than in the case of the lipid headgroup or glycerol moiety. This is in a qualitative agreement with recent experimental NMR studies of DMPC lipids at different levels of hydration,53 which demonstrated an even larger increase of the magnitude of the tail order parameter when going from full to low hydration (by about 0.08-0.1 in comparison with a 0.040.05 increase in our simulations). Since the order parameters for the tails are negative, the interpretation here is that the lipid chains adopt a more extended and ordered conformation when the amount of water is decreased. This is also in correspondence with the overall tendency of the system to adopt a more ordered gellike phase at low hydration. However the effect of ordering seems to be weaker in the case of simulations, which may originate in some preference of the force field for the fully hydrated liquid-crystalline phase. The temperature also affects the order parameters. Figures 5a and 5b demonstrate an approximately 0.02 decrease of the absolute values of the tail order parameters with an increase of the temperature from 30 to 50 °C. The effect is quite expected, less ordering with higher temperature. A similar effect was observed in experimental studies by Petrache et al.37 with a slightly higher temperature effect. A detailed comparison with the experimental data of ref 37 also shows that while for a temperature of 50 °C and full hydration the agreement between the simulated and experimental results is almost perfect, for 30 °C the simulated order parameters are about 0.01-0.02 less negative than the experimental ones. Figures 7a and 7b show the first rank order parameters of water (that is, average cosine of the angle between water dipole and membrane normal) as a function of the distance from the bilayer center. The functions have a clear maximum at about 16-17 Å from the membrane mid-plane. This distance corresponds to the lipid headgroup area seen in the electron density profiles. This is also the region of the strongest electrostatic field caused by the charges at N and P groups and which in turn causes preferential orientation of water dipoles. The magnitudes of maxima are slightly higher at lower temperature and lower hydration. Another interesting observation concerns the behavior of the water order parameter in the region between the two bilayers. On the graph, the curves in each case are drawn up to the box border, which is just in the middle of the water layer separating bilayer and its periodic image. In the cases of low and medium hydrations, these curves decrease almost linearly down to zero, and one can imagine that they continue with the same slope to the negative values after crossing the middle point of the water layer, indicating the opposite orientation of water molecules in the neighboring bilayer. Only in the case of high hydration (29 H2O/lipid) the water orientational order parameter is flattened at zero values, indicating the presence of “bulk” water. Such

Ho¨gberg and Lyubartsev

Figure 7. First rank order parameters of water dipoles as a function of distance to the bilayer center for (a) 303 K and (b) 323 K: 29 H2O/ lipid, black; 10 H2O/lipid, red; 5 H2O/lipid, blue.

“bulk” water disappears already at moderate hydration (10 water/ lipid), when all water molecules become affected by the lipids. Spatial Distribution Functions. To obtain a more detailed three-dimensional picture of the spatial arrangement of the lipid headgroups, spatial distribution functions (SDFs) for the headgroup nitrogen and the phosphorus were constructed with respect to each other. The local basis for the headgroup SDF was determined by the phosphorus atom, the nitrogen atom, and the carbon atom in the junction of the two tails. The distributions (SDFs) are visualized by the surfaces drawn at the given intensity level, showing the ratio of local to average densities. The SDFs for nitrogen and phosphorus atoms in adjacent lipids are shown in Figure 8 for the cases of (a) 5 and (b) 29 water molecules per lipid. The SDFs are colored orange and purple for nitrogen and phosphorus, respectively. For the case of full hydration the spatial distribution function clearly indicates that the nitrogens have a large preference for the phosphate group and vice versa. The scenario is the same for the 5 water molecules per lipid simulation. The nitrogen and phosphorus groups are located on the both sides of the lipids with which they are interacting. This is in accordance with the experimentally observed and theoretically described phenomena of lipid charge pairing.54 For the simulation with the lowest water content another maximum of the spatial distribution for nitrogen can be seen from the top of the phosphate group. This contribution results from the periodic boundary conditions and the fact that when the water content in the system is low, the distance between

MD Investigation of a DMPC Bilayer

J. Phys. Chem. B, Vol. 110, No. 29, 2006 14333

Figure 9. Lipid mean square displacements for the determination of diffusion coefficients: 29 H2O/lipid, black; 10 H2O/lipid, red; 5 H2O/ lipid, blue. Simulations at 303 K are marked by solid lines, and simulations at 323 K by dashed lines.

Figure 8. Spatial distribution functions for headgroup nitrogen and headgroup phosphorus at T ) 303 K. The distribution of both P (purple) and N (orange) are plotted with a density level 15 (a) for 5 H2O per lipid and (b) for 29 H2O per lipid. (The SDF density level is the ratio of the local to the average densities of the corresponding atoms.)

adjacent membrane layers is strongly reduced, and we see a contribution to the nitrogen SDF from the nitrogens of the adjacent membrane bilayer. This is an interesting observation showing that at lower water content there is significant correlation between the headgroups of neighboring layers. Physically, these correlations are caused by stronger electrostatic interactions due to smaller separations between layers. Diffusion. To gain information about the diffusion of the lipids in the simulated bilayers, the diffusion coefficient was calculated from the Einstein relation

D ) lim tf∞

1 d 〈|∆rxy(t)|2〉 2n dt

(3)

where 〈|∆rxy(t)|2〉 ) 〈|rxy(t0 + t) - rxy(t0)|2〉 is the mean square displacement (MSD) in the x-y-plane during time t starting from initial time t0, n is the dimensionality of the system (here for the lateral diffusion n ) 2), and D is the diffusion coefficient. The average is taken over all the molecules and over all initial time moments t0 acceptable for the given time t. It was pointed out in a number of papers30,55 that, while calculating diffusion coefficients of lipids in a bilayer, one

should calculate the MSD of lipids relative to the center of mass motion of each monolayer separately, since monolayers can move relative to each other, maintaining a fixed total center of mass. The motion of the center of mass may increase the observed diffusion by a factor of 3-4, and since its contribution is artificial, it must be subtracted. The MSD time dependencies, obtained after compensation for the monolayer center of mass motion, are displayed in Figure 9 for the six simulations. Initial times t0 in eq 3 were taken in the range from 100 ns to 200 t ns in 30 °C simulations and from 5 ns to 50 - t ns in simulations at 50 °C. As can be seen from the figure, the dependencies can be fairly well fitted to straight lines at t > 5 ns yielding the diffusion coefficient. Deviations from the linear behavior seen at large t and 50 °C are due to low statistics of averaging over acceptable t0. Analogous dependencies for water show practically ideal linear behavior (data not shown). In Table 2 the values of the diffusion coefficients for both water and lipids are given. For convenience, lipid diffusion coefficients obtained without correction for the monolayer center of mass motion are also shown, though in the discussion below we refer to the corrected diffusion coefficients. Statistical uncertainty is evaluated from the fluctuations of the slopes of MSD versus time plots. For the fully hydrated system at 30 °C the diffusion coefficient is calculated to be 5.8 × 10-8 cm2/s. It is in a perfect agreement with the experimental result of 6 × 10-8 cm2/s obtained by fluorescence recovery measurements56 but somewhat below the experimental NMR result57 of 9 × 10-8 cm2/s at the same temperature. For lower hydration (20% weight H2O) another experimental result is available58 3 × 10-8 cm2/s, which is also in agreement with our result of 3.2 × 10-8 cm2/s for the simulation at 10 H2O/lipid and 30 °C. As can be seen from Table 2, further lowering of the water content in the bilayer causes the diffusion coefficient to decrease drastically, and for the system with the lowest water content the diffusion coefficient is about one order lower than that at the full hydration. The increase in temperature from 30 to 50 °C causes a noticeable increase of the diffusion coefficient. This increase is more pronounced at full hydration (a factor of 3.1) than those in the cases of medium (2.3) and full hydration (2.1). Assuming Arrhenius behavior of the diffusion coefficient, the corresponding activation energies for full, medium, and low hydration are 46, 34, and 30 kJ/M. This activation energy is in perfect agreement with work by Oradd et al.57 for full hydration (49 kJ/M) while for medium hydration it is about 2 times lower in

14334 J. Phys. Chem. B, Vol. 110, No. 29, 2006

Ho¨gberg and Lyubartsev

comparison with the results of Kuo et al.58 (63 kJ/M). Oppositely, Ameida et al.56 report a slower increase of diffusion, by factor 2.2 with temperature, changing from 30 to 50 °C at full hydration. One of the simple theoretical models for a qualitative description of diffusion is the free volume theory of Cohen and Turnbull,56,59 which recently attracted much interest for understanding diffusion in membranes and the interconnection of diffusion with structural properties.60,61 According to the free volume theory,56 a diffusing particle must find a vacancy (empty space or void) next to it and have enough energy to overcome interaction with its neighbors and jump to this vacancy. While the latter factor leads to Arrhenius-like behavior of the diffusion coefficient, the former one depends on packing effects and often is used to describe a change in diffusion based upon a change in the membrane composition, such as, for example, the addition of cholesterol.56,60 For the lateral diffusion in membranes the free volume theory predicts the factor exp(-a0/aF) in the diffusion coefficient, where a0 is the cross-sectional molecular area and aF is the free area per molecule (with a0 + aF giving the total area per molecule).56 Accurate evaluations of free area along the membrane normal made in works of Falck et al.60 and Kupainen et al.61 showed that area profiles (that is, average areas occupied by lipid or solvent) follow very closely the corresponding mass or electron density profiles. Since the size and the mass of water molecules are much less than those of lipids and water molecules move on much faster time scale, one can assume that the volume occupied by water is “free” for lipids. Our simulations show that electron densities of lipids increase noticeably upon dehydration in the region of headgroups (black and blue lines in Figure 4). Thus, at low hydration, less free area remains for the lipid headgroup to move in, which, according to the free area theory, should lead to slower diffusion. The strong dependence of lipid diffusion on the degree of hydration is related to the fluid character given to the hydrated bilayer by the water molecules. When the water is removed, the area per lipid decreases, imposing steric hindrances on the lipid molecules and causing decreased diffusion. Another factor that may have an effect is charge pairing as described above between the positive nitrogen and the negative phosphorus groups, including that between adjacent layers, which causes the lipids to be more locked in their positions. Rotational Dynamics. To obtain further insight into the dynamics of lipids, reorientatioinal time correlation functions (TCFs) of lipids and water were examined. The time correlation function is determined as

c(t) )

〈n b(t0 + t)n b(t0)〉 〈n b(t0)2〉

(4)

where b n(t) is some vector fixed in the molecular frame determining the orientation of the molecule at time t. As in the case of computation of the diffusion coefficients, the average is made over all molecules and all acceptable t for initial time moments t0 > 100 ns for 30 °C and t0 > 5 ns for 50 °C. For lipids, we calculated two different reorienational TCFs. For the first one the vector b n(t) was defined by the x-y projection of the phosphorus to the nitrogen vector. This TCF characterizes the dynamics of the lipid headgroups. The second TCF was defined by the x-y projection of the CH vector at the g2 glycerol carbon. The g2 carbon is a junction of the lipid headgroup and the two tails; that is why rotation of the corresponding CH bond may characterize rotation of the lipid as a whole. In both cases, the z-component of the vector was removed because it has a

Figure 10. Reorientational time correlation functions for headgroup rotation in the x-y-plane: 29 H2O/lipid, black; 10 H2O/lipid, red; 5 H2O/lipid, blue. Simulations at 303 K are marked by solid lines, and simulations at 323 K by dashed lines.

Figure 11. Reorientational time correlation functions for lipid (defined by the Cg2-H vector) in the x-y-plane: 29 H2O/lipid, black; 10 H2O/ lipid, red; 5 H2O/lipid, blue. Simulations at 303 K are marked by solid lines, and simulations at 323 K by dashed lines.

nonzero average value related to the fact that no “flip-flop” motions of lipids take place on the time scale of the simulations. The TCFs for the lipid headgroups are shown in Figure 10 and for the lipids as a whole in Figure 11. On a logarithmic scale, the curves have an almost linear character on times larger than 5 ns. The relaxation time, given in Table 2, was obtained by linear fitting of these curves for times of 5-15 ns. (For large times the functions become too noisy due to their low values and low statistics.) Analysis of the TCF shows that rotational dynamics of lipids was affected strongly by both the degree of hydration and the temperature. This is most strongly pronounced for the system with the lowest water content, 5 H2O/lipid and 30 °C, where the function shows a very high orientational correlation even after 40 ns. (Tracking of the TCF for a longer time period is difficult because of the low statistics.) A slower relaxation in this case can be attributed to the drastically hindered rotational motion of the headgroups when the water content is decreased. When the water is removed, there is an increased tendency for nitrogen-phosphorus pairing, seen in the SDFs above, forming a crystal-like structure. In this structure the headgroups have less orientational freedom resulting in a very slowly decaying TCF.

MD Investigation of a DMPC Bilayer

Figure 12. Reorientational time correlation functions of water: 29 H2O/lipid, black; 10 H2O/lipid, red; 5 H2O/lipid, blue. Simulations at 303 K are marked by solid lines, and simulations at 323 K by dashed lines.

It is interesting to note that rotational dynamics of the lipid headgroups in all cases presented in Table 2 are still faster than the rotation of the whole lipid. That is why slowing down of the overall rotation should not be attributed to the hindering rotations of the headgroup only. In fact, lower hydration and lower temperature lead also to decreases of the average area per lipid, which in turn hinder rotation in the tail area of lipids due to steric effects. For the 5 H2O/lipid system, both translational and rotational motions of lipids turned out to be extremely slow. Water Dynamics. Water self-diffusion coefficients are given in Table 2. They follow the same trends as the diffusion coefficient of lipids, slower diffusion at low water content and lower temperature. The decrease in the water content causes however much stronger reduction of water diffusion than that of lipids. The explanation of this fact is simple: While at full hydration there is a large fraction of water molecules moving almost freely as in bulk water, at low hydration almost all water molecules are bound to lipid headgroups that strongly retard their diffusion. The reorientational TCFs of water molecules are presented in Figure 12. The vector b n(t) defining the TCFs in eq 2 was determined as normal to the HOH plane. The water reorientational TCFs are decaying on a time scale of hundreds of picoseconds, with longer relaxation times for lower temperature and lower hydration. One can also note that the water relaxation time increases more at low water content (by a factor about of 10) than that of lipids (a factor of 3-5). NMR studies of water dynamics in DOPC and POPC bilayers at different hydrations62 also show about a 10-fold increase of water reorientation correlation time while going from full to about 5 H2O/lipid hydration. As in the case of diffusion, this fact is explained by the lower content of the “bulk” water at lower hydration. It is also worth noting that even in the case of the lowest hydration at 30 °C water dynamics (both translational diffusion and rotational relaxation) are much faster, by 2-3 orders of magnitude, than the dynamics of lipids. While the lipids are almost freezing, the water molecules are relatively freely moving between them. Conclusions In the present work, structural and dynamical properties of a DMPC lipid bilayer at different levels of hydration and temperatures have been studied. A summary of the most

J. Phys. Chem. B, Vol. 110, No. 29, 2006 14335 essential results is presented in Table 2. The area per lipid, electron density, lipid diffusion, and order parameters show good agreement with available experimental data. All examined properties were found to be affected by reduced water content and temperature. The simulations give an overall picture of the bilayer going to a more ordered state when the hydration level is reduced. This is reflected in the lower and more stable area per lipid, sharper features of the electron density distributions and angular distributions of selected CH vectors, and slower diffusion and reorientational dynamics. According to the phase diagram for DMPC, a bilayer at the lowest water content (5 H2O per lipid) and 30 °C should correspond to the β-gel phase.25 Our simulation shows some indications to approach the gel phase when the water content is reduced. This is clear from the area per lipid for the different simulations, behavior of the order parameters, and overall slower dynamics. However, no tendency to form a “tilt” structure of the tails, which is a typical feature of the gel phase, has been noted. Also the area per lipid is too high to correspond to the gel phase. This indicates that the results may correspond to some mixed phase between the liquid-crystalline and gel phase in the phase diagram.26 The observed “mixed” phase is however not of the type of the ripple phase consisting of fluidlike and gellike domains,42 since no such domains were observed in our simulations. Also, most of the effects that are observed at 30 °C are also seen in the 50 °C simulations, where the gel phase is not formed at any hydration level. That is why we still consider the simulation at 30 °C and low hydration as carried out in a liquid-crystalline phase and interpret all of the observed changes as caused by the change in hydration. The fact that the force field used yields results resembling the liquid-crystalline phase at 30 °C and 5 H2O/lipid hydration level, a somewhat higher area per lipid for the fully hydrated system, and somewhat smaller tail order parameters at 30 °C and their slower growth with dehydration may suggest that this force field has a preference for the liquid-crystalline phase. It is however a very difficult, often impossible, problem to calibrate a force field in such a way to reproduce all the properties within a few percent precision at varying thermodynamical conditions. Absence of explicit hydrogens in the GROMOS force field may also play a role, especially in light of a recent IR spectroscopy and Hartree-Fock study63 suggesting the hydrogen bonding of water with headgroup methyl and methylene moieties. Still, the trends observed in our simulations due to change in hydration and the temperature follow the same trends observed in the experimental studies. Our results support the previously made hypothesis that the lipid headgroups orient theirselves more into the x-y-plane of the lipid bilayer at lower water contents. A decrease of the normal to membrane component of the lipid dipole moments leads to a decrease of the electrostatic repulsion between them and favors a lower area per lipid. Another interesting observation, from both a theoretical and experimental point of view, is that strong correlations between the headgroup nitrogen and phosphorus atoms on adjacent leaflets arise when the water content is low. For experimentalists this should be an important point to take into consideration. Comparative analysis of the presented simulation results shows that while the structural properties of the bilayer only moderately depend on the degree of hydration (most of the observed changes are within 10%), the dynamics of the system are affected very strongly, showing sometimes more than 1 order of magnitude change in such parameters as the diffusion coefficient or reorientational correlation time.

14336 J. Phys. Chem. B, Vol. 110, No. 29, 2006 The hydration properties of bilayers should not only be of theoretical interest but also have important implications for cell membranes in vivo where on some occasions the water content is relatively low. This occurs, for example, in cell fusion and other biological relevant phenomena. Acknowledgment. We are grateful to Dick Sandstro¨m and Sergey Dvinskih for many discussions. This work has been supported by the Swedish Research Council (Vetenskapsrådet). We are also thankful to the Center for Parallel Computing at the Royal Institute of Technology and to the HPC2N Computer Center, Umeå, for granting computer facilities. References and Notes (1) Pastor, R. W. Curr. Opin. Struct. Biol. 1994, 4, 486-492. (2) Marrink, S. J.; Berendsen H. J. C. J. Phys. Chem. 1994, 98, 41554168. (3) Huang, P.; Perez, J. J.; Loew, G. H. J. Biomol. Struct. Dyn. 1994, 11, 927-956. (4) Tieleman, D. P.; Berendsen H. J. C. J. Chem. Phys. 1996, 105, 4871-4880. (5) Tieleman, D. P.; Marrink, S. J.; Berendsen, H. J. C. Biochim. Biophys. Acta 1997, 1331, 235-270. (6) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 20022013. (7) Zubrzycki, I. Z.; Xu, Y.; Madrid, M.; Tang, P J. Chem. Phys. 2000, 112, 3437-3441. (8) Scott, H. L. Curr. Opin. Struct. Biol. 2002, 12, 495-502. (9) Feller, S. E.; MacKerell, A. D. J. Phys. Chem. B 2000, 104, 75107515. (10) Chandrasekhar, I.; Kastenholz, M.; Lins, R. D.; Oostenbrink, C.; Schuler, L. D.; Tieleman, D. P.; van Gunsteren, W. F. Eur. Biophys. J. 2003, 32, 67-77. (11) Lindahl, E.; Edholm, O Biophys. J. 2000, 79, 426-433. (12) Marrink, S. J.; Mark, A. E. J. Phys. Chem. B 2001, 105, 61226127. (13) Pandit, S. A.; Bostick, D.; Berkowitz, M. L. Biophys. J. 2003, 85, 3120-3131. (14) Wohlert, J.; Edholm, O. Biophys. J. 2004, 87, 2433-2445. (15) Falck, E.; Patra, M.; Karttunen, M.; Hyvonen, M. T.; Vattulainen, I. Biophys. J. 2004, 87, 1076-1091. (16) Lopes, F. L.; Nielsen, S. O.; Klein, M. L. J. Phys. Chem. B 2004, 108, 6603-6610. (17) Ulrich, A. S.; Watts, A. Biophys. J. 1994, 66, 1441-1449. (18) Hristova, K.; White, S. H. Biophys. J. 1998, 74, 2419-2433. (19) Binder, H.; Gawrisch, K. J. Phys. Chem. B 2001, 105, 123781239. (20) Killian, J. A. Biochim Biophys. Acta 1998, 1376, 401-416. (21) Ueda, I.; Yoshida, T Chem. Phys. Lipids 1999, 101, 65-79. (22) Feller, S. E.; Yin, D.; Pastor, R. W.; MacKerell, A. D. Biophys. J. 1997, 73, 2269-2279. (23) Armen, R. S.; Uitto, O. D.; Feller, S. E. Biophys. J. 1998, 75, 734744. (24) Mashl, R. J.; Scott, H. L.,; Subramaniam S.; Jakobsson, E. Biophys. J. 2001, 81, 3005-3015. (25) Janiak, M. J.; Small, D. M.; Shipley, G. S. J. Biol. Chem. 1979, 254, 6068-6078. (26) Carlson, J. M.; Sethna, J. P. Phys. ReV. A 1987, 36, 3359-3374. (27) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Harmans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; pp 331-342. (28) Lindahl, E.; Hess, B.; van der Spoel. D. J. Mol. Model. 2001, 7, 306-317.

Ho¨gberg and Lyubartsev (29) Hess, B.; Bekker, H.; Berendssen, H. J. C.; Fraaije, J. G. E. J. Comput. Chem. 1997, 18, 1463-1472. (30) Anezo, C.; de Vries, A. H.; Holtje, H. D.; Tieleman, D. P.; Marrink, S. J. J. Phys. Chem. B 2003, 107, 9424-9433. (31) Hoover, W. G. Phys. ReV. A 1986, 34, 2499-2500. (32) Parinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196-1199. (33) Darden, T. A.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089-10092. (34) De Vries, A. H.; Chandrasekhar I.; van Gunsteren, W. F.; Hunenberger, P. H. J. Phys. Chem. B 2005, 109, 11643-11652. (35) Patra, M.; Karttunen, M.; Hyvo¨nen, M. T.; Falck., E., Lindqvist, P.; Vattulainen, I. Biophys. J. 2003, 84, 3636-3645. (36) Patra, M.; Karttunen, M., Hyvo¨nen, M. T.; Falck, E.; Vattulainen, I. J. Phys. Chem. B 2004, 108, 4485-4494. (37) Petrache, H.; Dodd, S.; Brown, M. Biophys. J. 2000, 79, 31723192. (38) Costigan, S. C.; Booth, P. J.; Templer, R. H. Biochim. Biophys. Acta 2000, 1468, 41-54. (39) Nagle, J. F.; Tristam-Nagle, S. Biochim. Biophys. Acta 2000, 1469, 159-195. (40) Kucerka, N.; Liu, Y.; Chu, N.; Petrache, H. I.; Tristam-Nagle, S.; Nagle, J. F Biophys. J. 2005, 88, 2626-2637. (41) Parsegian, V. A.; Fuller, N.; Rand, R. P. Proc. Natl. Acad. Sci. U.S.A. 1979, 76, 2750-2754. (42) de Vries, A. H.; Yefimov, S.; Mark, A. E.; Marrink, S. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 5392-5396. (43) Benz, R. W.; Castro-Roman, F.; Tobias, D. J.; White, S. H. Biophys. J. 2005, 88, 805-817. (44) Petrache, H. I.; Tristam-Nagle, S.; Nagle, J. F. Chem. Phys. Lipids 1998, 95, 83-94. (45) Bechinger, B.; Seelig, J. Chem. Phys. Lipids 1993, 58, 1-5. (46) A° man, K.; Lindahl, E.; Edholm, O.; Håkansson, P.; Westlund, P.O. Biophys. J. 2003, 84, 102-115. (47) Macdonald, P. M.; Leisen, J.; Marassi, F. M. Biochemistry 1991, 30, 3558-3566. (48) Raghavan, K.; Reddy, M. R.; Berkowitz, M. L. Langmuir 1992, 8, 233-240. (49) Aussenac, F.; Laguerre, M.; Schmitter, J.-M.; Dufourc, E. Langmuir 2003, 19, 10468-10479. (50) Gross, J. D.; Warschawski, D. E.; Griffin, R. G. J. Am. Chem. Soc. 1997, 119, 796-802. (51) Otten, D.; Brown, M. F.; Beyer, K. J. Phys. Chem. B 2000, 104, 12119-12129. (52) Douliez J.-P.; Leonadr, A.; Dufourc E. J. J. Phys. Chem. 1996, 100, 18450. (53) Dvinskikh, S.; Castro, V.; Sandstro¨m, D. Phys. Chem. Chem. Phys. 2005, 7, 3255-3257. (54) Pasenkiewicz-Gierula, M,; Takaoka, Y.; Miyagawa, H.; Kitamura, K.; Kusumi, A. Biophys J. 1999, 76, 1228-1240. (55) Lindahl, E.; Edholm, O. J. Chem. Phys. 2001, 115, 4938-4950. (56) Almeida, P. F. F.; Vaz, W. L. C.; Thompson, T. E. Biochemistry 1992, 31, 6739-6747. (57) Oradd, G.; Lindblom, G.; Westerman, P. W. Biophys. J. 2002, 83, 2702-2704. (58) Kuo, A.-L.; Wade, C. G. Biochemistry 1979, 18, 2300-2308. (59) Cohen, M. H.; Turnbull, D. J. Chem. Phys. 1959, 31, 1164-1169. (60) Falck, E.; Patra, M.; Karttunen, M.; Hyvo¨nen, M. T.; Vattulainen, I. Biophys. J. 2004, 87, 1076-1091. (61) Kupiainen, M.; Falck, E.; Ollila, S.; Niemela¨, P.; Gurtovenko, A. A.; Hyvo¨nen, M. T.; Patra, M.; Karttunen, M.; Vattulainen, I. J. Comput. Theor. Nanosci. 2005, 2, 401-413. (62) Volke, F.; Eisenbla¨tter, S.; Galle, J.; Klose, G. Chem. Phys. Lipids 1994, 70, 121-131. (63) Pohle, W.; Gauger, D. R.; Bohl, R.; Mrazkova, E.; Hobza, P. Biopolymers 2004, 74, 27-31.