1042
Langmuir 1991, 7, 1042-1044
Computer Simulation of a Water/Membrane Interface Max L. Berkowitz’ and K. Raghavan Department of Chemistry, University of North Carolina, Chapel Hill, North Carolina 27599 Received January 28, 1991. In Final Form: March 13, 1991 In an effort to understand the origin of the hydration force, we performed a molecular dynamics computer simulation on a system composed of water molecules embedded between surfaces of dilauroylphosphatidylethanolamine (DLPE)molecules. The simulations show that the orientational polarization of water displays an oscillatory decay and its behavior is independent of thermal motion of the head groups of the membrane molecules. Since the Landau-type theory proposed by Marcelja and Radic predicts a monotonic decay for the order parameter, we conclude that the orientational polarization is not a good order parameter. Instead, we propose a new quantity to measure the order parameter, “hydrogen bond deficiency”, and calculate its behavior from the simulations. It agrees well with the theoretical predictions.
Introduction In 1976 LeNeveu, Rand, and Parsegian provided the first accurate data on strong repulsive forces observed to act between electrically neutral lecithin bilayers.’ Later it was shown that the same kind of force dominates the interaction between surfaces of most neutral or charged phospholipid membranes a t distances below -2-3 nm.2 Since it was argued that the force is derived from the work required to remove water of hydration, the force was termed the hydration force. The hydration force decays exponentially with the distance between membrane surfaces and the value of the decay constant is usually in the range of -0.1-0.25 nm. This decay constant is independent of the electrolyte concentration but depends on the properties of the surfacese2 Immediately after the publication of the first data on the hydration force, Marcelja and Radic proposed an elegant theory to explain the nature of the observed strong force.8 According to their theory the force is due to the modification of water structure near the membrane/water interface, where the water molecules are different from the water molecules in the bulk; they are more “ordered”. To describe this “order” an order parameter ~ ( x )is introduced and a Landau expansion of the free energy density is performed; i.e. the free energy density of the system is written as g = go
+ a?+ + ... + C [ d T / d X ] 2
(1)
Assuming that the interfaces are positioned a t x = h / 2 and x = - h / 2 and that ~ ( h / 2=) -q(-h/2) = 70, the minimization of free energy given by eq 1 results in the differential equation d2a(X)/dx2- ( a / ~ ) ~ = ( x0)
(2) The solution of eq 2 subject to the boundary conditions is ~ ( x= ) T~ sinh [(x/X)]/sinh [h/2X]
(3)
where X = ( ~ / a ) ’ / It ~ .is then simple to show that if the free energy is given by eq 1 and the order parameter is given by eq 3, the repulsion force follows an exponential decay.l In the 1976 paper Marcelja and Radic did not specify what quantity plays the role of the order parameter, but (1) LeNeveu, D. M.; Rand, R. P.; Parsegian,
V. A. Nature 1976,259,
601-603. (2) Rand, R. P.; Parsegian, V. A. Biochim. Biophys. Acta 1989,988, 351-376. (3) Marcelja, S.; Radic, N. Chem. Phys. Lett. 1976, 42, 129-130.
\,
I
.
b
X
Figure 1. A snap shot of a part of the DLPE system used in the simulation. The water molecules are omitted for clarity. The vertical lines denote the location of the molecular surfaces, as described in the text. in later publications the orientational polarization of water molecules was considered to be the order parameter in the problem.43 To learn more about the behavior of water next to membranes, Kjellander and Marcelja performed a molecular dynamics computer simulation of water embedded between two lecithin surface^.^*^ Since these studies were performed five years ago, the researchers were limited by the computer power then available to them. Thus the dynamics of water molecules was performed in the presence of fixed lecithin head groups and the runs were rather short (8 P S ) . ~From the simulations of Kjellander and Marce1ja6v7it was observed that the orientational polarization propagated in a damped oscillatory manner. The findings from the Kjellander-Marcelja simulation questioned the wisdom of the choice of polarization as an order parameter. Moreover, since the head groups of the membrane molecules in the Kjellander-Marcelja simulations were fixed, the following question was asked: ”Could the monotonically repulsive hydration force between bilayers be a smeared-out oscillatory force, due to the thermal motion and roughness of the mobile head groups?”.8
Methods To answer the previous question and other questions about water/membrane interface, we have recently performed a simulation of a bilayer composed of 32 phospholipid molecules and (4) Gruen, D. W. R.; Marcelja, S. J. Chem. Soc., Faraday Trans. 2 1983, 79, 211-223. ( 5 ) Gruen, D. W. R.; Marcelja, S. J. Chem. Soc., Faraday Trans. 2 1983, 79, 225-242. (6) Kjellander, R.; Marcelja, S. Chem. Scr. 1985, 25, 73-80. (7) Kjellander, R.; Marcelja, S. Chem. Phys. Lett. 1985,120,393-396. (8) Israelachvili, J. N.; Wennerstrom, H. Langmuir 1990,6,873-876.
0743-746319112407-1042$Q2.50/0 0 1991 American Chemical Society
Letters
Langmuir, Vol. 7, No. 6, 1991 1043
..f-----";i +
Density
4.00
2.00
0.00
1.6
2.4
3.2
4.0
4.8
1.6
3.2
2.4
4.8
4.0
X (nml X (nml Figure 2. (a) Oxygen (solid line) and hydrogen (dashed line) densities of water (in units of g/cm3) between DLPE surfaces shown along the x coordinate. (b) Number of nearest neighbors per molecule (solid line) and the number of hydrogen bonds per molecule (dashed line). The horizontal lines show the bulk values. The locations of the DLPE surfaces are indicated by the arrows. (c and d) The densities (in arbitrary units) of the atoms comprising the backbone of the head groups in the DLPE molecules. The results are from the simulation involving flexible head groups. 437 water molecules at T = 285 K. We have chosen the phos(DLPE) as pholipid 1,2-dilauroyl-~~-phosphatidylethanolamine the object of our study since the structure of the molecule has already been determined by using X-ray crystallography and the X-ray coordinates are available in the literat~re.~In our simulations, we have used 32 molecules of DLPE to construct the bilayer. We arranged 16 molecules (four layers, each containing four DLPE molecules) in such a way that the molecular surface formed by the head groups was parallel to the yz plane of the coordinatesystem and the hydrocarbon chains were pointing toward + x direction. The other 16moleculeswere placed such that their molecular surface was facing the previous one, but shifted along the 1: axis. The hydrocarbon chains of these molecules were pointing toward -x axis. Then the two surfaces were slightly shifted with respect to one another along they axis. With this arrangement, the x axis is perpendicular to the molecular surfaces. The distance between the molecular surfaces was about 1.6 nm. This is roughly the distance between the PO, groups on opposite sides. Figure 1 displays a cross section of a layer of DLPE molecules in our simulation. Also shown in the figure are the positions of the molecular surfaces. The choice of these positions will be explained below. The DLPE molecules were then immersed into a large box of water molecules and the water molecules beyond the glycerol carbon C1 were removed. Also, to avoid any bad contact,water moleculeswithin the distance of 0.13 nm from any atom of the DLPE molecule were removed from the system. The number of water molecules was adjusted to give the density of bulk water in the middle of the simulation box. This yielded 437 water molecules in the system. The dimensions of the simulation box for repeating the bilayer in three dimensions were chosen to be 6.24nm in the x direction, and 3.20 and 2.05 nm in y and z directions, respectively. The dimensions for y and z directions were chosen such that the area per lipid molecule is 0.41 nm2, which corresponds to the value of the area per lipid in the gel phase.10 The phospholipid (9) Elder, M.; Hitchcock, P.; Mason, R.; Shipley, G . G. Proc. R. SOC. London, A 1977,354, 157-170. (10)McIntosh, T. J.; Simon, S . A. Biochemistry 1986,25,494&4952.
ka) '
2.4
Polarization
2.0
3.2
3.6
' I
nt
4.0
X (nml
Figure 3. Polarization of water in the direction normal to the DLPE surface from the simulations using (a)flexible head groups and (b) rigid head groups. molecules were chosen to be in a gel phase because the hydrocarbon chains of the molecules were fixed during the simulations. Periodic boundary conditions were used in all three directions. The force field used in the simulations was taken from AMBER software package," properly modified to include membrane-type molecules. The detailed description of the (11) Singh,U.C.;Weiner,P.K.;Caldwell,J.K.;Kollman,P.K.AMBER,
University of California: San Francisco, 1987; Version 3.1.
1044 Langmuir, Vol. 7, No. 6,1991
Letters 2 -0 .
0.2
> 0
c 0
; 0.
o.1
4
.rl
0
4 ID
2
.rl
y.
0 .o
m
0
I ID
0.
U
rl
C
0
0 m -0.
n -0 .1
I
I
-0 . 2
2.4
2.8
3.2
3.6
4.0
X (nm)
-0. X (nml
Figure 4. Polarization and the "hydrogen bond deficiency" for water in the region between the DLPE surfaces. The circles represent the data obtained from the simulation. The dashed line shows the fit to eq 3 for polarization and for "hydrogen bond deficiency". simulation will be given elsewhere.12 Two molecular dynamics runs, each of 100 pa length were performed. In the first run the head group atoms of the membrane molecules were flexible,while in the second run the head groups were fixed in some random configuration. As we have mentioned, the hydrocarbon chains of the molecules were fixed in both runs.
Results and Discussion The results from the simulations are displayed in Figures 2 to 4. Figure 2a shows the density of water along the x axis (axis perpendicular to the bilayer surface), and Figure 2b displays the profiles for the number of hydrogen bonds per water molecule and the number of nearest neighbors along the x axis. Two water molecules were considered to be nearest neighbors if the distance between their oxygen atoms was less than 0.33 nm. In addition, for a given nearest neighbor pair, consider the angles formed by OH and 00 vectors, originating from the same oxygen. If any of these angles was less than 35O, the pair was considered to be hydrogen bonded. Parts c and d of Figure 2 display the density and the location of the atoms in the head group of the phospholipid molecules. One of the most difficult problems to be resolved when dealing with rough surfaces is to decide how to measure the intersurface distance. In the present simulation, guided by Figure 2b we positioned the left boundary of the surface a t x = 2.32 nm, while the right boundary was positioned a t x = 3.89 nm. At these points the curves for the number of nearest neighbors and for the number of hydrogen bonds change abruptly their character. Our choice is also consistent with the density profiles displayed in parts a, c, and d of Figure 2. As we can see from these figures, the head group atoms slightly protrude into the water as expected.1° Figure 3 shows the polarization profiles along the x direction. Figure 3a displays the profile obtained from the simulation with the flexible head groups, while Figure 3b displays the polarization profile from the simulation with fixed head groups. No significant differences are observed between these two profiles. Although the oscillations in the profiles are reduced compared to the ones observed in ref 6 they are still present in our simulation. From our results we conclude that the thermal motion of head groups does not smooth out the polarization oscillations. Since the polarization profile does not display the anticipated monotonic decay, we concluded that it cannot serve as an order parameter. In search for its replacement, (12)Raghavan, K.; Rami Reddy, M.;Berkowitz, M. In preparation.
we considered the following quantity q=v
No - N ( x ) No
(4)
In this equation N ( x ) is the number of hydrogen bonds per water molecule a t point x and NOis the number of hydrogen bonds per water molecule in the bulk. This choice allows us to apply the same basic equations of Marcelja-Radic theory (eqs 1-3) to a new order parameter defined by eq 4. Figure 4 displays the behavior of both choices for the order parameter obtained from the simulation vs the theoretical predictions. As the figures show,the "hydrogen bond deficiency" defined by eq 4 is probably a better choice. From the fit to the simulation data, we estimate that the decay constant X for the new order parameter is -0.15 nm, a value in the range of the experimental data.2 To summarize, our simulations show that the orientational polarization of water displays an oscillatory behavior that does not depend on thermal motion of the head groups of the membrane molecules. We also conclude that the Marcelja-Radic theory of the hydration force is still applicable but the order parameter in this theory might be the "hydrogen bond deficiency" as defined in eq 4. Therefore, the observed repulsive hydration force is due to the work required to break the hydrogen bonding network present in water. Finally, a word of caution about the results from these simulations. The numerical data obtained in the simulations should not be taken literally, but rather used heuristically, since the potential functions in the simulations are known only approximately, and the run times are still relatively short (although an order of magnitude larger than in the previous simulations). These words of caution should apply to any simulations performed on complex systems, such as biomolecules and their assemblies.
Acknowledgment. This work was supported from a grant from the Naval Research Office. The simulations were performed on a Cray-YMP a t the North Carolina Supercomputing Center and CONVEX supercomputer a t UNC. We enjoyed fruitful conversations with Drs. J. J. Hermans, L. Pratt, S. Simon, T. McIntosh,andS. Leikin. The insightful comments of the referees are gratefully acknowledged.
