ABSTRACT We present the results of a 10-ns molecular dynamics simulation of a dipalmitoylphosphatidylcholine/water system. The main emphasis of the present study is on the investigation of the stability over a long time and the dynamic properties of the water/membrane system. The motion of the lipid molecules is characterized by the center of mass movement and the displacement of individual atom groups. Because of the slow movement of the headgroup atoms, their contributions to the dipole potential vary slowly and with a large amplitude. Nevertheless, the water molecules compensate the strong fluctuations and maintain an almost constant total dipole potential. From the lateral displacement of the center of masses, we calculate the lateral diffusion coefficient to be D^sub lat^ = (3 +/- 0.6) X 10^sup -7^ cm^sup 2^/s, in agreement with neutron scattering results. The rotational motion is also investigated in our simulations. The calculated value for the rotational diffusion coefficient parallel to the molecular long axis, D^sub ll^ = (1.6 +/- 0.1) X 10^sup 8^ s^sup -1^, is in good agreement with the experiment.
INTRODUCTION
Over the last 15 years an increasing number of computer simulation studies have contributed to the understanding of the molecular structure of membranes (for reviews see Pastor, 1994; Mouritsen and Jorgensen, 1997; Tieleman et al., 1997; Tobias et al., 1997). Because the main focus of these studies is on structural properties, dynamical properties of the membranes have usually been neglected. This is unfortunate, because the fluid mosaic model of Singer and Nicholson strongly suggests that membrane fluidity plays a major role in the biological properties of the membrane (Singer and Nicholson, 1972; Jacobson et al., 1995). A notable exception is the work of Pastor and Feller (1996), in which the authors discuss the time scales of lipid dynamics and the relationship of these time scales to molecular dynamics. The range of time scales for different types of molecular motions in membranes is impressively large. On the shortest time scales are intermolecular vibrations, with a period on the order of 10^sup -14^ s, and on the other end of the scale are transbilayer flip-flop processes across the bilayer, with a typical time scale of minutes to hours (Blume, 1993). Clearly, with detailed atomistic simulation techniques, only the faster motions are currently accessible. Another limitation to the study of dynamical processes arises from the fact that many slower processes are collective excitations involving large spatial domains. Current molecular simulations, however, rarely exceed 100 Angstrom in any direction.
In the present study we attempt to investigate the dynamical properties of individual lipid molecules. To cover a broad spectrum of modes we extended the total simulation time to 10 ns. As we shall see, this allows the calculation of experimentally observed transport properties such as lateral and rotational diffusion coefficients. The large time scale also allows us to investigate the stability of the membrane system. It has recently been questioned whether lipid-water systems are mechanically stable in constant volume simulations (Shinoda et al., 1997). The claim was that only constant pressure simulations can lead to stable trajectories. However, as we shall see below, a reasonable choice of box parameters leads to stable trajectories in a constant volume ensemble.
Another issue that can be investigated by performing long time simulations is the issue of the validity of certain results obtained previously from relatively short time simulations on membranes. Because of the great internal flexibility of the lipid molecules, individual groups, for example, the headgroup or the lipid chains, move with large amplitudes. This is of particular concern for the headgroups because the headgroups produce an electric dipole potential of several volts. Therefore the question arises whether the large movements of the headgroup lead naturally to large fluctuations of the dipole potential, thereby rendering the results of shorter simulations of several hundred picoseconds meaningless.
METHODOLOGY
Because the computational protocol of the present calculation is similar to that of our previous study (Essmann et al., 1995a) we repeat here only a few characteristics of this protocol. Sixty-four dipahnitoylphosphatidylcholine (DPPC) molecules and 1315 water molecules are arranged in the primary simulation cell to form a bilayer/water system. By applying three-dimensional periodic boundary conditions, a multilamellar system is generated. This choice seems to be reasonable because many experiments have been performed on multilayer systems. The choice of the number of water molecules corresponds to 20.5 water molecules/lipid or 33.5% water fraction. Recent x-ray experiments on the structure of DPPC bilayers found that important physical characteristics such as bilayer thickness and area per headgroup remain conserved between a water content of 15.5 waters/lipid and the hydration limit of 28 waters/lipid (Nagle et al., 1996).
Recently the concern has been raised that membrane systems cannot be simulated with a fixed box geometry (Shinoda et al., 1997). The present result, however, indicates that for a reasonable choice of physical parameters, the stability of the system is preserved. The disadvantage of a simulation with a fixed box geometry is the need to specify the three box lengths. In general the three box angles have to be specified too; however, except for crystal simulations, they are usually set at 90 deg. Assuming that the system is isotropic within the plane of the membrane, one can set the two box lengths in the lateral dimension equal to each other. Furthermore, setting the box angles at 90 deg, there are still two parameters left, the area per headgroup and the bilayer repeat distance. In the case of pure DPPC membranes in water these two parameters can be obtained, with reasonable accuracy, from experimental data. However, in many interesting cases, such as the incorporation of peptides or proteins into membranes, the total membrane area is not known and constant pressure simulations are necessary (Feller et al., 1995).
FLUCTUATIONS IN THE ELECTROSTATIC MEMBRANE POTENTIAL
The membrane dipole potential is thought to play an important role in regulating not only the structure and function of membranes (Brockman, 1994), but even more complex processes such as the adsorption of proteins (Nordera et al., 1997). Experimentally this quantity is rather difficult to measure (Brockman, 1994). In general for DPPC bilayers, a value between -0.2 V (Gawrisch et al., 1992) and -0.4 V (Pickar and Benz, 1978) is obtained from the measurements. Membrane dipole potentials have also been calculated in molecular dynamics simulations. In our previous study the dipolar potential was calculated from a 300-ps simulation. Because of the zwitterionic nature of the headgroup, a dipole potential is associated with these headgroups. This dipole potential, in turn, leads to an ordering of adjacent water molecules. Because the conformations of the lipid headgroups vary strongly, one might expect that the dipole potential varies strongly over time as well. This could also explain the differences between reported values (Marrink et al., 1993; Essmann et al., 1995a).
To study the character of the fluctuations in the dipole potential we divided the total run into blocks of 300 ps length and calculated the contribution of the headgroups, the contribution of the water molecules to the dipole potential, and the total potential according to the formula
where psi(z) denotes the dipole potential and p(z) the charge density. The z axis is chosen to be parallel to the membrane normal.
In Fig. 2 the contributions from the lipid headgroups as well as from the water molecules are shown together with the total potential. The electrostatic potential of the headgroup fluctuates by +/- 0.24 V around an average value of 2.24 V, and similarly, the fluctuations of the electrostatic potential of the water are +/- 0.23 V around the average value of -3.20 V, whereas the fluctuations of the total electrostatic potential around the average value of -0.96 V are only +/- 0.03 V and therefore are one order of magnitude smaller.
This result can be rationalized as followings. Because of the slowness of the motion of the headgroups, the electrostatic potential of the headgroups fluctuates on the time scale of several nanoseconds. In response to these fluctuations, the water molecules reorient themselves and compensate for the fluctuating field. However, this reorientation of the water molecules happens on a much faster time scale, so that the potential is balanced overall. It should be emphasized that the fact that the total potential stays constant over time is not a result of the use of the Ewald summation, because the individual components fluctuate very strongly.
The values for the diffusion coefficients we obtained from our simulations are in reasonable agreement with the values obtained from the neutron scattering experiments given above (Pfeiffer et al., 1989; Konig et al., 1992; Sackmann, 1995). We compare our results with results from neutron scattering experiments, because the time scale as well as the length scales match those of molecular simulations. Most other techniques cover either much larger time scales and/or much larger length scales. The larger time or length scales could possibly cover phenomena not included in our simulation.
To get a better feel for the dynamics on different time scales, we looked at "stroboscopic pictures" of the trajectories of the COM of the phospholipid molecule every 4 ps for a time interval of 100 ps and a time interval of 1 ns and every 40 ps for a time interval of 10 ns. This is done in Fig. 6, where trajectories are shown for the time interval of 100 ps, 1 ns, and 10 ns.
As we can see from the top panel of Fig. 6, over 100 ps, the COM are not displaced very far from the initial positions, such that the motion of the COM in the lateral direction over this time period is reminiscent of a glassy state. Even during the 1-ns time period we do not observe a large displacement of COMs from their initial positions (Fig. 6, middle panel). The picture changes quite substantially in the bottom panel of Fig. 6, where the trajectories of the COMs clearly show that the membrane layer is liquid. This view is supported by Fig. 5, which shows that some lipids moved quite substantially during the course of the simulation, whereas others remained close to their initial position, which is the behavior expected for Brownian motion. On average we observed that the particles moved ~12 Angstrom over the course of the simulation, a distance larger than the next-neighbor distance.
The displacements in the z directions are predicted by the so-called protrusion model of Israelachvili and Wennerstrom (1990). This model is an attempt to explain the strong repulsive interaction between membranes over short distances. Experimentally it has been observed (LeNeveu et al., 1976; Rand and Parsegian, 1989) that at small separations, a repulsive interaction between opposing bilayer surfaces varies exponentially with distance. In the model of Israelachvili and Wennerstrom (1990) the movements of individual lipid molecules out of the membrane planes, socalled protrusions, are responsible for this force. If two membranes approach each other, the movement of individual lipid molecules out of the membrane plane is restricted by the presence of the lipid molecules of the approaching membrane. The presence of the opposing membrane therefore limits the configurational space of the lipid molecules. Assuming that the free energy of a protrusion increases linearly with its height and then using a Boltzmann ansatz for the height distribution lead to an exponentially varying repulsive force.
The above-mentioned numbers for the movements out of the membrane plane look rather short at first. However, as has been shown in a previous publication (Essmann et al., 1995a), even at a hydration level of 20.5 water molecules/ lipid, individual headgroups come rather close. At this hydration level the width of the water slab according to the definition of McIntosh and Simon (1986) is ~12 Angstrom, and a detailed analysis of the shortest distance between lipid molecules of opposing bilayers yields a similar length. At this distance the movement of the headgroup is sufficient to bring the nitrogen atoms in close proximity. Therefore the headgroups of the apposing bilayers start lo restrict each other's available conformational space. At this point it is not possible to decide whether this effect is the cause of the hydration force, because the free energies involved in the hydration force are very small compared to the free energies measurable in molecular simulations.
OVERALL ROTATIONAL DYNAMICS
Another very important characteristic of the dynamics of lipid molecules in membranes is the rotational motion. This has been measured by NMR (Blume, 1993), electron spin resonance (Lange et al., 1985), and dielectric spectroscopy (Klosgen et al., 1996). Because of the nature of the lipid bilayer, the overall motion can be characterized by a rotation around an axis parallel to the bilayer normal and a so-called wobbling motion in which the long axis changes its orientation.
Because of the great internal flexibility of the lipids, it is not possible to define unambiguously the body fixed axes along which we can calculate the rotation of lipid molecules. To define the axes we first calculated the moments of inertia and the principal axes for each lipid molecule (Goldstein, 1980). The principal axis that makes the smallest angle with the bilayer normal is called the z axis. If the lipid molecules have two parallel Snl and Sn2 chains as often sketched in textbooks, then this would also be the axis associated with the smallest moment of inertia. However, in practice we observed that sometimes one of the other axes was associated with a smaller moment of inertia. We attribute this to conformations in which the two alkane chains were "spread apart." Because of the great internal flexibility we did not use the principal moments to define the y axis, but rather defined them as the normalized perpendicular component of the vector from the center of the Snl to the center of the Sn2 chain. The x axis is the unit vector forming a right-handed axis system.
SUMMARY
The long-time run of the present study demonstrates that the constant volume ensemble leads to stable simulations of membrane-water systems. As long as the geometry of the simulated system is chosen properly, the strong tendency of the water molecules to form hydrogen bonds is sufficient to maintain a stable water-membrane interface. The limitation of this approach is that only those systems can be studied for which enough experimental data are available to determine the geometry of the simulation cell.
Even though the internal flexibility leads to intramolecular fluctuations on a time scale of nanoseconds, certain quantities such as the dipole potential are rather stable and can be computed with good accuracies from much shorter simulations.
We thank Prof. A. Geiger for helpful discussions, Dr. T. Darden of National Institute of Environmental Health Science for providing the PME code, and Dr. M. Crowley for the T3DT3E version of the PME code. Furthermore, we thank Dr. K. Schweighofer for carefully reading the manuscript. The simulations were performed at the North Carolina Supercomputer Center and the HLRZ Supercomputing Center, Germany. Parts of this work have been conducted within the GMD research group Computational Methods in Chemistry of the High-Performance Computing Center (HLRZ).
This work was supported in part by the National Science Foundation and the National Institutes of Health.
[Reference]
REFERENCES
[Reference]
Berne, B. J. 1971. Time-dependent properties of condensed media. In Physical Chemistry, An Advanced Treatise, Vol. VIIlb. H. Eyring, D. Henderson, and W. Jost, editors. Academic Press, New York. 539-716.
Blume, A. 1993. Dynamic properties. In Phospholipids Handbook. G. Cevc, editor. Marcel Dekker, New York. 455-509.
Brockman, H. 1994. Dipole potential of lipid membranes. Chem. Phys. Lipids. 73:57-79.
Cevc, G. 1993. Appendix B: thermodynamic parameters of phospholipids. In Phospholipids Handbook. G. Cevc, editor. Marcel Dekker, New York. 939-956.
[Reference]
Darden, T., D. York, and L. Pedersen. 1993. Particle mesh Ewald: an N-log(N) method for Ewald sums. J. Chem. Phys. 98:10089-10092.
Essmann, U., L. Perera, and M. L. Berkowitz. 1995a. The origin of the hydration interaction of lipid bilayers from MD simulation of dipalmitoylphosphatidylcholine membranes in gel and liquid crystalline phases. Langmuir. I 1:4519-4531.
Essmann, U., L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. Pedersen. 1995b. A smooth particle mesh Ewald method. J. Chem. Phys. 103:8577-8593.
Feller, S. E., Y. Zhang, and R. W. Pastor. 1995. Computer simulation of liquid/liquid interfaces. II. Surface tension-area dependence of a bilayer and monolayer. J. Chem. Phys. 103:10267-10276.
Frisch, M. J., G. W. Trucks, M. Head-Gordon, P. M. W. Gill, M. W. Wong, J. B. Foresman, B. G. Johnson, H. B. Schlegel, M. A. Robb, E. S. Replogle, R. Gomperts, J. L. Andres, K. Raghavachari, J. S. Binkley, C. Gonzalez, R. L. Martin, D. J. Fox, D. J. Defrees, J. Baker, J. J. P. Stewart, and J. A. Pople. 1992. Gaussian 92, Revision E.2. Gaussian, Pittsburgh, PA.
[Reference]
Gawrisch, K., D. Ruston, J. Zimmerberg, V. A. Parsegian, R. P. Rand, and N. Fuller. 1992. Membrane dipole potentials, hydration forces, and the ordering of water at membrane surfaces. Biophys. J. 61:1213-1223.
Goldstein, H. 1980. Classical Mechanics, 2nd Ed. Addison-Wesley, Reading, MA.
Israelachvili, J. N., and H. Wennerstrom. 1990. Hydration or steric forces between amphiphilic surfaces? Langmuir. 6:873-876.
Jacobson, K., E. D. Sheets, and R. Simson. 1995. Revisiting the fluid mosaic model of membranes. Science. 268:1441-1442.
Klosgen, B., C. Reichle, S. Kohlsmann, and K. D. Kramer.1996. Dielectric spectroscopy as a sensor of membrane headgroup mobility and hydration. Biophys. J. 71:3251-3260.
Konig, S., W. Pfeiffer, T. Bayerl, D. Richter, and E. Sackmann. 1992. Molecular dynamics of lipid bilayers studied by incoherent quasi-elastic neutron scattering. J. Phys. II. 2:1589-1615.
Lange, A., D. Marsh, K.-H. Wassmer, P. Meier, and G. Kothe. 1985. Electron spin resonance study of phospolipid membranes employing a comprehensive line-shape model. Biochemistry. 24:4383-4392.
LeNeveu, D. M., R. P. Rand, and V. A. Parsegian. 1976. Measurement of forces between lecithin bilayers. Nature. 259:601-603.
[Reference]
Lynden-Bell, R. M., and A. J. Stone. 1989. Reorientational correlation functions, quaternions and Wigner rotation matrices. Mol. Simulation. 3:271-281.
Marrink, S.-J., M. L. Berkowitz, and H. J. C. Berendsen. 1993. Molecular dynamics simulation of a membrane/water interface: the ordering of water and its relation to the hydration force. Langmuir. 9:3122-3131.
Mayer, C., G. Grobner, K. Muller, K. Weisz, and G. Kothe. 1990. Orientation-dependent deuteron spin-lattice relaxation times in bilayer membranes: characterization of the overall lipid motion. Chem. Phys. Lett. 165:155-161.
McIntosh, T. J., and S. A. Simon. 1986. Hydration force and bilayer deformation: a reevaluation. Biochemistry. 25:4058-4066.
Mouritsen, O. G., and K. Jorgensen. 1997. Small-scale lipid-membrane structure: simulation versus experiment. Curr. Opin. Struct. Biol. 7:518-527.
[Reference]
Nagle, J. F., R. Zhang, S. Tristram-Nagle, W. Sun, H. Petrache, and R. M. Suter. 1996. X-ray structure determination of fully hydrated La phase DPPC bilayers. Biophys. J. 70:1419-1431.
Nordera, P., M. Dalla Serra, and G. Menestrina. 1997. The adsorption of Pseudomonas aeruginosa exotoxin A to phospholipid monolayers is controlled by pH and surface potential. Biophys. J. 73:1468-1478.
Pastor, R. W. 1994. Molecular dynamics and Monte Carlo simulations of lipid bilayers. Curr. Opin. Struct. BioL 4:486-492.
Pastor, R. W., and S. E. Feller. 1996. Time scales of lipid dynamics and molecular dynamics. In Biological Membranes: A Molecular Perspective from Computation and Experiment. K. M. Merz and B. Roux, editors. Birkhauser, Boston. 3-29.
Pearlman, D. A., D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, D. M. Ferguson, G. L. Seibel, U. C. Singh, P. K. Weiner, and P. A. Kollman. 1995. AMBER, Version 4.1. Department of Pharmaceutical Chemistry, University of California, San Francisco.
[Reference]
Pfeiffer, W., T. Henkel, E. Sackmann, W. Knoll, and D. Richter. 1989. Local dynamics of lipid bilayers studied by incoherent quasi-elastic neutron scattering. Europhys. Lett. 8:201-206.
Pickar, A. D., and R. Benz. 1978. Transport of oppositely charged lipophilic probe ions in lipid bilayer membranes having various structures. J. Membr. Biol. 44:353-376.
Rand, R. P., and V. A. Parsegian. 1989. Hydration forces between phospholipid bilayers. Biochim. Biophys. Acta. 988:351-376.
Ruocco, M. J., and G. G. Shipley. 1982. Characterization of the subtransition of hydrated dipalmitoylphosphatidylcholine bilayers. Biochim. Biophys. Acta. 691:309-320.
Ryckaert, J.-P., and A. Bellemans. 1975. Molecular dynamics of liquid n-butane near its boiling point. Chem. Phys. Lett. 30:123-125.
Ryckaert, J.-P., G. Ciccotti, and H. J. C. Berendsen. 1977. Numerical integration of the Cartesian equation of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comp. Phys. 23: 327-341.
Sackmann, E. 1995. Physical basis of self-organization and function of membranes: physics of vesicles. In Handbook of Biological Physics, Vol. lA, Structure and Dynamics of Membranes. R. Lipowsky and E. Sackmann, editors. Elsevier, Amsterdam. 213-304.
[Reference]
Saffman, P. G., and M. Delbruck. 1975. Brownian motion in biological membranes. Proc. Natl. Acad. Sci. USA. 72:3111-3113.
Shinoda, W., N. Namiki, and S. Okazaki. 1997. Molecular dynamics study of a lipid bilayer: convergence, structure and long-time dynamics. J. Chem Phys. 106:5731-5743.
Singer, S. J., and G. L. Nicolson. 1972. The fluid mosaic model of the structure of cell membranes. Science. 175:720-731.
Szabo, A. 1984. Theory of fluorescence depolarization in macromolecules and membranes. J. Chem. Phys. 81:150-167.
Tieleman, D. P., S. J. Marrink, and H. J. C. Berendsen. 1997. A computer perspective of membranes: molecular dynamics studies of lipid bilayer systems. Biochim. Biophys. Acta. 1331:235-270.
Tobias, D. J., K. Tu, and M. L. Klein. 1997. Atomic-scale molecular dynamics simulations of lipid membranes. Curr. Opin. Colloid Interface Sci. 2:15-26.
Zannoni, C. 1994. An introduction to the molecular dynamics method and to orientational dynamics in liquid crystals. In The Molecular Dynamics of Liquid Crystals. G. R. Luckhurst and C. A. Veracini, editors. Kluwer, Dordrecht. 139-169.
[Author Affiliation]
Ulrich Essmann* and Max L. Berkowitz#
*GMD-German National Research Center for Information Technology, D-53754 St. Augustin, Germany, and #Department of Chemistry, University of North Carolina, Chapel Hill, North Carolina 27599 USA
[Author Affiliation]
Received for publication 27 April 1998 and in final form 5 January 1999.
Address reprint requests to Dr. Ulrich Essmann, GMD-SCAI, German National Research Center for Information Technology. Schloss Birlinghoven, D-53754 St. Augustin, Germany. Tel.: 49-2241-14-2795; Fax: 49-2241-14-2656; E-mail: essmann@gmd.de.
