INTRODUCTION

TOP ABSTRACT INTRODUCTION A MESO-SCALE MODEL FOR... MESO-SCALE SIMULATION OF LIPID... THERMAL OSCILLATIONS RESULTS CONCLUSIONS APPENDIX REFERENCES

Many large biological assemblies inherently possess multiple length and time scales, resulting from the disparity in the dimensions of their structure. In the case of lipid bilayers, the width, h, of the bilayer exists in microscopic domains (nanometers), whereas the area, A, can persist up to nearly macroscopic length scales (micrometers) (Lipowsky and Sackmann, 1995; Tieleman et al., 1997a,b; Bagatolli and Gratton, 2000; Bagatolli et al., 2000; Forrest and Sansom, 2000). Thus, to completely model the structure and dynamics of such assemblies, it is necessary to span the entire regime from atomistic to macroscopic length scales. Concurrently, it is also necessary to examine larger and larger time scales as the fundamental length scale is increased. For example, orders of nanoseconds are required to examine the dynamics of a single molecule in liquid water; however, macroscopic times are required to model its continuum fluid flow properties.

We have developed a multiscale simulation method (Ayton et al., 2001a,b) whereby atomistic-level simulations can be bridged to continuum-level models by calculating the evolving macro-scale material properties from atomistic models. With this method, Green-Kubo Theory (see e.g., Evans and Morriss, 1990) defines the relationship between time-averaged correlations of microscopic quantities and the corresponding macroscopic transport coefficients. Thus, the basis of the technique relies on accurately calculating the specific transport coefficients or other material properties (depending on the nature of the material) from detailed molecular models. This information is then used within a constitutive relationship valid at longer spatial and temporal scales. In essence, beyond atomistic length and time scales, we abandon the notion of using a molecular representation as the fundamental unit of description, and instead we use more coarse-grained and time-averaged material properties as the governing representation of the system. An important component of this methodology beyond that which is usually used in nonequilibrium statistical mechanics is the feedback loop that is constructed from the macro-scale back to the atomistic-level simulations.

Our multiscale simulation method has been developed within the context of lipid bilayers (Ayton et al., 2002b). In this case, the bulk modulus of a dimyristoylphosphatidylcholine (DMPC) bilayer, , was calculated using detailed atomistic-level nonequilibrium molecular dynamics (NEMD) (Ayton et al., 2002a), and this quantity was then transferred to a continuum-level model of a giant unilamellar vesicle (GUV) (Ayton et al., 2002b). In formulating the continuum-level equations of motion for this application, it was found that only the bulk modulus was required to resolve the desired continuum-level dynamics. However, to correctly bridge these two length and time scales, it was necessary to constrain the system to be under a prestressed state. Only then could the modulus, as determined from the microscopic simulation, be related directly to the continuum representation. To be more specific, with our formulation of NEMD, system sizes on the order of the membrane thickness are modeled under periodic boundary conditions. The small system size implicitly inhibits any long-wavelength thermal-bending fluctuations. At the continuum level, undulation-free states are only obtained in prestressed states with a significant surface tension.

To examine flaccid membranes or membranes subject to other conditions, a simulation method must be adopted that operates in length and time scales where thermal fluctuations can occur and be accurately modeled. Experimentally, for example, it is found that the stress versus strain curve for a GUV is nonlinear in the low strain regime (Rawicz et al., 2000). The explanation for this behavior is proposed to be related to the existence of subvisible thermal ripples or bending modes (Evans and Rawicz, 1990, 1997). In this low-strain regime, the apparent area is less than the actual area because of thermal oscillations. The small value of the observed modulus is therefore a result of "undoing" these soft bending modes. After these modes have been removed, the apparent area becomes close to the actual area, and the response becomes elastic. In the spirit of our multiscale model, it is the bulk modulus under these various conditions that must then be propagated up to the continuum level, so, accordingly, some method of calculating the low-strain effective bulk modulus becomes imperative.

Large-scale bilayer atomistic molecular dynamics (MD) simulations with length scales on the order of 20 nm using the GROMACS force field (Spoel et al., 1996) have been performed (Lindahl and Edholm, 2000; Marrink and Mark, 2001). However, to obtain computational feasibility, electrostatic interactions were handled with a cut-off. Although the structural properties are not severely altered by an incomplete treatment of the electrostatics, there is strong evidence that the material properties (shear viscosity, bulk viscosity) are quite sensitive to long-ranged interactions. For example, in Feller et al. (1996) and Wheeler et al. (1997), the shear viscosity for water and methanol, respectively, were calculated using an Ewald summation (de Leeuw et al., 1980; Essmann et al., 1995) compared with a cut-off. Although the viscosity calculated with full Ewald treatment of the electrostatics under periodic boundary conditions was in excellent agreement with experiment, the shear viscosity calculated via a cut-off was overestimated by almost an order of magnitude. Reductionist lipid MD models as in Shelley et al. (2001) have also been developed, which can successfully reproduce key equilibrium structural properties, for example the membrane thickness. However, then the ability to reproduce other properties such as the bulk modulus is unknown. Other meso-scale models for lipid membranes have been developed (Goetz and Lipowsky, 1998; Goetz et al., 1999), as well as mesoscale models for the cross-linked actin filament networks found in the human erythrocyte cytoskeleton (Boey et al., 1998; Discher et al., 1998). Recently, Brownian dynamics have also been used (Noguchi and Takasu, 2002).

Meso-scale simulations with dissipative particle dynamics (DPD) (Hoogerbrugge and Koelman, 1992; Koelman and Hoogerbrugge, 1993; Espanol and Warren, 1995; Espanol, 1995, 1996; Groot and Warren, 1997; Marsh et al., 1997) have been used to model a lipid bilayer in Groot and Rabone (2001), where individual lipid molecules were modeled using the characteristic "soft" DPD potential (Espanol and Warren, 1995; Groot and Warren, 1997) using a Flory-Huggins treatment to parameterize different "molecular" interactions. The soft conservative force used in this version of DPD gives the static equilibrium pressure of the system, and it has the feature of having no discontinuity as the interparticle distance approaches zero. In other words, DPD particles can "pass" through one another. This feature is an important part of the mesoscopic interpretation associated with DPD; the "particles" are to be thought of as small dynamically evolving clusters of molecules, i.e., mesoscopic hydrodynamic entities. By modeling a bilayer with an approach such as DPD as in Groot and Rabone (2001), larger systems can, in principle, be examined with significant time-step and length-scale increases. However, in their particular meso-scale membrane model, there are two potential concerns. First, the abstraction of using soft potentials to model actual lipid molecules is problematic: What does it mean physically that one lipid molecule can pass through another due to the soft conservative force? Second, and more importantly, the soft interaction, along with the viscous dissipative forces used in DPD, means that the response of this bilayer model to uniform changes in area is generally not elastic, but viscous. The model therefore fails to reproduce one of the most important physical features of lipid bilayers, i.e., the existence of an elastic bulk modulus.

Lipid bilayers do in fact possess the property of a solid-like elastic bulk modulus along with a fluid-like viscous shear viscosity, and they are defined by a state of zero shear modulus (Evans and Needham, 1987; Hallet et al., 1993). This interesting behavior results from the fact that the lipids diffuse within the plane of the bilayer, yet they respond elastically under uniform area dilations. The value of the bulk modulus is required to resolve the continuum-level equations of motion. However, the shear viscosity will only be needed if shear components are apparent.

To build on the multiscale simulation method using the strategy described in Ayton et al., (2001a,b, 2002b), we have concluded that it is crucial for the intermediate meso-scale representation of the bilayer to be related to information obtained from atomistic MD simulations and then propagated to longer spatial and temporal domains via the relevant material properties. In the case of a lipid bilayer, the key information obtained from atomistic-level MD is the "high-strain" bulk modulus, as noted earlier. Previous mesoscopic models fail to capture this essential quantity by virtue of either approximated electrostatics (Lindahl and Edholm, 2000; Marrink and Mark, 2001) or soft potentials (Groot and Rabone, 2001). We propose in this paper an alternative coarse-grained mesoscopic model for biological assemblies such as lipid bilayers. Green-Kubo Theory formally relates transport coefficients to time-averaged correlations of microscopic quantities. In the case of the shear viscosity, for example, correlations in the P_xy components of the pressure tensor are required. Thus the effect of detailed electrostatic and molecular interactions are averaged in time to result in a material property. Herein lies the basis of our proposed model: rather than attempt to parameterize the meso-scale model using simplified molecular models, we instead use material properties as the fundamental "interaction." In doing so, we correspondingly imply a time-averaged picture. As the description of the system moves out in length scale, from atomistic to mesoscopic and beyond, a corresponding integration and averaging of microscopic interactions into material properties is implicitly performed. In the case of linking microscopic to continuum-level representations, this averaging is formally expressed via the constitutive relations at the continuum level (Ayton et al., 2002a). As mesoscopic properties become important, the procedure becomes more complicated.

To model a bilayer at the mesoscopic level, the bulk modulus is a core property around which the model can be built. That is, the fundamental coarse-grained unit in such an approach is not designed to represent a reductionist model of a lipid, but rather to model the response of a small volume of bilayer to plane stress. Thus the fundamental "particle" is significantly abstracted from an atomic representation to instead represent, essentially, a small region of bilayer whose properties are averaged over a microscopic time.

For this paper, we correspondingly propose a mesoscopic model that captures the solid-like elastic bulk behavior of a lipid membrane as a baseline property. To accomplish this goal, DPD is recast in an elastic membrane formulation. A key feature of the model is that it will obtain key parameterizations from a detailed atomistic description, where the parameterizations for a specific membrane system include not only the bulk modulus, but the membrane thickness, area, and density. With this atomistically determined information in hand, a DPD particle can be constructed that models how a small area of membrane on length scales of the collective microscopic model responds to various in-plane and out-of-plane deformations. We emphasize that this meso-scale membrane model, which we call the elastic membrane DPD (EM-DPD) model, will eventually become a key component of a larger multiscale simulation methodology that ranges from atomistic to continuum levels (Ayton et al., 2001a,b).

The organization of this paper is as follows: In the next section, details of the meso-scale EM-DPD model will be explained. The following section then gives specifics of the algorithm and implementation of the EM-DPD simulation, and, in the following section, the method is used to explore the effects of meso-scale thermal perturbations on lipid bilayers. The results are presented and conclusions are given in the last two sections.

	A MESO-SCALE MODEL FOR A LIPID BILAYER

TOP ABSTRACT INTRODUCTION A MESO-SCALE MODEL FOR... MESO-SCALE SIMULATION OF LIPID... THERMAL OSCILLATIONS RESULTS CONCLUSIONS APPENDIX REFERENCES

Our meso-scale model of a bilayer can be constructed by considering a small area of the bilayer and its response to plane stress. The constitutive relation relating stress to strain in a membrane is given by (Hallet et al., 1993)

(1)

where  = 1/2(P_xx + P_yy), A = A  A₀, and A₀ is the initial area. It can also be expressed in terms of an energy as

(2)

where h is the thickness of the membrane. Eq. 1 can be found from Eq. 2 by evaluating dE/dA, noting the thermodynamic relationship dE = PdV, where dV = hdA, and then using Newton's First Law, P = .

With the initial area density defined as ₀ = N/A₀, Eq. 2 can be written as

(3)

To construct the EM-DPD model, the area A, which is, in general, time dependent in a simulation, is discretized into a number of small elements. These discretized elements will eventually be represented by EM-DPD-like particles. To first order, A/A₀ = 2, where the strain is given in terms of the components 2 = _x + _y.

The energy of A can, in turn, be written as a sum of N small discretized elements that interact via a pair-wise additive interaction. The energy of the ith element is E_i = _ji,rijrcut E_ij where E_ij = (2_ij)²,  is a constant to be determined, and (2_ij)² is related to the local strain between two elements i and j as a function of the interparticle distance given by r_ij = |r_i  r_j|. The constant  can be found by writing the total energy in terms of the average square of local strains as

(4)

where N_cut is the average number of j particles included within the cut-off distance r_cut over N. For small deformations, (_ij)² ~ ², and large N, equating Eq. 3 and Eq. 4 gives

(5)

To this point, the geometry of A₀ is arbitrary but in the case of a circular area, A₀ = r. To evaluate the force on an element i at the center of A, F_i, we must evaluate F_i = E_i. Consider two points in A with relative separation given by r_ij, where i is originally at the origin of A, and r_ij  r_cut. To first order, local strain in terms of the original distance r can be written as

(6)

In this way, the force on i due to j is then given by,

(7)

In the next sections, the specific details of the implementation of this model in an EM-DPD simulation will be presented.

	MESO-SCALE SIMULATION OF LIPID BILAYERS

TOP ABSTRACT INTRODUCTION A MESO-SCALE MODEL FOR... MESO-SCALE SIMULATION OF LIPID... THERMAL OSCILLATIONS RESULTS CONCLUSIONS APPENDIX REFERENCES

An EM-DPD membrane simulation was constructed to model a DMPC membrane under periodic boundary conditions in the xy plane. The present EM-DPD membrane formulation was designed to model the bulk material properties as calculated from atomistic-level MD, where it was observed (Ayton et al., 2002a) that the bulk response to uniform area dilations was elastic, in accord with experimental observations. The bulk modulus was determined using NEMD (Evans and Morriss, 1990) by relating the time derivative of the stress response to an artificially imposed strain rate,  = /2, in the limit that the strain rate was small. Details of the NEMD calculation of the bulk modulus for DMPC are explained in detail in Ayton et al. (2002a).

The aim of the present mesoscopic model is to bridge the microscopically determined material properties to the continuum level. In that regard, the length and time scales for the EM-DPD model are chosen to model material properties and not molecular properties, as stressed earlier. This is in contrast to the DPD model membrane in Groot and Rabone (2001) where DPD was used to model actual coarse-grained lipid molecules.

EM-DPD simulation of lipid bilayers

With EM-DPD, the dynamics of small mesoscopic clusters of molecules are modeled by three pair-wise forces: a conservative force F, a dissipative force, F, and a random force F. Details of the statistical mechanics behind these interactions can be found in Espanol and Warren (1995), Espanol (1996), Groot and Warren (1997), and Marsh et al. (1997), and, here, we only briefly describe them in the context of our lipid-bilayer model. The choice of the conservative force for a lipid bilayer is critical, because it is this interaction that will determine the membrane's elastic bulk properties. Our choice for the conservative interaction is based on the results in Smondyrev and Berkowitz (2001) and Ayton et al. (2002a), where detailed atomistic MD and NEMD simulations were performed on a DMPC bilayer. The bulk expansion modulus, , density, area per lipid, and membrane thickness of DMPC were calculated for an atomistic MD system size in which no long-wavelength bending modes could develop. The parameterization of the current EM-DPD model obtains directly from the results in Ayton et al. (2002a), and thus the pairwise additive conservative force is as derived in Eq. 7 and expressed as

(8)

where  is as defined in Eq. 5. The thickness of the DMPC membrane was that determined in Smondyrev and Berkowitz (2001) with h = 3.4 nm, and the prestressed state modulus was calculated in Ayton et al. (2002a) as 32.7 amu/nm ps². The conservative force as defined here is similar in spirit to bonded forces used in polymer networks as in Groot and Warren (1997), Groot and Madden (1998), and Groot et al. (1999), except that the current formulation bonds EM-DPD particles in a two-dimensional (2D) network rather than a linear spring.

The original configuration for our EM-DPD model was obtained from an isotropic 2D EM-DPD fluid. Also, because the original equilibrium atomistic MD simulation conditions were under zero-stress conditions, the traditional linear conservative force as defined in Espanol and Warren (1995) and Groot and Warren (1997) was not included. The present EM-DPD formulation is designed to model deviations from a zero-stress state. In the case that a solvent, or another membrane, was to be explicitly included in the EM-DPD simulation, this repulsive interaction would have to be specified.

The dissipative and random forces are as defined in Groot and Rabone (2001), but they are now expressed explicit in units of amu nm/ps² as

(9)

(10)

In these expressions,  is the usual strength parameter for the random DPD force (Groot and Warren, 1997), and has the value of  = 4 amu nm/ps^3/2, k_B is Boltzmann's constant, T = 308 K, v_ij = v_i  v_j, where the velocity of particle i is defined in terms of its momenta p_i, p_i = m_iv_i, and the mass of the EM-DPD particle is m_i = _mN_DPDA. The time-step used was set at t = 0.04 ps, and  is a random variable assigned for each pair of interacting particles. Details of the requirements on  and  are described in Groot and Warren (1997).

The mass density _m for DMPC was that found from an equilibrium MD simulation (Smondyrev and Berkowitz, 2001), _m = 595.8 amu/nm³, and N_DPD = 6920. Because this is a pure system, dimensionless units as in Espanol and Warren (1995) and Groot and Warren (1997) are possible. However, to focus on the relevant time and length scales, retaining the fundamental units of mass, length, and time makes comparison with microscopically determined quantities and experiment easier. The weighting function in Eqs. 9 and 10 is given as f(r_ij) = (1  r_ij/h) for r_ij < h, and is zero otherwise.

In this model, the dissipative and random forces model the heat-dissipating viscous fluid-like properties of the membrane, and the viscous interaction of the bilayer with the solvent. Because the parameters used in the EM-DPD model arose from detailed atomistic-level MD simulations with full hydration and long range electrostatics, the effect of the solvent has been collapsed down onto the observed material properties of the membrane itself. The calculation of the bulk modulus at the atomistic level includes the surrounding solvent, and thus incorporates the effects of lipid-water interactions. The dissipative and random interactions now model only the heat transfer between membrane and solvent.

Algorithm

The EM-DPD simulation was constructed via the following algorithm, starting from our atomistic-level MD as done in Ayton et al. (2002a).

STEP 1  Determine the expansion bulk modulus for a DMPC bilayer via NEMD from an atomistic-level MD simulation. This modulus describes the stress response to uniform area dilations upon expansion from an initial state. The system sizes used in the microscopic level calculation are on the order of N ~ 15,000 atoms, with length scales on the order of 3-5 nm. For these small system sizes, thermal fluctuations are dampened via the periodic boundaries. The zero-stress conditions used in the atomistic-level MD correspond to a larger length-scale state of high surface tension. To achieve a "flat" membrane under real conditions, a large surface tension must be applied to "pull out" thermal oscillations.

STEP 2  Construct the meso-scale EM-DPD particle. From NEMD, the "undulation free" modulus originating from a microscopic zero-stress state confirms elastic behavior,  = A/A₀, where A₀ is the initial zero-stress area. Thus, deviations from this initial area will result in an increase in energy given by Eq. 2. As the atomistic-level response for a specific length-scale has been determined, a reasonable choice for the fundamental EM-DPD particle is the atomistic simulation itself. In contrast to Groot and Rabone (2001), the EM-DPD particles here do not represent lipids; rather they are designed to model the bulk response of a meso-scale region of membrane due to area dilations.

STEP 3  Construct the meso-scale EM-DPD membrane. Constructing the EM-DPD model of the bilayer requires both the specification of the conservative force and the geometry of the membrane. A lipid bilayer, due to its elastic bulk response, differs from other systems that have been modeled with DPD. To model the solid-like elastic response, an algorithm was developed that results in an isotropic elastic membrane. First, a 2D EM-DPD fluid in the xy plane at * = N/A = 5 was constructed. In this case the fundamental unit of length, _DPD, was chosen to be h, the membrane thickness. A brief equilibration in two dimensions was performed using parameters as prescribed in Groot and Warren (1997) to generate a 2D EM-DPD fluid. At the end of this equilibration, a replica of the system was created and located a distance h' in the z direction. The two EM-DPD systems are designed to model each leaflet of the bilayer. About each EM-DPD particle, a cut-off radius r_cut = h was used to "tag" all other particles within that radius, including particles in the lower bilayer. With * = 10, r_cut = h = 3.4 nm, and N_DPD = 6920, each EM-DPD particle was linked with, on average, 16 ± 8 other EM-DPD particles within the cutoff. The initial length of a bond, r, between a pair of EM-DPD particles is as specified from the distance determined from the initial configuration. In this way, the isotropic fluid correlations between particles remains intact. Local fluctuations in density will result in some regions of the EM-DPD membrane being more rigid, whereas some regions are less strongly networked. For example, lower initial densities resulted in regions of the membrane where EM-DPD particles were not sufficiently bonded to others, resulting in "holes." The EM-DPD density is not related to the actual membrane density, but is more akin to the numerical resolution of the model. Higher densities resulting in densely cross-linked networks gave similar final results. The chosen density is a trade off between computational speed and membrane stability. For the chosen final density of * = 10, for one monolayer, each particle is bonded with 6-8 particles, near the coordination number of a 2D particle at high density.

This EM-DPD model does not attempt to capture any viscous fluid material behavior. As long as shear stresses are not present, or diffusion within the membrane is not required, this model should reasonably represent a membrane at meso-scale dimensions. Also, as stated earlier, no explicit solvent is incorporated. The EM-DPD dissipative force described in detail in Groot and Warren (1997) and Groot and Rabone (2001) acts in this situation to model the viscous fluid-membrane interaction. The exclusion of the explicit solvent may affect the time dependence of membrane oscillations, but the steady-state properties should not be affected as long as the dissipative force is included. As such, the viscosity of the surrounding fluid is not required, and the magnitude of the dissipative force is chosen to satisfy conservation of kinetic energy.

	THERMAL OSCILLATIONS

TOP ABSTRACT INTRODUCTION A MESO-SCALE MODEL FOR... MESO-SCALE SIMULATION OF LIPID... THERMAL OSCILLATIONS RESULTS CONCLUSIONS APPENDIX REFERENCES

The EM-DPD model as previously described uses small cross-linked domains that are parameterized to respond to area dilations in accord with a detailed atomistic model. The resulting membrane model is free to buckle, and, for undulations out of the plane, the resulting magnitude of the oscillations will arise from frustrations within the random "material bond" network. So, to validate that the model has the correct flaccid behavior, the zero surface tension oscillation spectrum can be used to not only determine the wave-vector dependence of the resulting thermal oscillations, but also the bending modulus, k_c, itself. In contrast to the original atomistic MD model, the EM-DPD model is not used to calculate explicitly the bending modulus, rather the resulting bending modulus arising from the bulk modulus parameterization is used as a qualitative diagnostic to examine the behavior of the resulting thermal undulations. A detailed discussion of the treatment of membrane undulations can be found in Sackmann (1994), Lindahl and Edholm (2000), and Marrink and Mark (2001), and here we will only summarize some important points.

A membrane with area A in the xy plane under zero surface tension will exhibit oscillations in the z direction with a displacement at r = r_x î + r_y given by u(r). This can be expressed as an expansion over Fourier modes u(q) as

(11)

where, in the situation that A = L², the wave vector is q = (2/L)(n_x î + n_y ), and n_x, n_y are integers. The total undulation energy over A is given by E_und = 1/2k_c  |²u(r)|² dr, where k_c is the bending modulus. This expression holds in the case that q h¹ and simplifies in Fourier space to

(12)

On average, in the case of zero surface tension,  = 0, equipartition gives (Lindahl and Edholm, 2000, Marrink and Mark, 2001)

(13)

where k_B is Boltzmann's constant. This expression predicts the undulation spectrum with a 1/q⁴ wave-vector dependence that is inversely related to the bending modulus k_c. Typical values for the bending modulus for lipid bilayers are 24-30 amu (nm/ps)² (or 4-5 × 10²⁰ J) (Lindahl and Edholm, 2000; Marrink and Mark, 2001). We will report the bending modulus in amu (nm/ps)² to reflect the fundamental units of mass, length, and time that govern the meso-scale simulation.

	RESULTS

TOP ABSTRACT INTRODUCTION A MESO-SCALE MODEL FOR... MESO-SCALE SIMULATION OF LIPID... THERMAL OSCILLATIONS RESULTS CONCLUSIONS APPENDIX REFERENCES

The results presented here will be divided into two sections. First, the high surface tension regime will be examined, and then the zero surface tension regime will be studied to quantify the oscillation behavior of the membrane.

High surface tension regime

We begin by examining the initial EM-DPD membrane at * = 10, originally constructed in the xy plane at a state of zero stress and energy (Eq. 2). The initial dimensions of the membrane in the x and y directions were 90 nm, and the separation between EM-DPD membrane layers was set at h' = 3. Because no explicit solvent is included, the length of the z direction of the simulation cell (L_z) was chosen such that the membrane was free to undulate without periodic boundary effects. The simulations were performed under constant NVT conditions. However, given that the density of the membrane is unaffected by altering L_z (beyond periodic effects), a better description of the simulation state parameters is constant NAT, where A is the effective area of the membrane. The actual area includes the thermal undulations.

An initial equilibration run of 40 ns (1 × 10⁶ time steps with a time step of t = 0.04 ps) was performed, followed by a subsequent production run of the same duration. As the membrane was allowed to thermally oscillate in all directions, the final equilibrium stress was found to be  = 4.23 ± 0.2 amu/nm ps² with a corresponding surface tension of 12.19 ± 0.01 amu/ps². The inclusion of thermal modes at wavelengths beyond the original atomistic-level simulation thus resulted in a substantial surface tension. A snapshot of the simulation is shown in Fig. 1, where inspection reveals that the planar profile of the membrane is almost intact. However, small oscillations clearly persist even under this high surface tension. This suggests that the zero-stress microscopic state corresponds to a meso-scale membrane under significant surface tension, but that there is still sufficient thermal energy for small thermal undulations to persist. The persistence of such modes even at high surface tension is in agreement with the experimental observations in Evans and Rawicz (1990, 1997).

View larger version (94K): [in this window] [in a new window]   FIGURE 1   A snapshot of the initial EM-DPD membrane under  = 12.19 ± 0.01 amu/ps2. Small thermal undulations still persist at this high surface tension.

In Fig. 2, the stress as a function of various EM-DPD simulation areas is shown. The membrane as previously described was subjected to either an expansion or contraction in area to different states. Upon either a dilation or contraction, the membrane was allowed to equilibrate under constant NAT to the new geometry for 40 ns, and then production runs of similar lengths were performed. The total stress was found from the negative value of (P_xx + P_yy)/2 for that specific area. The solid square corresponds to the initial flat membrane that was allowed to thermally undulate, and now lies near the onset of the linear stress-versus-area regime. Below this point, the stress-versus-area curve is no longer linear as thermal buckles begin to develop. At absolute areas lower than 6500 nm², thermal undulations dominate the behavior (Evans and Rawicz, 1990, 1997).

View larger version (15K): [in this window] [in a new window]   FIGURE 2   Absolute stress versus area for the EM-DPD fluid. The solid circle corresponds to results using a constant NT algorithm with  = 0. The solid square corresponds to the initial perfectly planar starting state of the membrane.

A more direct examination is found by using the zero surface tension area A to evaluate ln() versus (A  A)/A. The determination of the  = 0 area will be discussed in the next section. The slope of this plot in the low tension regime is related to the elastic bending modulus k_c, and, as tension increases, a crossover regime is found where the slope approaches the elastic expansion modulus K_A (Evans and Rawicz, 1990; Rawicz et al., 2000) where K_A ~ h. From Fig. 3, the EM-DPD membrane exhibits both the thermally dominated low strain regime at (A  A)/A < 0.05, and a linear elastic regime above (A  A)/A > 0.1. With enough simulation points around  = 0, the bending modulus k_c could, in principle, be calculated from this slope. Alternatively, k_c can also be calculated via the undulation spectrum and Eq. 13, and it is this latter method that is used in the present work.

View larger version (13K): [in this window] [in a new window]   FIGURE 3   ln() versus strain (A  A)/A starting from  = 0. The low strain regime, (A  A)/A < 0.05, is governed by thermal undulations, whereas the high strain regime (A  A)/A > 0.1 is dominated by linear elastic behavior.

To calculate the high-strain bulk modulus, it is better to examine the stress  versus strain (A  A₀)/A₀ for an initial area, A₀, that corresponds to a prestressed state. An obvious choice for the prestressed state is the initial starting point (the solid square in Fig. 2). From this prestressed starting point, thermal undulations have been severely dampened but not erased. So, it is predicted that the corresponding stress versus strain behavior will include small perturbations due to the presence of soft thermal modes. Keeping in mind that the initial parameterization for the EM-DPD unit had  = 32.7 amu/nm ps², as found from Ayton et al. (2002a), the stress-versus-strain behavior about the initial starting area (shown as the solid symbols in Fig. 4) yields a modulus slightly less than the input value, with _DPD = 31.8 ± 0.3 amu/nm ps².

View larger version (18K): [in this window] [in a new window]   FIGURE 4   Expansion and contraction stress and strain for the EM-DPD fluid around the initial state point. Absolute values of  and A/A0 are shown for the compressed states.

The above result explains the consistently overestimated expansion moduli as calculated from small atomistic-level simulation using NEMD. The small system sizes used in Ayton et al. (2002a,b) by definition cannot exhibit thermal undulations, and the calculated modulus reflects that ideal planar geometry. When that information is bridged to the meso-scale, thermal effects introduce soft undulation modes, and the resulting modulus is less. In this way, the EM-DPD simulation acts as a thermal-reservoir that includes thermal-undulation perturbations to the original nonperturbed modulus.

In Ayton et al. (2002a) it was also observed that the compression modulus _cont was over a factor of three larger than the corresponding expansion modulus _exp. This was in apparent disagreement with experiment (Koenig et al., 1997), where the two moduli were found to be nearly equal. The explanation for the discrepancy at that time was attributed to the fact that, in real experiments, it is essentially impossible to compress a membrane without exciting the soft bending modes, as opposed to directly compressing the membrane within in a perfect plane. However, with the atomistic MD-scale system sizes used in Ayton et al. (2002a), the compression modulus reflects the latter scenario. The same situation can be examined using EM-DPD. Upon compression from the initial state, the corresponding meso-scale modulus (as found from the slope of stress versus strain) is, in fact, slightly less than the expansion modulus, despite the EM-DPD parameterization containing the modulus determined from atomistic MD. At the meso-scale, the explanation is clear: under compression, the EM-DPD particles will not compress in the plane according to Eq. 2 (with  = _cont upon contractions and  = _exp upon expansions), rather they will respond by buckling out of plane to create an undulation. This deformation is not as energetically costly as a direct compression.

A series of simulations where Eq. 2 was modified to include both  = _exp (for local expansions) and  = _cont (for local contractions) was used to test this hypothesis, and the resulting high-strain expansion EM-DPD modulus was indeed larger than the original atomistically MD-determined input value. This result makes sense: the membrane will never elect to directly compress within the plane, even locally. Even under prestressed states, compressive local fluctuations persist due to the persistence of high-strain thermal undulations. In these regions, the response of the membrane to the compressive stress will take the form of undulation, protrusion, and peristaltic bending modes (Marrink and Mark, 2001). Thus the inclusion of the compression modulus in Eq. 2 does not correctly model the bilayer's compressive behavior, and thus will overestimate the stiffness of the membrane.

It is noteworthy that the original experimental set-up used to calculate the compression modulus involved multi-lamellar systems (Koenig et al., 1997), and not a single bilayer. An obvious question to be raised involves the effects of small thermal undulations in multi-lamellar systems. The existence of small out-of-plane thermal undulations, even in a multi-lamellar case, might alter the membrane's compressive behavior. Upon compression, a bilayer within the multi-lamellar system could respond not only by compression within the plane of the membrane, but also by buckling out of plane. One can even imagine a scenario where collective bucklings of adjacent bilayers might occur.

As a first probe into the effects of membrane buckling in multi-lamellar systems, we have performed EM-DPD simulations of a multi-lamellar system including solvent interactions. The inter-bilayer spacing was set to match that found in the experimental system at ~50 Å. Three EM-DPD membranes (with the original parameterization as used in the Algorithm) under periodic boundaries with a DPD solvent as described in Groot and Warren (1997) were used to construct a meso-scale multi-lamellar system. An identical simulation protocol as described in High Surface Tension Regime was used. We present these results here only to indicate the possibility of undulatory modes in a multi-lamellar system. The exact EM-DPD-solvent interaction has not been tuned to explicitly represent the solvent-lipid interaction; rather it acts only to propagate undulations from one bilayer to the next via DPD particle interactions.

In Fig. 5, we show the stress-strain compression plot, where the slope gives the compression modulus of 27.9 amu/(nm ps²) (cf. 32.7 amu/(nm ps²) is the NEMD result), demonstrating that small thermal undulations persist even in multi-lamellar systems, and that the behavior of the EM-DPD system is considerably different from the completely dampened atomistic-level bilayer simulation. We note that the system used in the microscopic NEMD calculations was also under periodic boundaries and, as such, represents a perfectly flat multi-lamellar system. The slightly smaller modulus obtained from the stress-strain plot is due to thermal undulation modes.

View larger version (15K): [in this window] [in a new window]   FIGURE 5   Compression stress versus strain for a multilamellar EM-DPD membrane.

Small undulatory modes are visible in Fig. 6 where a cut-away snapshop of the multi-lamellar EM-DPD system is shown. Panel a shows the full DPD solvent (white) and EM-DPD membrane (grey), whereas panel b shows another cut-away view of the same system where only the EM-DPD membrane is shown. These results are shown as a qualitative demonstration that undulation modes can persist in multi-lamellar systems. We present these new results to reinforce our viewpoint that the value of the compression modulus obtained from NEMD and multi-lamellar osmotic stress experiments is different due to the drastically different boundary conditions. The small system size used with atomistic NEMD results in an ideal compression modulus that only contains the effects of in-plane membrane compression. Real systems can contain small out-of-plane undulations, even in multi-lamellar systems, and the present EM-DPD multi-lamellar results support this point.

View larger version (88K): [in this window] [in a new window]   FIGURE 6   Snapshots of the EM-DPD multilamellar system where (a) includes the solvent (white) and (b) shows the same system but with only the EM-DPD membrane. Small bilayer undulations are observed.

Zero surface tension regime

The zero surface tension regime can be used to examine the thermal undulation behavior via Eq. 13. A planar membrane under  = 0 conditions is expected to have a 1/q⁴ undulation spectrum behavior for q < q₀, where q₀ is a critical frequency. Above the critical frequency, a 1/q² behavior is expected. Simulation studies by Lindahl and Edholm (2000) and Marrink and Mark (2001) observed this behavior and found q₀ ~ 1 nm¹. The system sizes used in that study were on the order of 20 nm for a typical box length. So, q₀ ~ 1 nm¹ corresponds to an undulation wavelength of ~6 nm. Above this frequency, protrusion and peristaltic modes dominate.

The current EM-DPD model consists of linked particles with a fundamental length scale of 3.4 nm, almost an order of magnitude in linear dimension larger than atomic length scales. Thus at high frequencies, the resolution of the model reaches its coarse-grained limit. For that reason, the undulation spectrum of the EM-DPD simulation is expected to break down above some critical frequency, where the detailed atomistic model spectrum must be resolved. However, at low frequencies, the undulation spectrum of the EM-DPD model may be examined.

The EM-DPD membrane as constructed in the Algorithm section, and exactly as used in the previous high-tension results was now examined under a state of zero surface tension. Again, the EM-DPD bilayer was composed of two leaflets with a separation of h' = 3 nm, and cut-off at r_cut = 3.4 nm. The combination of the degree of cross-linking of the EM-DPD bonds, along with the separation of the EM-DPD leaflets results in a resistance to bending undulations. Thus, calculating the bending modulus from an examination of the undulatory mode spectrum can be used as a test of the parameterization of the model.

To examine the zero surface tension regime, a variation of the EM-DPD algorithm in the spirit of Elliot and Windle (2000) was implemented. Nose-Hoover feedback was used to keep the average surface tension  = h[P_zz  1/2(P_xx + P_yy)] equal to zero on average. Specifically, this was accomplished by augmenting the equations of motion in the x and y directions, the details of which are described in the Appendix. To test the algorithm, results from detailed constant NAT runs ~ = 0 were compared to the results from the constant  = 0 algorithm. The resulting areas from the NT simulations were in excellent agreement with the NAT values. This is shown in Fig. 2, where the solid circle is the resulting stress and area under  = 0 conditions.

Under these dynamics, a state of  = 0 is maintained, and Eq. 13 can be used to calculate the bending modulus of the membrane. For this model, the calculation of the bending modulus is used as a quantitative measure of the low-strain thermal oscillations for the EM-DPD model. Recall that the parameterization of the EM-DPD model contains the membrane thickness, density, and bulk modulus. By examining the oscillation spectrum for the EM-DPD model, not only can the q⁴ regime be examined, but the resulting value of the bending modulus gives an indication whether the current EM-DPD parameterization is too soft or too rigid.

In Fig. 7, the raw undulation spectrum u²(q)A/k_BT for  = 0 is shown. Fitting this curve to 1/q⁴ in the regime where 1/q⁴ undulations dominate gives the bending modulus as k_c = 27.4 ± 1 amu nm ps² (4.6 ± 0.2 × 10²⁰ J), which is in excellent agreement with simulation (k_c = 4 × 10²⁰J) (Lindahl and Edholm, 2000) and experiment (k_c = 5 × 10²⁰ J) (Evans and Rawicz, 1990). A more detailed examination of the high-frequency regime is shown in Fig. 8, where the data is represented logarithmically. A cross-over regime between q⁴ and q² behavior is observed around q ~ 0.2 nm¹. The location of the cross-over point occurs at a lower frequency than is observed in atomic-level simulation (Lindahl and Edholm, 2000; Marrink and Mark, 2001), where the cross-over point occurs at q ~ 1 nm¹. The reason for the shift is due to the size of the EM-DPD particle (3.4 nm) versus typical atomistic length-scales (0.3 nm) as discussed earlier. With meso-scale particles, high-frequency undulations do not exist for length scales below _DPD. Correspondingly, undulations at frequencies just below q ~ 2/_DPD are also perturbed. The resulting effect of increasing the fundamental coarse-grained-size unit is to shift the cross-over between q⁴ and q² to lower frequencies.

View larger version (15K): [in this window] [in a new window]   FIGURE 7   The undulation spectrum for  = 0. The dotted line corresponds to the 1/q4 curve fit yielding kc = 27.4 ± 1 amu nm ps2.

View larger version (10K): [in this window] [in a new window]   FIGURE 8   The undulation spectrum for  = 0. The high-frequency regime exhibits a cross-over regime at q ~ 0.22 nm1. The dotted line designates the region dominated by undulations with a 1/q4 behavior.

The low-frequency undulations, however, are consistent with a DMPC bilayer under zero surface tension. The accessible length and time scales used in the current EM-DPD model make examinations over length scales on the order of 100 nm and time scales on the order of 50 ns feasible. The low-frequency oscillations that result indicate that, under zero surface tension, significant undulations exist. A snapshot of the  = 0 system after 50 ns of equilibration is shown in Fig. 9. In contrast to the original state (Fig. 1), significant undulations now persist. Comparing these two figures clearly demonstrates the effect of what a microscopic zero-stress state corresponds to when longer length scales are examined. The true zero surface-tension state is one where significant thermal undulations are present, and, to resolve these undulations, the minimum system size seems to be around that used in Lindahl and Edholm (2000) and Marrink and Mark (2001). With the present EM-DPD model, system sizes at least an order of magnitude larger than those in the previous references are feasible. Furthermore, by virtue of the larger EM-DPD time step and short-ranged effective interaction, the accessible time scales are easily on the order of 100-1000 ns or even greater.

View larger version (67K): [in this window] [in a new window]   FIGURE 9   A snapshot of the  = 0 system. This image was acquired after over 50 ns of equilibration.

Alternative EM-DPD membrane parameterizations

The current parameterization for the EM-DPD membrane has a reduced density of * = 10, a EM-DPD cut-off of r_cut = h = 3.4 nm, and a membrane leaflet separation of h' = 3 nm, resulting in each EM-DPD particle being linked with, on average, 16 ± 8 other EM-DPD particles within the cutoff. It is possible that other alternative parameterizations could be constructed. Moreover, the particular choice of the cut-off, along with the parameterization set, may have an effect on the structure of the EM-DPD membrane.

The EM-DPD parameterization and implementation, as discussed in detail in the Algorithm, begins with the calculation of the bulk modulus at the microscopic level, originating from a small system under periodic boundary conditions. The small system size at the MD level inhibits any long-wavelength undulations, and, as such, the calculated bulk modulus is representative of an idealized perfectly flat bilayer. If a larger MD cell were chosen, it is possible that thermal buckles may begin to develop and the calculated modulus will differ. The choice of the cut-off of 3.4 nm is representative of the size of the MD cell from which the bulk modulus was obtained; thus the cut-off is restricted to length scales similar to those of the original the MD cell. Using much larger cut-offs would require a new set of MD-level simulations to calculate a new bulk modulus. We constructed two additional parameterizations within the allowable length scales. For each, the separation between EM-DPD leaflets was decreased to h' = 2.5 nm (from h' = 3 nm) and the cut-off was decreased to r_cut = 2.72 nm. However, the density of the EM-DPD points was varied from N_DPD²/A = 10 (resulting in each EM-DPD particle being bonded, on average, to 6 other particles) to N_DPD²/A = 20 (resulting in each EM-DPD particle being bonded, on average, to 12 other particles). In Fig. 10, the oscillation spectrum is shown for the two different parameterizations. The solid circles correspond to N_DPD²/A = 10, whereas the solid squares correspond to the N_DPD²/A = 20 system. The open circles are the oscillation spectrum for the original parameterization.

View larger version (12K): [in this window] [in a new window]   FIGURE 10   The oscillation spectrum for two different parameterizations, where h' = 2.5 nm and rcut = 2.72 nm. The solid circles correspond to NDPD2/A = 10, whereas the solid squares correspond to the NDPD2/A = 20 system. The open circles are the oscillation spectrum for the original parameterization.

Qualitatively, the high-density system (solid squares) exhibits an oscillation spectrum similar to the original parameterization, whereas the low-density system (solid circles) exhibits larger thermal oscillations. These results suggest that it is the combination of the cut-off and density that determines the bending stiffness of the membrane. Interestingly, the original parameterization had each EM-DPD particle being bonded, on average, to 16 other particles, which is similar to the present N_DPD²/A = 20 high-density system.

It should be noted that the bending modulus was not an initial input to the model, rather it was measured in the simulation to "tune" the model. The resultant bending-mode amplitudes are determined, roughly speaking, by the degree of cross-linking within the membrane. A high-density cross-linked membrane will exhibit significant bending resistance, whereas a low-density structure will have greater bending-mode amplitude.

As a rough guide to obtain alternative EM-DPD parameterizations in terms of cut-off and density, the resulting bond network should contain, on average, at least 12 to 16 particles. The chosen values for the cut-off should be similar to the original dimensions of the MD simulation.

Finally, one important point is to discuss the effect of not including the shear viscous behavior associated with real bilayers. The consequence of a bilayer being defined by a state of zero shear modulus, and instead possessing a shear viscosity, is that the membrane cannot support a shear. If a particular deformation resulted in a shear stress, the lipids in a real membrane would diffuse, and the shear would not be supported. For the present model, where we have parameterized the EM-DPD interaction according to the bulk modulus, we can examine only those deformations that result in no shear stress on average. A nonzero value of the off-diagonal component of the pressure tensor P_xy will indicate the presence of a shear stress, and if this is observed, then the present model cannot be used. For the high surface tension regime  ~ 12 amu ps², P_xy = 0.045 ± 0.03 amu/nm ps² as averaged over 40 ns. Likewise, for  = 0 under constant NT dynamics, P_xy = 0.03 ± 0.01 amu/nm ps², indicating that, within simulation error, the membrane is not under shear stress, and that the bulk modulus parameterization is valid.

	CONCLUSIONS