Molecular Dynamics Study on Electrostatic Properties of a Lipid

Furthermore, to investigate an anisotropy of the rotational motion of water .... In Figure 13, a schematic picture of the direction of dipoles in the ...
0 downloads 0 Views 123KB Size
Download PDF
J. Phys. Chem. B 1998, 102, 6647-6654

6647

Molecular Dynamics Study on Electrostatic Properties of a Lipid Bilayer: Polarization, Electrostatic Potential, and the Effects on Structure and Dynamics of Water near the Interface Wataru Shinoda,† Masanori Shimizu, and Susumu Okazaki* Department of Electronic Chemistry, Tokyo Institute of Technology, 4259 Nagatsuta, Midori-ku, Yokohama 226-8502, Japan ReceiVed: March 12, 1998; In Final Form: June 3, 1998

Electrostatic properties of the dipalmitoylphosphatidylcholine bilayer in water have been investigated on the basis of molecular dynamics calculations in the isothermal-isobaric ensemble. Analysis of the long-time trajectory of partial charges on the atoms showed that polarization of the total system has an anisotropic fluctuation; the component along the bilayer normal was 1 order of magnitude smaller than that parallel to the bilayer plane, the polarization of lipid being canceled out mostly by that of water along the bilayer normal. A positive electrostatic potential was observed inside the membrane. The potential was found to vary in two steps across the interface; the first gradient comes from the preferential orientation of the ester group, and the second one is attributed to asymmetrical solvation of water around the phosphate group. Slower reorientation of the dipole projected on the bilayer normal was found than that on the bilayer plane. Dynamics of water near the headgroup was strongly influenced by the hydrogen bonding with the phosphate group.

I. Introduction Lipid molecules assemble spontaneously in water to form a bilayer with their headgroup facing toward water and with the tails of alkyl chains in contact with those in the counterpart layer. Then, charges on the headgroups produce an electric field, which is thought to play an essential role in the physical properties of the membrane.1 Further, the electric field much influences the permeability of polar small molecules and ions across the membrane as well as the conformation of proteins in the membrane. Thus, electrostatics of the membrane and its interface with water has been attracting much interest of both theorists and experimentalists.2 However, experimental approaches are, in general, not straightforward with respect to the microscopic electrostatics of the membrane, and available experimental data are insufficient to elucidate it. In this sense, computer simulation, which figures out structure and dynamics of the membrane directly on the atomic level, may be the most promising approach. In recent years, in fact, molecular dynamics (MD) simulation studies have been intensively done for several kinds of lipid bilayers.3-25 Analyses of the electric potential profile along the bilayer normal based upon the MD trajectories have also been made by several authors.3-9 In particular, a detailed analysis concerning the electric potential profile along the bilayer normal has been reported by Zhou and Schulten for dilaurylphosphatidylethanolamine (DLPE) bilayer based upon a 200 ps long stochastic boundary MD calculation.4 They demonstrated that the electric potential inside the bilayer is positive relative to the water layer and that it is caused mostly by dipoles of the ester group. However, it must be a little different for * To whom correspondence should be addressed: E-mail okazaki@ echem.titech.ac.jp; Fax +81-45-924-5489. † Present address: Yokohama Research Center, Mitsubishi Chemical Corp., 1000 Kamoshida-cho, Aoba-ku, Yokohama 227-8502, Japan. E-mail [email protected].

the case of dipalmitoylphosphatidylcholine(DPPC) bilayer of the present study, because DPPC has the bulky choline group instead of the ethanolamine. Further, except for the electric potential profile above, little has been clarified for other electrostatic properties of the system. For example, dielectric properties of the whole system have not been discussed yet despite its great importance in various model theories. In the present study, thus, dielectric properties, electric potential profile, and structure and dynamics of water at the interface have been investigated, in detail, for the DPPC bilayer in water. A rough analysis of water concentration dependence of the electrostatic properties is also reported. Five independent long-time MD simulations have been performed in order to obtain necessary statistics of the analysis. II. Calculations Five independent MD calculations of the DPPC bilayer system at two water concentrations have been performed in the liquid crystal phase at 353 K. The first three calculations are the calculations that have been reported in our previous paper;21 32 DPPC molecules together with 434 D2O molecules (13.5 water/lipid or 25 wt % water) were contained in the fully flexible cell in the periodic boundary condition. The calculations successfully reproduced structure and dynamics of the bilayer in the liquid crystal phase,21 although the size effect on the calculated physical properties is not clear at present. To examine the dynamics of water, 0.15 ns MD calculation of the DPPC bilayer in H2O has been performed at 25 wt % water, starting from one of the final configurations of the above three. In this calculation, a time step as short as 0.5 fs was needed to perform a stable MD run, whereas ∆t ) 1 fs was enough for the above D2O system. Further, to investigate the role of water concentration in the bilayer system, 0.6 ns long MD calculation for the membrane consisting of 32 DPPC and 1047 H2O molecules (32.7 water/

S1089-5647(98)01480-1 CCC: $15.00 © 1998 American Chemical Society Published on Web 08/05/1998

6648 J. Phys. Chem. B, Vol. 102, No. 34, 1998

Shinoda et al.

TABLE 1: Conditions of the MD Calculations system water concentration (wt %) simulation lengtha (ns) ∆t (fs) pressure control

32 DPPC + 434 D2O

32 DPPC + 434 H2O

32 DPPC + 1047 H2O

864 H2O

25

25

44

2.0 × 3

0.15

0.6

0.11

(1.5 × 3) 1.0 PRb

(0.1) 0.5 PRb

(0.5) 0.5 PRb

(0.1) 0.5 ANc

a The value in parentheses represents time length of the production run. b The Parrinello-Rahman method.27,28 c The Andersen method.30

Figure 1. Calculated electron density profiles of DPPC bilayer along the bilayer normal with the center of the bilayer at z ) 0. Thin solid line, 25 wt % water; dashed line, 44 wt % water; thick solid line, the experimental profile.35 The plot for the membrane at 25 wt % water is averaged over three MD calculations. Error bars in the figure for the 25 wt % water system represent the maximum and minimum values among three MD runs.

lipid) at 44 wt % water has been performed. This corresponds to the saturated or fully hydrated state. The initial configuration for this was constructed by adding the water slab, which had been equilibrated separately at the same thermodynamic condition, to the center of the water lamella of the final configuration of the 25 wt % water calculation. An MD calculation for the pure water, 864 H2O molecules, has also been done for 0.11 ns to obtain its dynamics as reference. All MD calculations were performed in the isothermalisobaric ensemble; the Nose´26-Parrinello-Rahman27-29 and Nose´-Andersen30 methods were used for the calculation of the bilayer and the pure water, respectively. Temperature of the thermal bath and the external pressure were set at 353 K and 0.1 MPa, respectively. The OPLS31 and AMBER32 force field were adopted for the DPPC molecule and water molecules were treated as the rigid body, the TIP4P model.31 The Gear predictor-corrector method33 was used to integrate the equations of motion. Long-range Coulombic interaction was calculated by the Ewald method.34 Lennard-Jones (LJ) interaction and the real-space Ewald terms were truncated with a cutoff radius of 11 Å. Correction of the LJ potential and virial with respect to the contribution from the outer region were also taken into account assuming random structure. MD trajectories were saved every 25 fs for the first three runs and every 10 fs for the others. The MD conditions are summarized in Table 1.

shows a satisfactory convergence and statistics of the calculations. No difference was found, of course, in the calculated structure of the membrane between the two kinds of solvent, i.e., H2O and D2O. Further, little difference is observed within the present calculations between the plots for the system at 25 and 44 wt % water concentrations, although the simulation time of 0.6 ns for the latter concentration might be too short to test the concentration dependence. This suggests that the structure of the membrane is almost independent of the water concentration, corresponding to a recent X-ray diffraction study35 where little structural change was measured for the DPPC bilayer in the liquid crystal phase at several lamella repeat spacings. The experimental electron density profile is also plotted in the figure. The MD calculations are in satisfactory agreement with the experiment. In particular, a good agreement is found for the headgroup region, i.e., |z| ) 17 ∼ 22 Å, where the electrostatic properties of the system are mostly determined. Peak to peak distance of the profile was calculated to be 38.9 ( 0.6 and 39.1 Å for the 25 and 44 wt % systems, respectively, the experimental value being 39.6 Å.35 This guarantees the quality of the present analysis for the electrostatics of the DPPC bilayer in the liquid crystal phase. The calculated molecular membrane area was 53.9 ( 0.6 and 53.6 Å2 for 25 and 44 wt % systems, respectively. These values are slightly smaller than the experimental one of 62.9 Å2.35 B. System Polarization. To characterize the bilayer in water as a dielectric substance, it is interesting to analyze the fluctuation of polarization of the system. The polarization P is defined by the dipole moment M of the system per unit volume:

Pt

(3.1)

where µi is the dipole moment of molecule i. The z-axis of the polarization vector P was chosen such that it is parallel to the bilayer normal, with the x- and y-axes then lying on the bilayer plane. Polarization rather than dipole moment is favorable to the present analysis, because the system volume fluctuates due to the flexible MD cell. To examine the polarization of lipid and water separately, the lipid bilayer region was defined by the space between phosphorus atoms in the upper and lower layers. Then, the volume of the lipid domain is

Vlipid t hxxhyydP-P

(3.2)

where h denotes the cell matrix, and dP-P is the instantaneous distance between the phosphorus atoms on each side averaged over the bilayer. Note that we adopt the triangular cell matrix in the MD run.21 Remaining volume in the cell is, thus, the water domain. This might be a rough definition for the volume. However, even if we take the membrane thickness to be dP-P ( 4 Å instead of the present dP-P, the change in the calculated susceptibilities χ was as small as (8%. According to the linear response theory,36 the susceptibility or specific dielectric constant of the system may be calculated by

III. Results and Discussion A. Electron Density Profile. Figure 1 shows the calculated electron density profiles along the bilayer normal (z-axis) for the systems at 25 and 44 wt % water. A good agreement of the plot among the three calculations for the membrane at 25 wt % as well as the almost completely symmetrical distribution

∑

iµi M ) V V

〈|∆M|2〉 ) β〈|∆P|2〉〈V〉 〈V〉

χ0 ) β

(3.3)

where 0 is the dielectric constant of vacuum, V is the system volume, and β )1/kBT. The averages in the equation should be taken over a very long time covering all membrane dynamics.

Molecular Dynamics Study of a Lipid Bilayer

J. Phys. Chem. B, Vol. 102, No. 34, 1998 6649

Figure 2. Time evolution of x- and z-components of the polarization of the total system, lipid, and water regions for the system at (a∼f) 25 wt % water and (g∼l) 44 wt % water.

However, since our calculation time is too short to include all of them, the susceptibility evaluated here is an apparent one for about 10-9 s. Furthermore, the susceptibility calculated for TIP4P water cannot reproduce the experimental one satisfactorily, the former and the latter being 53 ( 237 and 80, respectively, at 293 K. Thus, we restrict our present discussion to a qualitative one. Time evolutions of the polarization Px and Pz for the total system at two kinds of water concentrations are plotted in Figure 2 together with the contributions from the lipid, Plx and Plz, and water, Pwx and Pwz. Plots of Py are not presented in the figure because they are essentially the same as those of Px. The most striking finding in the figure is that the fluctuation of Pz for the total system is much smaller than that of Px. This implies that there is a large anisotropy in the dielectric constant of the membrane + water system. The large fluctuation of Px for the total system seems to have a period of several hundred picoseconds. This must come from the lipid polarization. In fact, fluctuation of the lipid polarization Plx is of almost the

same shape as that of the total system Px. Although the x-component of water polarization Pwx fluctuates very frequently with a considerably large amplitude, no correlation is found between lipid and water. In contrast to this, it is clearly shown that the z-component of water polarization Pwz fluctuates such that it cancels out the lipid polarization Plz. Then, the fluctuation of Pz for the total system is very small as shown in Figure 2, panels d and j. These result in anisotropic behavior of the susceptibility as shown in Table 2. The value of χzz is more than 1 order of magnitude smaller than that of χxx and χyy. The difference is evident despite large variations of the calculated susceptibilities among the MD runs caused by the present short time calculations compared with the slow dynamics of the membrane. Second, it is notable that the fluctuation of Plz and Pwz, itself, is much smaller than that of Plx and Pwx, respectively, presenting the smaller regional susceptibilities in the z direction than in the x and y directions as listed in Table 2. The small χzz for the lipid may arise from the very small z-component of the

6650 J. Phys. Chem. B, Vol. 102, No. 34, 1998

Shinoda et al.

TABLE 2: Diagonal Components of Susceptibility χ for the DPPC Bilayer at 44 and 25 wt % Water Concentrationsa 44 wt % 25 wt %

χxx χyy χzz χxx χyy χzz

total

lipid

water

20 19 1.7 29 ( 21 46 ( 32 1.2 ( 0.1

8.4 9.2 2.4 35 ( 25 45 ( 37 3.9 ( 0.8

34 29 3.9 27 ( 1 27 ( 2 8.5 ( 1.3

a The susceptibility was calculated separately from the fluctuation of polarization of the total system, lipid, and water regions using eq 3.

Figure 3. Electrostatic potential profile across the membrane. Center of the water lamella is chosen to be at z ) 0. (a) 44 wt % water; (b) 25 wt % water. Solid line, total system; dashed line, lipid; dotted line, water. Error bars in the figure for the 25 wt % water system represent the same as in Figure 1.

dipole moment in the headgroup; a vector made by the P and N atoms inclines at about 70°, on average, with respect to the bilayer normal, and further, the dipole moments in the opposite directions on both layers cancel out each other. The small fluctuation of Pwz, on the other hand, may be caused by the restricted motion of water by a large electric potential gradient along the bilayer normal. This will be discussed in detail in the next subsection. C. Electrostatic Potential Profile across the Membrane. The electrostatic potential profile along the bilayer normal was calculated based upon Poisson’s equation:

φ(z) - φ(0) ) -

1 0

∫0z ∫0z′ F(z′′) dz′′ dz′

(3.4)

where φ(z) and F(z) are the electrostatic potential and the charge density at z, respectively. The origin of the z-axis was chosen at the center of water lamella. In Figure 3, the results are plotted for the system at two kinds of water concentrations. The potential curves due to the charge distribution of lipid and water were evaluated separately and plotted in the figure, too. A slight discrepancy of the potential difference between two interfaces, i.e., the asymmetry of the potential profile of water and lipid, is a measure of the statistical error. However, the error is satisfactorily small for the both systems. The two plots show the same tendency that lipid dipoles contribute to make a negative potential inside the bilayer. This is canceled by water polarization excessively. Consequently, a positive potential as large as about 1 V with respect to the water region is found in the bilayer. The system at the other water concentration shows almost the same value of the potential difference. This implies that there is no significant difference in the structure at the interface due to the water concentration. Thus, the discussions below will be based upon the analysis for the system at 25 wt % water concentration. The positive potential in the bilayer interior indicates that permeation of anions across the membrane is easier than that of cations. This is in good correspondence

Figure 4. Density profiles of polar atoms across the membrane. Center of the water lamella is chosen to be at z ) 0. (a) 44 wt % water; (b) 25 wt % water. Solid line, phosphorus; dashed line, nitrogen; dotted line, carbonyl carbon in the sn-1 chain; dotted-dashed line, carbonyl carbon in the sn-2 chain; thin solid line, water. Error bars for the 25 wt % water system represent the same as in Figure 1.

with the experimental observation of the permeation of hydrophobic ions through lipid bilayers.1 Furthermore, this large potential gradient must restrict the reorientation of water molecules with respect to the bilayer normal at the interface. This may be one of the origins for the small fluctuation of the polarization of water in z direction as stated in the previous subsection. A stepwise variation of total electrostatic potential along the bilayer normal was found; two steps are likely to be found in the figure. It is clear from Figure 4, atom distributions along the bilayer normal, that the first potential gradient is located at the position of the ester group in the bilayer, ∆φ ∼ 300 mV, and the other is found over the range of about 10 Å from phosphate to water, ∆φ ∼ 700 mV. The averaged position of the relevant atoms in the lipid molecule as well as its standard deviation summarized in Table 3 may be helpful to identify this. The large potential gradient around the phosphate group is characteristic of phosphatidylcholine (PC), since it was not found in phosphatidylethanolamine (PE).4 The potential gradient near the ester group is produced in a similar way to the case of DLPE.4 Distribution of the angle between the CO vector in the ester group and the bilayer normal is plotted in Figure 5. Both types of CO vectors, one from the carbonyl C to the carbonyl O and another to the ether O, have a preferential orientation pointing toward the outside of the bilayer. These are the origin of the potential gradient found in this region. Figure 6 presents the radial distribution functions of water around the carbonyl and ether oxygens. Because of the heterogeneity of the membrane system, these functions do not reach 1 in the range of distances shown in the figures. The water molecules hydrated to the carbonyl oxygen were found to have a preferential orientation, i.e., the water hydrogen is located closer to the carbonyl oxygen than the water oxygen, resulting in two sharp peaks in the distribution functions. Although these water molecules might reduce the dipole moment produced by the carbonyl dipoles, the number of the hydrated water molecules is very small. The calculated number of water molecules in the first hydration shell, i.e., within the distance of 3.8 Å, was about 0.6. Further, the water molecules that are hydrogen-bonded to the ether oxygen were rarely found as shown in Figure 6b; little preferential orientation is observed. Thus, the dipole moment formed by the hydrated water molecules is not enough to cancel the ester group dipoles. In fact, we can find in Figure 3 that the total potential gradients near the ester group are formed solely by the lipid dipole potential; no contribution from water is observed. On the other hand, with respect to the second slope around the phosphate group, it is clearly found in Figure 3 that the electric potential from water exceeds that from the lipid in this

Molecular Dynamics Study of a Lipid Bilayer

J. Phys. Chem. B, Vol. 102, No. 34, 1998 6651

TABLE 3: Average Position and Its Standard Deviation of Polar Atoms Along the Bilayer Normala 44 wt % 25 wt % a

upper layer lower layer upper layer lower layer

zP (Å)

zN (Å)

zC1γb (Å)

ZC1βc (Å)

20.5 ( 1.9 -20.5 ( 2.2 9.0 ( 2.0 -9.0 ( 2.0

19.4 ( 2.9 -19.1 ( 3.1 7.8 ( 3.0 -7.7 ( 2.8

25.4 ( 2.1 -25.9 ( 2.3 13.9 ( 2.1 -14.2 ( 2.0

24.5 ( 1.9 -24.4 ( 1.9 12.8 ( 2.1 -12.7 ( 1.8

The center of water lamella was chosen to be z ) 0. b The ester carbons in the sn-1 chain. c The ester carbons in the sn-2 chain.

Figure 5. Distribution function of the angle θ between CO vectors of the ester group and the bilayer normal for the system at 25 wt % water. Solid line, vector from carbonyl carbon to carbonyl oxygen (dO); dashed line, vector from carbonyl carbon to ether oxygen (-O-). Error bars represent the same as in Figure 1.

Figure 7. Radial distribution function g of water molecule around the headgroup: (a) phosphorus and (b) nitrogen. Solid line, water oxygen; dashed line, water hydrogen; dotted line, coordination number n of water oxygen. Error bars represent the same as in Figure 1.

Figure 6. Radial distribution function g of water molecule around the ester oxygens in the system at 25 wt % water. (a) Water molecule around the carbonyl oxygen (dO); (b) water molecule around the ether oxygen (-O-). Solid line, water oxygen; dashed line, water hydrogen; dotted line, coordination number n of water oxygen. Error bars represent the same as in Figure 1.

region. In other words, the potential gradient by water is steeper than that by the lipid. To investigate water structure in this region, radial distribution functions of oxygen and hydrogen of the water molecule around the phosphorus and nitrogen atoms were calculated and plotted in Figure 7. Very sharp first peaks in the distributions of P-water (water molecules hydrating the phosphorus atom) show that the hydrogen atom is located closer to the phosphorus by about 1 Å than the oxygen atom. This represents a strong linear hydrogen bond between the phosphate group and water. In contrast to this, g(r)s of N-water (water molecules hydrating the nitrogen atom) show broad distribution. This implies that the water molecules are weakly bonded to N compared with P. Difference of position of the first peaks between N-hydrogen and N-oxygen is very small, indicating little preferential orientation of water. Thus, the orientation of

Figure 8. Polarization profile of water along the bilayer normal. Center of the water lamella was chosen to be at z ) 0. Solid line, 25 wt % water system; dashed line, 44 wt % water system. Error bars for 25 wt % water system represent the same as in Figure 1.

the water molecule near the interface is dominated mostly by the phosphorus atom. Collective orientational order of water dipoles along the bilayer normal gives rises to a potential gradient. The calculated polarization profile defined by the integration of charge density along the z axis

P(z) )

∫0z F(z′) dz′

(3.5)

is shown in Figure 8. Water polarization is found over the range of about 10 Å at the interface. This is too rough to find the origin of potential gradient in this region, because the membrane/ water interface is intrinsically bumpy. Here, it is desirable to investigate the interface for each lipid molecule separately. Consider a cylindrical coordinate (r, ψ, z) whose z axis is perpendicular to the bilayer plane and the origin is chosen on the position of the phosphorus atom of the lipid molecule.

6652 J. Phys. Chem. B, Vol. 102, No. 34, 1998

Shinoda et al.

Figure 10. Residence time of water molecules as a function of z along the bilayer normal, where the position of P was chosen to be at z ) 0. (O) 25 wt % water system; (b) 44 wt % water system.

correlation function: Figure 9. Orientational distribution and density profile of water in the slab at z ) z in the cylinder. For details, see text. (a) Water dipole orientational distribution, where θ is the angle between water dipole and the bilayer normal. (b) Density profile of the water oxygen. Solid line, 25 wt % water system; dashed line, 44 wt % water system. Error bars for 25 wt % water system represent the same as in Figure 1.

C(z,t) )

Figure 9a shows distribution of cos θ, where θ is the angle between the water dipole and the bilayer normal, for the water molecules at πr2 < 53.9 Å2, averaged membrane area per lipid molecule, as a function of z. It is clearly found that the water near the phosphorus tends to direct its dipole toward the phosphorus. The water dipole orientation depends, thus, solely on its relative position with respect to the phosphate group. Probability distribution of the water molecule in the cylinder as a function of z is also plotted in Figure 9b. The figure clearly shows an anisotropic hydration around the phosphorus; the number of water molecules is very small at negative z compared with positive z. A significant increase of water molecules is found for z from -5 to 5 Å. This demonstrates that contribution of the water dipole to the electrostatic potential in the z direction is predominantly made by a large amount of water molecules at z ) 0 ∼ 5 Å. In this case, cancellation of the dipole moment usually found in the symmetric hydration is not observed. This must be the main origin of the second potential slope shown in Figure 3. A peak of the density profile of the water oxygen at about 3.5 Å corresponds to the position of the first hydration shell; see Figure 7. The density profile shows an almost constant value only in the range of z > 10 Å, because the choline groups occupy a space at z < 10 Å. However, the choline group does not have a significant effect on the water orientation along the z axis. In the case of PE, the situation may be different to some extent, because the highly charged hydrogen on the ethanolamine groups will make a stronger hydrogen bond with water. In the MD calculation of the DLPE bilayer, the membrane dipole potential was produced mainly by the ester groups (about 450 mV) and the contribution by the PE headgroup region is small (about 150 mV).4 D. Dynamics of Water. In the previous section, it is clearly observed that the collective motion of water with respect to total polarization is strongly affected by the lipid dipole potential, which results in the anisotropy of the susceptibility. In this subsection, we focus our attention on the dynamics of the single water molecule near the interface. To do this, here, trajectories of H2O were analyzed. 1. Position Dependence. First, an averaged residence time of water molecules at z ) z was investigated by the following

where hi(z,t) was defined as

)

hi(z,t) )

{

〈hi(z,t)hi(z,0)〉 〈hi(z,0)hi(z,0)〉

1 〈h (z,t)hi(z,0)〉 Nz i

(3.6)

1 if the water molecule i is found at z ) z ( ∆z from t ) 0-t (3.7) 0 otherwise

and Nz is the average number of water molecules in the region. Here, the value of ∆z was chosen to be 2.5 Å. Relaxation times, τ, obtained by assuming an exponential form of the correlation function at large t and integrating over t, are plotted in Figure 10 as a function of z. The origin of the abscissa was chosen at the position of the phosphorus. A systematic change of the residence time along the bilayer normal was found. The residence time of water in the interfacial region, |z| < 5 Å, is about three times as long as that in the bulk water, z > 10 Å. A very large value of τ for water at z < 5 Å may be explained by the hydrogen bonding between water and the phosphate groups. In any case, even in the interfacial region, the residence time of the water molecule is about 2 orders of magnitude shorter than the time scale of the headgroup motion, the time constant of the reorientational motion of the headgroup being about 1 ns.21 Reorientational motions of water molecules in the water slab at z are characterized by the following autocorrelation function of the water dipole moment:

P1(z,t) ) 〈cos θ(t)〉z ) 〈µˆ (t)µˆ (0)〉z

(3.8)

where µˆ (t) is the unit dipole vector at time t. Sampling was done for the molecules found at z ) z ( ∆z at t ) 0. Relaxation time was estimated by the same procedure as the case of the residence time. The results are plotted in Figure 11a. It is clearly found that the water molecule closer to the membrane has a slower reorientational motion. The water at z > 10 Å, the bulk water, shows the fastest reorientation, 1.2 ps, which is the same as that estimated from the MD run of pure water. For the system at lower water concentration, reorientation of all water molecules is restricted to some extent; no bulk water is found. Furthermore, to investigate an anisotropy of the rotational motion of water molecule, the rotational correlation functions for the unit dipole vectors projected on the x-y, y-z, and z-x planes were evaluated. The relaxation time is also plotted in

Molecular Dynamics Study of a Lipid Bilayer

J. Phys. Chem. B, Vol. 102, No. 34, 1998 6653

Figure 11. Rotational relaxation time of water as a function of z along the bilayer normal, where the position of P was chosen to be at z ) 0. (a) Three-dimensional reorientation; (b) two-dimensional reorientation. (4, 2) τx-y; (0, 9) τy-z; (3, 1) τz-x. Open symbols, 25 wt % water; closed symbols, 44 wt % water system.

Figure 11b. Water molecules at the boundary region, z < 5 Å, have shorter τx-y than τy-z and τz-x. The anisotropic rotation in these layers may be caused by the large electric potential gradient along the z axis; it is hard for the molecule to rotate freely with the changing z component of the dipole moment. The large value of the relaxation time in these regions is due to the rigid hydration to the phosphate groups. In contrast, the water molecules at z > 5 Å show isotropic rotation that is almost the same as that of the bulk water. Mean square displacements (MSD) of the water molecules at z were also calculated to estimate the diffusion coefficient D(z):

D(z) )

1 1 lim 〈|r(t) - r(0)|2〉z 2d tf∞ t

(3.9)

where d denotes the dimension, and MSD sampling was carried out for the water found at z ) z ( ∆z at t ) 0. Since the water molecules will migrate out of the home slab at large t, the diffusion coefficient was estimated using short time MSD from 2 to 7 ps. Water molecules at the center of water lamella for the system of 44 wt % water were found to have a large D value, 8.8 × 10-5 cm2 s-1, which is very close to that of the bulk water; the D value of pure TIP4P water was calculated to be about 8.5 × 10-5 cm2 s-1. The calculated D for the usual three-dimensional diffusion as well as one-dimensional diffusion are plotted in Figure 12. The abscissa is the same as that in Figure 11. In the systems at both concentrations, slower diffusive motion of the water molecules near the interface was found than for the bulk water. Furthermore, the lower translational mobility of water along the bilayer normal, i.e., smaller Dz than Dx and Dy, was found in almost all layers, though the difference is small at the interfacial region. 2. Water Molecules Coordinated to the Headgroups. As shown in the previous section, the water molecules bound to the headgroups, in particular to the phosphate group, play an important role in producing electrostatic potential gradient at the interface. In this subsection, we investigate the dynamics of the water molecules in the first hydration shell of the phosphorus and nitrogen atoms. The residence correlation function of the water hydrated to the phosphate and choline groups were calculated on the basis of a function similar to eq 3.6, separately. The residence time

Figure 12. (a) Three-dimensional and (b) one-dimensional diffusion coefficients of water as a function of z along the bilayer normal. (4, 2) Dx; (0, 9) Dy; (3, 1) Dz. Open symbols, 25 wt % water; closed symbols, 44 wt % water system.

TABLE 4: Rotational Relaxation Time τ and Diffusion Coefficient D of Water Hydrated to the Head Group P-water τ (ps) τx-y (ps) τy-z (ps) τz-x (ps) D (× 10-5 cm2 s-1) Dx (× 10-5 cm2 s-1) Dy (× 10-5 cm2 s-1) Dz (× 10-5 cm2 s-1)

N-water

25 wt %

44 wt %

25 wt %

44 wt %

16.2 11.6 16.9 16.5 1.6 1.7 1.7 1.6

17.0 12.1 17.6 17.7 1.7 1.8 1.9 1.6

8.6 6.1 9.0 9.0 2.8 3.0 3.1 2.3

9.9 6.4 9.2 9.4 3.4 3.6 3.7 3.0

of P-water and N-water were estimated to be about 8.3 and 5.9 ps, respectively, for the system at 44 wt % and 7.9 and 5.8 ps, respectively, for the system at 25 wt %. The difference of the residence time between the two systems at different water concentrations is within the fitting error that is as large as about 1 ps. However, a longer residence time of P-water than N-water is clearly found. This may be caused by the strong hydrogen bonding of the water hydrogen with the branched oxygen (not with the backbone oxygen) of the phosphate group. In any case, the residence time of the hydrating water is 2 orders of magnitude shorter than the characteristic time scale of the headgroup reorientational motion (∼1 ns).21 Thus, the hydrating water molecules are exchanged by the bulk molecules several tens of times while the headgroups change their conformation once. Reorientational and translational motions of the water bound to the headgroups were also investigated in the same manner as in the previous section. The calculated rotational relaxation time τ as well as diffusion coefficient D is listed for P- and N-water molecules in Table 4. A lower mobility of P-water than of N-water was clearly found. A slower reorientational motion of water dipole in the x-y plane than in the y-z and z-x planes was also observed for P-water as well as for N-water. This rotational anisotropy must be caused by the dipole ordering along the z axis near the headgroups (cf. Figure 9). IV. Conclusions Molecular dynamics calculations of the DPPC bilayer in the liquid crystal phase have been performed in the Nose´Parrinello-Rahman NPT ensemble at two water concentrations. No significant difference in the bilayer structure between two

6654 J. Phys. Chem. B, Vol. 102, No. 34, 1998

Shinoda et al. References and Notes

Figure 13. Schematic picture of the direction of dipoles at the interface. The signs and arrows in the figure represent the signs of partial charges on the atoms and the dipole moments, respectively.

systems was found, showing a good correspondence to a recent experimental result.35 Polarization of the system showed an anisotropic fluctuation; the component along the bilayer normal was 1 order of magnitude smaller than that in the bilayer plane. This results in the anisotropy of the susceptibility. It is found, too, that collective polarization of water is strongly correlated with the membrane structure such that it cancels out the polarization made by the lipid giving a very small total fluctuation along the bilayer normal. A positive electrostatic potential inside the membrane was also found. The potential was found to vary in two steps across the interface between the membrane and water; the first gradient comes from the preferential orientation of the ester group, the second one from water dipole contribution around the phosphate groups. In Figure 13, a schematic picture of the direction of dipoles in the bilayer system is presented. The ester dipole moment contributes to the total potential due to the lack of the water in the region. The water asymmetrically coordinated to the phosphate group presents a large potential gradient, too, canceling the dipole of the headgroups excessively. Water molecules in the potential gradient showed anisotropic rotation; a slower reorientation of the dipole vector projected on the x-z and y-z plane was found. Water dynamics near the headgroups was strongly influenced by the hydrogen bonding with the phosphate group, too. Acknowledgment. We thank the computer centers of the Institute for Molecular Science, Tokyo Institute of Technology, and Japan Atomic Energy Research Institute for the use of supercomputers. The work was supported in part by a Grantin-Aid for Scientific Research (09440197) from the Ministry of Education, Science and Culture, Japan. W.S. has been granted a fellowship of the Japan Society for the Promotion of Science for Japanese Junior Scientists.

(1) Gennis, R. B. Biomembranes; Molecular Structure and Function; Springer-Verlag: New York, 1989. (2) Cevc, G.; Marsh, D. Phospholipid Bilayers; Wiley-Interscience: New York, 1987. (3) Marrink, S. J.; Berkowitz, M.; Berendsen, H. J. C. Langmuir 1993, 9, 3122. (4) Zhou, F.; Schulten, K. J. Phys. Chem. 1995, 99, 2194. (5) Chiu, S. W.; Clark, M.; Balaji, V.; Subramaniam, S.; Scott, H. L.; Jakobsson, E. Biophys. J. 1995, 69, 1230. (6) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871. (7) Berger, O.; Edholm, O.; Ja¨hnig, F. Biophys. J. 1997, 72, 2002. (8) Feller, S. E.; Pastor, R. W.; Rojnuckarin, A.; Bogusz, S.; Brooks, B. R. J. Phys. Chem. 1996, 100, 17011. (9) Cascales, J. J. L.; Berendsen, H. J. C.; De la Torre, J. G. J. Phys. Chem. 1996, 100, 8621. (10) Damodaran, K. V.; Merz, K. M., Jr.; Gaber, B. P. Biochemistry 1992, 31, 7656. (11) Venable, R. M.; Zhang, Y.; Hardy, B. J.; Pastor, R. W. Science 1993, 262, 223. (12) Stouch, T. R. Mol. Simul. 1993, 10, 335. (13) Heller, H.; Schaefer, M.; Schulten, K. J. Phys. Chem. 1993, 97, 8343. (14) Essex, J. W.; Hann, M. M.; Richards, W. G. Philos. Trans. R. Soc. London B 1994, 344, 239. (15) Robinsona, A. J.; Richards, W. G.; Thomas, P. J.; Hann, M. M. Biophys. J. 1994, 67, 2345. (16) Egberts, E.; Marrink, S. J.; Berendsen, H. J. C. Eur. Biophys. J. 1994, 22, 423. (17) Essmann, U.; Perera, L.; Berkowitz, M. L. Langmuir 1995, 11, 4519. (18) Perera, L.; Essmann, U.; Berkowitz, M. L. Langmuir 1996, 12, 2625. (19) Shinoda, W.; Fukada, T.; Okazaki, S.; Okada, I. Chem. Phys. Lett. 1995, 232, 308. (20) Tu, K.; Tobias, D. J.; Klein, M. L. Biophys. J. 1995, 69, 2558. (21) Shinoda, W.; Namiki, N.; Okazaki, S. J. Chem. Phys. 1997, 106, 5731. (22) Merz, K. M., Jr.; Roux, B. Biological membranes. A molecular perspectiVe from computation and experiment; Birkha¨user: Boston, MA, 1996. (23) Feller, S. E.; Zhang, Y.; Pastor, R. W. J. Chem. Phys. 1995, 103, 10267. (24) Feller, S. E.; Pastor, R. W. Biophys. J. 1996, 71, 1350. (25) Feller, S. E.; Pastor, R. W. In Proceedings of the Pacific Symposium in Biocomputing; World Scientific: Singapore, 1997; pp 142-150. (26) Nose´, S. J. Chem. Phys. 1984, 81, 511. (27) Parrinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196. (28) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182. (29) Nose´, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055. (30) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (31) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (32) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G., Jr.; Profeta, S.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (33) Gear, C. W. ANL Report 1966, No. ANL-7126. (34) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1989. (35) Nagle, J. F.; Zhang, R.; Tristram-Nagle, S.; Sun, W.; Petrache, H. I.; Suter, R. M. Biophys. J. 1996, 70, 1419. (36) Kubo, R. J. Phys. Soc. Jpn. 1957, 12, 570. (37) Neumann, M. J. Chem. Phys. 1986, 85, 1567.