INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

Biological membranes are an essential part of life processes in living organisms. Membrane function goes beyond its obvious role as a physical barrier to contribute to important life processes. Membranes are semipermeable, highly selective barriers containing ion channels and pumps to modulate and maintain balance as required. These channels and pumps are involved in signal transduction and response to stimuli from the environment. In addition, important energy conversion and storage processes take place at the membrane level. Consequently, a fundamental understanding of bilayers and membrane proteins from the atomic point of view is of great biochemical, biophysical, and medical interest (Mouritsen and Jørgensen, 1997; Pastor and Feller, 1996; Jacobson et al., 1995). A better understanding of the structure and dynamics of membranes and membrane proteins can contribute to the development of pharmaceuticals, anesthetics, and drug-delivery agents. In addition, studies of pure component-membrane structural and dynamical properties at the atomic level can enhance our understanding of more complicated biological membrane functions and their environmental interactions.

In this light, we turn our focus to pure single-component membranes that have been studied both experimentally and theoretically for their dynamic and structural characteristics. The time scales of biochemically relevant fluctuations in membranes can span anywhere from femtoseconds (10¹⁵ s), for intramolecular vibrations, to minutes or hours (10³ s), for lipid molecule transbilayer flips (Blume, 1993). Lipid molecules in large membranes are believed to assemble and move collectively as aggregates (so called "rafts"), which can span several hundred angstroms of the bilayer surface (Schütz et al., 2000).

Important transport properties in membranes, some believed to exhibit anomalous behavior, are currently poorly understood even for single-component membranes at the atomic level (Granek, 1997; Granek and Pierrat, 1999; Shlesinger et al., 1999; Schütz et al., 2000; Schwille et al., 1999). These dynamic properties include lateral diffusion of lipids within a membrane, dynamical fluctuations within a lipid moiety (rotational, reorientation, and relaxation), gauche-trans equilibrium, and their respective interconversion to name a few. A quantitative understanding of such processes at the molecular level has yet to be achieved, and simulations of the type presented in this paper aim to provide atomistic detail, which will possibly lead to a better understanding of such properties.

Molecular Dynamics (MD) calculations can be a powerful tool to probe structural and dynamical properties in membranes at the atomic level. However, two limitations exist when one uses MD to probe membrane properties. The first limitation is system size. Typical MD simulations seldom span more than 100 Å in any direction of the simulation cell, making it feasible to only treat a few hundred lipid molecules with present-day resources. The second limitation involves the accessible time scales. Current MD simulations seldom span more than a few nanoseconds, with simulations of 10 or more nanoseconds being relatively rare for this type of system. The longest/largest simulation of a lipid that has been reported to date (Lindahl and Edholm, 2000) extends these factors by an order of magnitude. However, certain simplifying approximations were necessary to enable this calculation. Even so, the scales for time and space are still small compared to the time and length scales for some important lipid membrane phenomena. We therefore limit ourselves to study membrane properties, which are relatively fast and localized. Fortunately, many important properties such as isomerization, rotational relaxation, and lateral diffusion fall within the current limits. The present study is an attempt to capture and elucidate some of the structural and dynamical fluctuations of a pure component lipid on the order of several nanoseconds.

We have chosen to simulate a pure dimiristoylphosphatidylcholine (DMPC) bilayer, which has a phosphatidylcholine head group with two 14-carbon atom alkane chains (Fig. 1). DMPC has been studied experimentally by electron spin resonance (ESR), nuclear magnetic resonance (NMR), infrared absorption (IR), and fluorescence recovery after photobleaching (FRAP), among others for structure and dynamics (for a review, see Cevc, 1993). Recent studies of DMPC include x-ray diffraction (Petrache et al., 1998; Movromoustakos et al., 1990), NMR (Marsan et al., 1999; Trouard et al., 1999; Hong et al., 1996), and computer simulation (Zubrzycki et al., 2000; Kothekar, 1996; Damodaran and Merz, 1994; Pasenkiewicz-Gierula et al., 1999; Duong et al., 1999) for structural and dynamical properties. Similar lipids (such as DPPC (Essmann and Berkowitz, 1999; Lindahl and Edholm, 2000; Venable et al., 2000; Shinoda et al., 1997), DPhPC (Husslein et al., 1998), DLPE (Damodaran and Merz, 1994), POPC (Chui et al., 1999), and DOPC (Chui et al., 1999) to name a few) have also been studied by MD simulation, but, to our knowledge, no multinanosecond study of DMPC bilayers with an all-atom force field has been reported. It has recently become common to truncate the long-range electrostatic interactions to reduce computational effort. This truncation is done in spite of compelling evidence that long-range electrostatic interactions are very important for both the dynamics and structure of molecules (Feller et al., 1996; Norberg and Nelsson, 2000). Our calculations do not truncate the electrostatic interaction. Instead, we use the Ewald summation technique to incorporate the full electrostatic interactions.

View larger version (17K): [in this window] [in a new window]   FIGURE 1   The chemical bonding representation of the DMPC molecule. Carbon atom labeling is referenced throughout the text. The head group Zwitterion dipole is indicated by the P-N vector.

In the following sections, we present the key results of a 10-ns all-atom MD simulation of DMPC in the L (liquid crystalline) phase. The remainder of the present work is organized as follows. In the next section, we describe the methods used, followed by a discussion of the initial conditions and setup in the following section. We then focus on the stability of our calculations, the robustness of our results, and the average structural quantities of DMPC from our simulation. The next section turns to dynamical properties that have been calculated from the simulation, and the work concludes in the last two sections.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

We have developed and used a state-of-the-art program to perform classical MD simulations on any system for which a classical force field exits. The Center for Molecular Modeling Molecular Dynamics (CM³D) program allows for simulations to be run under a wide variety of ensembles, ranging from constant energy (NVE) to fully-flexible constant pressure and constant temperature (NPT) (Moore and Klein, 1997).

MD simulations can provide an accurate description of physical phenomena by integrating the differential classical equations of motion in a computer (Allen and Tildesley, 1987). The choice of time step to accurately reproduce the dynamics of the system being treated is governed by the fastest degrees of freedom, which, in most cases, are the intramolecular bond vibrations. This limitation has been relieved by the implementation of multiple time-step integrators (Tuckerman et al., 1992; Humphreys et al., 1994; Martyna et al., 1996; Cheng and Merz, 1999; Tuckerman, 2000). The reversible REference System Propagator Algorithm (RESPA) (Tuckerman et al., 1992), integrates the fast motions with small time steps and slow motions with larger time steps, resulting in an effective decrease of wall clock CPU time in the simulation. Our implementation of reversible integrators stems from an approximation to the Liouville operator (Trotter, 1959; Martyna et al., 1996). These integrators can handle the different classical simulation ensembles with extended system approaches (Evans and Holian, 1985; Tobias et al., 1993; Martyna et al., 1996).

We used periodic boundary conditions to minimize the finite size effects on the system, and no constraints to internal degrees of freedom were included. We have used the AMBER force-field (Cornell et al., 1995), which has been shown to provide a good description of lipid properties (Tieleman et al., 1997). Due to our choice of conditions, the results of this simulation should correspond to a DMPC lipid bilayer in the liquid crystalline (L) phase.

Taking advantage of the previously mentioned algorithms, we have used a 4-fs time step, which is 8 times greater than the time step required to integrate the fastest bond vibrations accurately (on the order of ~0.5 fs). As will be discussed in section after next, this allows for stable trajectories with no need for rescaling or resampling of velocities, and a great reduction in wall clock CPU time.

We truncate the short-range van der Waals interaction at 12 Å. It has been shown that long-range electrostatic interactions can be extremely important to the structural and dynamical properties of charged systems (Feller et al., 1996; Norberg and Nelsson, 2000). The long-range electrostatic interactions were evaluated using either Ewald (Brown and Neyertz, 1995; Nymand and Linse, 2000), or particle mesh Ewald (PME) (Darden et al., 1993; Essmann et al., 1995) summations techniques. We have also implemented a parallel scheme, the replicated data technique (Wilson et al., 1997), which effectively uses multiprocessor computers with 4-64 nodes. All the above techniques have been included within CM³D. This choice of integrator, parallelization, cut-off, Ewald and RESPA scheme allows for extremely stable trajectories over the full 10-ns trajectory, which we describe in more detail in the section after next.

The present simulation has been run at the Pittsburgh Supercomputing Center on 32 nodes of the Cray T3E, 16 Nodes of an SGI origin 2000 at NCSA, and on local machines. The total simulation time resulted in approximately 60,000 CPU hours of single-processor computer time.

	INITIAL CONDITION

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

The simulation of the lipid and water system was constructed by replicating a single-lipid molecule into an idealized eight-molecule lipid bilayer and then adding water to "cap" the top and bottom, making an effective bilayer with periodic boundary conditions. This eight-molecule lipid bilayer was then run with constant volume and constant temperature MD (NVT) for a few picoseconds. This small simulation was performed to relax the structure from nonphysical contacts created by the construction of the bilayer and water cap. The relaxed eight-lipid configuration was then replicated to produce 16 lipids, which was again allowed to relax for a few picoseconds. This procedure was performed once again for 32 lipids and finally for 64 lipids to obtain the initial bilayer coordinates, which consisted of 64 DMPC and 1792 water molecules for a total of 12,928 atoms. The ratio of water to lipid molecules (n_w) was set at 28, allowing for proper hydration of the lipid bilayer (Petrache et al., 1998).

This 64-lipid bilayer was then run for 500 ps using constant-volume NVT and then 500 ps more using constant pressure, constant temperature (NPT) to equilibrate the simulation cell. Both the NVT and NPT calculations were carried out at a temperature of 333 K. Our criterion for equilibrium of the simulation cell dimensions was the observation of fluctuations of each box side around some stable average value for more than 100 ps. Once an average size for the simulation cell was established, we transferred the system to the NVE ensemble, taking the average cell lengths as the fixed box lengths to collect dynamical data. The final relaxed cell dimension were 41.7 × 44.2 × 68.84 Å, and an average temperature of 333 K.

Because we want to extract dynamical data as well as structural data, we use the NVE ensemble. The dynamics of NVT or NPT ensemble will be influenced by the coupling of the necessary extended variable and could give rise to errors in the dynamics results. The error associated with the external variable is of order 1/, where N_f is the number of degrees of freedom to which the extended variables are coupled.

This method to build the simulation cell using NVT and NPT and then transfer to NVE has been shown to provide suitable results for dynamics of membranes in the liquid phase (Essmann and Berkowitz, 1999). All the quantities that we observed (energy, temperature, total pressure, individual components of the pressure tensor, etc.) had no detectable long timescale instabilities.

	SYSTEM STABILITY

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

In this section, we examine the question of multinanosecond stability of the trajectory. It has been shown in other studies that simulations involving membrane systems are extremely susceptible to starting conditions (Husslein et al., 1998; Tu et al., 1996; Tu, 1995). To ensure that our system is both at equilibrium and that it behaves properly, we take great care that conserved quantities fluctuate around a constant value with no drifts throughout the simulation.

It has been reported that long time-scale simulation can become unstable and, in particular, that electrostatic truncation can influence the stability of the simulation (Norberg and Nelsson, 2000). As we stated in the second section, we do not truncate the electrostatics. As seen in Fig. 2, the total energy is conserved, does not drift. The trajectory is stable throughout our entire simulation (the embedded graph in Fig. 2 shows the total energy fluctuates three orders of magnitude less than kinetic energy in a 1-ns time span). Furthermore, the temperature (proportional to the average kinetic energy) remains stable throughout the simulation at ~333 K. It is noteworthy that no constraints or temperature controls were used during the simulation, thus proving that a careful equilibration can lead to stable trajectories with lipid bilayer systems.

View larger version (29K): [in this window] [in a new window]   FIGURE 2   Kinetic energy, intramolecular potential energy (PEA), and intermolecular potential energy (PEI) components in the NVE ensemble. All components fluctuate less than 0.01% around well-defined average values and add up to the total energy (TE) line. The embedded graph shows an expansion of the variation of TE with respect to time for 1 ns, showing very slight variations in total energy.

The present results further confirm and extend those of Berkowitz and coworkers in a previous work concerning the stability of membrane systems in a fixed box geometry (Shinoda et al., 1997; Essmann and Berkowitz, 1999). Using a reasonable choice of initial parameters, the system is stable for many nanoseconds, and, with proper integration of the equations of motion, there is no further need for rescaling or temperature control.

	STRUCTURE

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

In this section, we examine some of the structural information obtained in the MD simulation. We calculate the electron density along the bilayer normal, and the alkyl chain, deuterium-order parameters and compare them to previous experiment and calculations. In general, we find good agreement in both of these quantities. The details of these calculations are discussed below.

Our procedure to simulate the system in the NVE ensemble involves, as discussed previously, relaxation in the NVT ensemble followed by further relaxation in the NPT ensemble. Once the system is transferred to the NVE ensemble, the volume of the box fixes the area per head group (A_H). The relaxed cell dimensions yielded A_H = 57.7 Å², which is in fair agreement with experimental results (Petrache et al., 1998) and other computations (Damodaran and Merz, 1994; Duong et al., 1999) summarized in Table 1. The value of A_H arises naturally from the relaxation of the system at fixed temperature, pressure, and number of molecules rather than a choice or fit to experiment, thereby showing that both the forcefield and the approach are valid. The average pressure tensor was ~1 atm along the diagonal and essentially zero in the off diagonal, indicating that there was no sheer stress on the system.

To test the robustness of our approach with respect to the starting conditions, we ran two additional 500-ps simulations using the NPT ensemble after 5 and 10 ns of simulation. However, these simulations merely reproduced the initial values, including box size, and energy, among others. We therefore conclude that our chosen NVE box was of suitable dimensions.

To further explore the equilibrium properties of our simulation, we have analyzed the lipid structure by calculating the gauche-trans distribution as a function of bond position along the tails, as shown in Fig. 3. The distribution varies along the length of the hydrocarbon chain, and it shows that the lipid tails are more likely to be in the trans conformation close to the head group, thus becoming more disordered further down the chain. This is in agreement with experiment and previous calculations (Urbina et al., 1995; Trouard et al., 1999; Tu, 1995; Shinoda et al., 1997).

View larger version (21K): [in this window] [in a new window]   FIGURE 3   Two representative torsional angle distributions for B1 (Torsion #1) and B11 (Torsion #11), where the bond labels correspond to those in Fig. 1.

We define the tilt angle as the angle between the lipid molecule principal moments of inertia and the bilayer normal. We find the average lipid tilt angle is 0 degrees. The distribution of tilt angles is fairly broad, demonstrating that the lipids move quite readily and can have tilt angles anywhere from 0 to 50 degrees, with similar profiles to that of Shinoda et al. (1997) for DPPC. We will expand on the dynamics of this "wobble" in the next section.

Electron density profile

We use the electron density profile along the bilayer normal, both as a comparison to experimental data and other calculations. Experimentally, the electron density profile (z) across an interface is related to the normalized x-ray reflectivity (R(q_z)/R_F(q_z)) by the equation,
where R_F(q_z) is Fresnel reflectivity, is the electron density of the semi-infinite bulk sub-phase and q'_z is the photon momentum transfer perpendicular to the interface surface (Zheng et al., 2001; Als-Nielsen and Pershan, 1983; Helm et al., 1991). Figure 4 reports the electron density profile calculated from our simulation, with contributions of individual selected atoms in the bottom part of the graph. The experimental electron density (at 27 atm) is also included in this graph. The zero of the x-axis has been set approximately at the center of the bilayer along the z-axis, where a decrease in electron density (the so called "methyl trough") is seen in the lipid tails. Slight statistical differences remain in the electron densities of each monolayer. We therefore symmetrized our results.

View larger version (28K): [in this window] [in a new window]   FIGURE 4   The electron density profile along the bilayer normal (z-axis) for individual atoms: choline nitrogen (×), the phosphate group phosphorus (+), the non-ester oxygen atoms both in the phosphate and in the carbonyl groups (), the carbon atoms in the carbonyl groups (- - -), water (), and the total lipid electron density (· · · · ·) are shown along with the total electron density (). The experimental values from Petrache et al. (1998) are displayed in as a thin solid line with no label.

The structure of the bilayer from the electron density is comparable to that reported previously (Nagle et al., 1996; Petrache et al., 1998). We also find good agreement with the corresponding x-ray diffraction values previously reported by Petrache et al. (1998) and those estimated by previous calculations (Zubrzycki et al., 2000).

For a direct comparison, we have digitized and overlayed the experimental x-ray diffraction by Petrache et al. (1998) in Fig. 4. The major discrepancy from our calculated electron density profile and that obtained by x-ray diffraction is that there is more electron density in the region of the head group, and a slight decrease in electron density in the middle of the bilayer. In the head group region, the extra width that we obtain could indicate that our interface is less ordered and more diffuse (given rise to a broader and less well-defined peak at the interface) than that in the experimental results. This might be due, in part, to the pressure difference between experiment (27 atm) and our calculation (~1 atm). These differences could also be caused by the potential energy functions we use, which will influence both the packing at the interface and the pressure that we calculate. Although there are slight differences in the magnitude of the electron density the bilayer thickness is the same. The D spacing that is derived from the electron density and those of other groups are quoted in Table 1 for comparison.

Order parameter

The deuterium-order parameter (S_CD) obtained from NMR gives a measure of the average methylene group orientation with respect to the bilayer normal (Seelig, 1977; Bloom et al., 1978; Meier et al., 1986; Urbina et al., 1995, 1998; Damodaran and Merz, 1994). The NMR experiments on phospholipids can be performed with site-selective labeling, thus providing detail on the chain-order parameter and a link to the present MD results. In simulations, S_CD is obtained as

(1)

where  is the angle between a vector normal of the plane formed by the Carbon (C) and Deuterium (D) atoms and a vector normal to the bilayer. The S_CD is obtained experimentally from the axially symmetric electric field splittings obtained from powder patterns. The quadrupole splittings are related to S_CD by

(2)

where e²qQ/h is the quadrupole coupling constant. Assuming that substitution of deuterium by hydrogen results in the same structure, we report the order parameter for each carbon along the chain in Fig. 5. From the 10-ns simulation, the uncertainty in the values corresponds to less than the thickness of the lines in Fig. 5.

View larger version (30K): [in this window] [in a new window]   FIGURE 5   Hydrocarbon deuterium order parameters, SCD as a function of the carbon position along the Sn-1 (- - -) and Sn-2 (· · · · ·) chain. The curve labeled with () denotes the average of both chains. The other lines represent experimental data. The curve labeled () is the Sn-2 order parameter determined by Urbina et al. (1995). The () label corresponds to the Sn-2 and the () to the Sn-1 order parameters determined by Trouard et al. (1999).

The present results agree fairly well with both experiments (Urbina et al., 1998; Trouard et al., 1999) and previous calculations (Damodaran and Merz, 1994). There are slight differences in each chain, suggesting that one of the chains is slightly more ordered than the other. This difference is also seen in the experiment of Trouard et al. (1999), where they see very similar order parameters for the Sn-1 and Sn-2 chains. Another interesting feature in the order parameters is a dip at carbon 3 on the Sn-1 chain, which agrees with experimental observations in DPPC (Seelig, 1977; Urbina et al., 1995, 1998). We find that the third carbon is, in fact, more disordered, in agreement with these experiments; a finding confirmed by the fraction of trans conformers displayed in Fig. 6. Not surprisingly, the S_CD trend is reproduced in the distribution of trans conformers, because the lipid tails are oriented with respect to bilayer normal; the more gauche defects the more disordered the chains will become. Figs. 5 and 6 confirm that the order parameter is linked to the gauche-trans distribution along the chain. This observation of the disorder of the third methylene carbon is, however, not observed in the data reported by Trouard et al. (1999) for their estimate of S_CD.

View larger version (20K): [in this window] [in a new window]   FIGURE 6   Fraction of trans conformers as a function of bond number as labeled in Fig. 1. The Sn-1 chain is denoted by () and the Sn-2 by ().

	DYNAMICS

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

In this section, we examine the calculated dynamical properties of the DMPC bilayer. We proceed by examining the lateral diffusion and the rotational relaxation. We begin with a qualitative picture of the dynamics of lipid motions given in Figs. 7 and 8. Figure 7 shows the "stroboscopic" picture of the movements of eight selected lipid center of masses (COM) projected onto the XY-plane. The COM trajectories of most lipids are akin to that of molecules in a 2d solid. To better understand the dynamics of the lipids, we have superimposed the motion of a single lipid on different time scales in Fig. 8. On the 100-ps time scale, one observes the intramolecular vibrations and few torsional gauche-trans bond flips. On the 1-ns time scale, we see that the lipid has increased amplitude motion. The "smearing" is due to translation and rotation of the lipids and torsional gauche-trans flips in the tail region. However, the head group has not moved very much and is still oriented in roughly the same direction. We will quantify this in the Rotational diffusion subsection. On the 10-ns time scale, we see that the lipid has started to diffuse laterally, there is large amplitude motion of the tail, and that the head group has undergone significant rotation.

View larger version (29K): [in this window] [in a new window]   FIGURE 7   The center of mass motion of seven selected lipid molecules projected onto the XY plane. The full 10-ns trajectory is shown. Only the lipid labeled with a circle exhibits a "jump."

View larger version (42K): [in this window] [in a new window]   FIGURE 8   Ten superimposed structure snapshots for different time intervals (left to right: 100 ps, 1 ns, and 10 ns).

Lateral diffusion

Lateral diffusion in lipid membranes has been studied for at least 40 years. However, the mechanism by which lipids diffuse is currently not well characterized or understood. Several theories exist (Shinoda et al., 1997; Devaux and McConnell, 1972) relating to their lateral diffusion. One might expect that this diffusion is similar to that of an ideal fluid on a two-dimensional surface. However, in the phase of interest (the L phase) the lipid is, in fact, a liquid crystal. Thus diffusion occurs by a combination of liquid-like 2d diffusion and cage hopping between locally confined regions in a liquid crystalline lattice (König et al., 1995). Previous simulation studies did not reveal two distinct time scales for the diffusion of lipids in a membrane bilayer (Essmann and Berkowitz, 1999). Even so, experiments suggest otherwise, and lipid diffusion rates have been reported ranging over two orders of magnitude, which suggest some type of anomalous effects (assuming they measure similar properties under similar conditions) (Blume, 1993; Wu et al., 1977; Schütz et al., 1997). The difference in these experimental results has been commonly accepted to be due to the localization of the relaxation and the time averaging that goes into calculating dynamical quantities from different experimental techniques. These theories also suggest different time scales for the lipid, which motivates a caging-type diffusion mechanism.

Figure 7 demonstrates the localization of the lipids, one of the trajectories (circled) is clearly different. This lipid trajectory has two distinct domains connected by a thin transfer path. A reasonable description of the motion is that the lipid exhibits a "jump" between the two sites. This type of diffusion will give rise to the reported caging-type mechanism. Unfortunately, we lack enough jump events in our simulation to quantify the diffusion mechanism. However, such events have not been previously reported in an MD simulation. Longer simulations could provide a clearer picture of the caging-type diffusion mechanisms.

From the previous discussion, there is at least qualitative evidence for both caging and 2d liquid-like diffusive behavior in our simulation. To quantify the lateral motion, we calculate the self-diffusion coefficient D, given by the slope of the average mean square displacement (MSD) at long times. Specifically, D is defined by

(3)

where d_f is the number of dimensional degrees of freedom in the system and r_i(t) is the COM position of the particle i at time t. We use Eq. 3 with d_f = 2 for lateral diffusion in a 2d space (the x-y plane).

The COM average MSD plots are shown in Fig. 9. After 2 ns, the MSD appears relatively constant and a linear fit in the 1.5-3.0-ns interval yields a diffusion coefficient of 6.5 ± 1.0 × 10⁸ cm²/s. This value is within the bounds of the experimental values previously reported using FRAP and NMR, namely D ranging from 2 to 8 [× 10⁸ cm²/s] (Cevc, 1993). If one uses the 300-500-ps interval slope, as has been done with previous simulations of other lipids (Essmann and Berkowitz, 1999), one finds a different diffusion coefficient of 23 ± 2.0 × 10⁸ cm²/s. Although the system has not yet reached the hydrodynamic limit, it is interesting to note that this value is comparable to results obtained with excimer techniques, which can probe diffusion on the subnanosecond time scale (see Table 2).

View larger version (24K): [in this window] [in a new window]   FIGURE 9   The in-plane (x-y) average mean square displacement for the center of mass of the lipids. The linear fits for 300-500 ps and 1.5-3 ns are shown with squares and circles, respectively.

From these results, we anticipate that cage hopping, as observed in the present work, occurs on a much longer time scale (5 ns or more). We only observe two jump events in our simulation, and longer simulations will be required to obtain a better description of this jump diffusion. The cage jumps and "domain" rearrangement, which we do not capture, could lead to anomalous subdiffusion, which has been seen experimentally (Shlesinger et al., 1999).

Rotational diffusion

To investigate the rotational diffusion of the DMPC molecule, we must first define a molecular fixed reference frame (Fig. 10). The principal axes of the moment of inertia tensor will be used to give us an idea of the overall rotation of the molecule. We also investigate a vector drawn from the P atom to the N atom of the head group (the P-N vector). This vector will give information about the local environment at the lipid-water interface. Similarly, we study the motion of a vector drawn from the COM of the Sn-1 tail to the COM of the Sn-2 tail. Finally we choose two vectors along the Sn-1 and Sn-2 chains drawn from the first to the last carbon of each chain. These last three vectors will provide information about the rotation properties inside the bilayer.

View larger version (22K): [in this window] [in a new window]   FIGURE 10   Illustration of the various vectors chosen to study rotational dynamics of DMPC. The Sn-1, Sn-2, P-N, and the intramolecular tail vectors characterize individual rotational components. The COM moment of inertia tensor provides a picture of the overall lipid rotation.

When using the moment of inertia tensor, we assign the eigenvector corresponding to the smallest eigenvalue as the long axis of the lipid (approximately along the bilayer normal). The other two eigenvectors are initially assigned arbitrarily, keeping a right-handed reference frame. These eigenvectors define the internal coordinate system of the lipid. At each consecutive time step we evaluate and overlap the eigenvectors such that the identity of each vector is preserved. In this way, we can follow the rotational character of the individual eigenvectors over time and analyze its correlations.

We have characterized the rotations by using the MSD of the various angular variables. To do this, we use Eq. 3 and replace r with . We also set d_f = 1 corresponding to the one degree of freedom. The  variable is defined as the angle between the projection of one of the vectors onto the xy-plane and the x axes. We choose the xy-plane projection to look at rotational properties that are perpendicular to the bilayer normal. Care must be taken in the treatment of  so that the variable samples the space from to + and not from and + (in radians) because this would give spurious results. Taking all these assignments into account, we can now calculate a rotational diffusion coefficient for the three principal axes of the inertia tensor, the two-tail vectors, the PN vector, and the intramolecular vector. From our experience, we have found this approach more robust than fitting the lipid molecules to models (for example, approximating the lipid to a cylinder) or using Wigner rotation matrices (for example, the C₂(t) correlation function (Pastor and Feller, 1996; Berne and Pecora, 1990; Lynden-Bell and Stone, 1989)), and fitting the results to a decaying exponential (which will yield a rotational relaxation rate). One reason for the uncertainty in using Wigner rotation matrices is that the molecules are not modeled well by a cylinder. The lipid molecule changes shape at the same time that it rotates, making it difficult to constrain the lipid to a geometrical construct. This can be seen in the variance of the moment of inertia, which fluctuates from a prolate to that of an oblate top, as well as an asymmetric ellipsoid.

We can see from Fig. 11 that each component of the lipid (as defined by Fig. 10) rotates at a different rate. We obtain excellent linear fits in the region of 0.5 to 3 ns; in Fig. 11, we truncate the graph at 2 ns for clarity (the rotational scales are not too different). The Sn-1 and Sn-2 tails have the same rotational diffusion of D_tail = 25 rad²/ns. The second fastest rotation D_axis = 12 rad²/ns is the vector corresponding to the long axis of the molecule, this is similar to a wobble projected onto the xy-plane. These two rotations are faster then the other rotation and have a small projection onto the XY-plane, thus it is uncertain that the projection onto the XY plane is a good description of these motions. However the fits do provide a qualitative picture that these motions are faster then the other rotational times.

View larger version (31K): [in this window] [in a new window]   FIGURE 11   Rotational diffusion. The mean squared rotational of the vectors illustrated in Fig. 10. Each vector is projected onto the XY-plane and the angle with respect to the X axis is correlated. The Drot values are based on the einstein relation. The lines are as follows: the P-N vector (· · · · ·), the lipid tails (- · -), the moment of initia tensor (MIT) (), the MIT wobble (- - -), the vector between C28 and C18 intra-tail (8) (- · -), and the vector between C214 and C114 intra-tail (14) (- · · -) (see Fig. 1 for labels).

The overall rotation of the lipid is determine by projecting the eigenvectors of the inertia tensor onto the xy-plane, and we obtain D_MIT = 3.2 rad²/ns. The rotation that corresponds to the P-N vector (D_PN = 2.2 rad²/ns) is only slightly less then the D_MIT, which implies that the head group and MIT are closely linked. This close association is expected because the P is close to the COM, and the PN vector is parallel to the MIT eigenvector corresponding to rotations of the lipid. However, the slightly faster PN vector demonstrates that there is a "crank-shaft" type motion of the headgroup, such that the overall rotation of the lipid is slightly slower than head group rotation.

The last rotational motion that we consider is that of a vector between the SN1 and SN2 chains of the tail. For the vector between the last carbons in SN1 and SN2, we obtain a rotational diffusion coefficient of D_IMT14 = 0.6 rad²/ns, whereas the vector between the tails at the eighth carbon has a rotation diffusion coefficient of D_IMT8 = 0.04 rad²/ns. These coefficients suggest that the ends of the tails are rotating faster with respect to one another than are the center of the tails. Both of these rotations are much slower then the MIT, which is expected because the direction of the eigenvectors is perpendicular to the tail's orientation. The slow tail vector rotation is due, in part, to frustration and entanglement of the tails.

Experimentally, rotational diffusion is difficult to measure, but, for DMPC, has been estimated to be in the nanosecond range for the relaxation time (Cevc, 1993). The diffusion coefficients of all vectors lie within this time scale. A direct measurement from Harms et al. (1999) reported a rotational diffusion coefficient of D_rot = .07 rad²/ns, for a fluorophore tagged in POPC molecule. This is in good agreement with the MIT rotation, but is 60 times slower than what we find in our simulation for the head group, which should be the proper comparison. The difference could be accounted for by the size of the (large) fluorophore (not the lipid itself) or the different lipid conditions, such as temperature and composition. Another discrepancy with experiments lies in the difference between the head group and the fluorophore rotation. The lipid head groups can rotate, twist, or arc, whereas only the rotational degree of freedom is accessible to the fluorophore. Others have also estimated the rotational relaxation time of lipids: _rot  1 ns using methods described by Pastor and Feller (1996) and _rot  6 ns using methods described by Essmann and Berkowitz (1999). If we take our average rotational velocity and calculate a relaxation time, these values fall between the wobble and the PN vector motion that we observe.

Our calculations provided a broad picture of the rotational motion of the lipid molecules. To sample the entire rotational configuration space (2 rad  40 rad²) with a rotational diffusion of 3.2 rad²/ns, the simulations time needs to be 10 ns or greater. As a result, only a minimum time scale for the rotational relaxation time of lipid molecules has been established. Future studies should focus on longer time scales to obtain the relaxation and rotational character of the lipid. Of course, the inclusion of other components into the membrane, such as other lipid molecules, proteins or cholesterol, will likely influence the overall rotational environment.

The present analysis, which projects the rotations onto a single angular coordinate, does not differentiate between pure rotation and diffusion involving rotation-translation coupling. However, the limited amount of translation diffusion contrasts with the extensive rotational motion. Thus, with the exception of the intramolecular vector, the lipids samples all the rotational phase space before they have time to translate or exchange positions with other lipids. Therefore, we conclude that the translation and most of the rotational motions are uncoupled under the conditions that we are studying. The intramolecular tail vector could show rotational-translational coupling, but we are unable to investigate this phenomenon because of the limited translational motion in our system.

	CONCLUSION

TOP ABSTRACT INTRODUCTION METHODS INITIAL CONDITION SYSTEM STABILITY STRUCTURE DYNAMICS CONCLUSION REFERENCES

We have demonstrated that a simulation greater than 10 ns can be stable in the NVE ensemble and hence the utility of MD to determine long-time physical properties of the lipid. Our calculated structural properties of the DMPC bilayer agree quantitatively with previous simulations and experiments. Longer simulations of 100 ns or more on larger systems could be insightful to investigate the lateral diffusion of lipids, unfortunately we do not have the computational resources at this time.

Our calculated translational diffusion coefficient is well within the experimental error, and we see evidence for caging and cage jumps of the lipids as well as liquid-like diffusion. The presence of both types of diffusion within a single lipid simulations has not been reported before. We also calculate the overall rotational diffusion and report a value of 3.2 rad²/ns, which suggests that, for the lipids to sample their cage surroundings, at least 10 ns are needed. The time scales for the rotation of individual lipid components (head and tail) span three orders of magnitude. This provides an interesting picture about the lipid motions; one where the lipid tails rotate faster than the head group, which in turn is rotating faster than the tail vector.