INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS A CONSTITUTIVE RELATION FOR... NONEQUILIBRIUM MOLECULAR... RESULTS CONCLUSIONS REFERENCES

The dynamics and structure of biological membranes can be described on time scales ranging from the microscopic (picoseconds to nanoseconds), to the macroscale (milliseconds and micrometers) (Pastor and Feller, 1996; Feller and Pastor, 1997; Tieleman et al., 1997; Forrest and Sansom, 2000). At the microscopic level, the individual motions of lipid and nonlipid molecules are resolved, resulting in diffusion within the bilayer (Essmann and Berkowitz, 1999) and conformational changes in individual lipid structures. The latter gives rise to microscopic fluctuations in mechanical thermodynamic properties such as the instantaneous area per head group (Smondyrev and Berkowitz, 1999a). The time scales for these motions can be on the order of nanoseconds (Pastor and Feller, 1996). For example, the frequency for trans-gauche isomerization occurs with frequencies of 10-20 ns¹ (Venable et al., 1993; Pastor and Feller, 1996). Rotations (wobbling) also occur on time scales of nanoseconds (Pastor and Feller, 1996), while examining lipid translations over "long" times is still not computationally feasible. Current state-of-the-art simulations are now at the point where even membrane "rafts" (Simons and Ikonen, 1997; Prallea et al., 2000) can be modeled and examined on time scales on the order of nanoseconds. However, these time and length scales are still orders of magnitude below a macroscale description of the system.

Currently, to model a membrane at macroscopic spatial and temporal scales can be accomplished via continuum level models (York et al., 1999; Olbrich et al., 2000) where all molecular detail has been removed and the membrane is represented as a smooth continuum. Standard finite-element or finite-difference methods can be used to model membranes numerically (Koivurova and Pramila, 1997), but they are limited to dealing with small deformations.

To solve the continuum level equations of motion, a constitutive relation is required that relates (in the case of an elastic material) stress to strain. The coefficients of the constitutive relation (the bulk modulus, for example) determine the nature and degree of the material's response to strain. However, the actual value of the particular constitutive coefficient is, in principle, an average over microscopic motions.

Experimentally, the material properties of bilayers (e.g., elasticity, bending modulus) have been calculated using, for example, micropipette pressurization of giant bilayer vesicles (Kwok and Evans, 1981; Schneider et al., 1984; Evans and Needham, 1987; Mui et al., 1993; Olbrich et al., 2000; Rawicz et al., 2000) and with nuclear magnetic resonance and x-ray diffraction (Koenig et al., 1997). A theoretical model for osmotic swelling and lysis was proposed by Hallet et al. (1993), where the membrane was modeled as a thin spherical shell, and the stress was related to changes in the surface area of the membrane.

Our particular interest in membrane material properties, such as the bulk modulus, involves a strategy to overcome the present computer simulation gridlock that is the case in large bioassembly modeling. To reach the macroscale is currently impractical with atomistic simulations (which are limited to the submicrosecond and submicron time and length scale). In two previous papers (Ayton et al., 2001a,b), we presented a new method for bridging time and length scales in membrane systems that involves direct interfacing of continuum and microscopic simulations via a coupled calculation of the membrane's constitutive coefficients. Essentially, macroscopically required constitutive coefficients for a particular membrane are calculated with a microscopic-level simulation. In turn, continuum-level densities are used as initial states for new microscopic-level membrane simulations. The result is a computer simulation technique that has the capability of jumping time and length scales, and effectively passes information across spatial and temporal scales. However, for computational efficiency purposes, highly simplified membrane models were used in the preliminary studies. To fully realize the possibility of micro-to-macro bridging in bilayer simulations, it is necessary to demonstrate that the required constitutive coefficients, in particular the bulk modulus , can be accurately calculated with current atomistic-level simulation. This is the topic of the present paper.

To date, the bulk modulus for dipalmitoylphosphadidylcholine (DPPC) has been calculated using volume fluctuations at equilibrium (Shinoda et al., 1997) using a small system of 32 lipids and 434 water molecules. There have also been efforts to calculate the area expansion modulus K_A for membranes using computer simulation (Feller et al., 1995; Zhang et al., 1995; Feller and Pastor, 1996, 1999) using a much larger system (72 lipids and 1376 waters). The area compressibility modulus (Feller and Pastor, 1999) is defined as

(1)

where  is the surface tension,  is the strain, and T is the temperature. A constant surface-tension equilibrium ensemble was used. In these studies, examinations of the effects of surface tension on membrane structure led to a determination of K_A for a DPPC lipid bilayer by examining the changes in area ratio, A/A, resulting from gradually increasing the surface tension, and are in good agreement with experiment. It was found that, using the fluctuations in membrane area as a means of calculating K_A, overestimated the experimental value by a factor of four. Interestingly, this effect was not observed in Shinoda et al. (1997), and may be linked with the unusually small system size.

However, in the context of linking microscopic and macroscale representations of a lipid bilayer, we are interested in calculating the bulk modulus , rather than the area compressibility modulus K_A, because it is  that appears in the continuum-level stress-strain constitutive relation. In this paper, we present a methodology whereby the bulk modulus for a bilayer can be calculated from a microscopic-level molecular dynamics simulation. The specific application is to dimyristoylphosthatidylcholine (DMPC). The basis of our technique relies on nonequilibrium molecular dynamics (NEMD) (Hoover, 1985; Evans and Holian, 1985; Evans and Morriss, 1990), in particular, a method known as cyclic compression (Hoover et al., 1980b). With cyclic compression NEMD, the bilayer is subjected to artificially induced area changes of a specified form. The bulk modulus is found from monitoring the stress (or pressure) response of the system due to the imposed strain.

This paper will be divided into a number of sections. In the next section, the details of the DMPC bilayer simulation will be discussed. The following section will involve a continuum-level discussion of the relevant constitutive relation for a membrane, and then the relationship of this formalism to the microscopic description will be presented in the section following. The results for the bulk modulus calculation for a DMPC membrane will be presented in the next to last section, and conclusions are then given in the final section.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS A CONSTITUTIVE RELATION FOR... NONEQUILIBRIUM MOLECULAR... RESULTS CONCLUSIONS REFERENCES

The molecular dynamics simulations were performed using the DL_POLY simulation package version 2.12 developed in Daresbury Laboratory, England (Smith and Forester, 1999). The DMPC bilayer was composed of 64 lipid molecules hydrated by 1312 water molecules, which corresponds to 20.5 waters per lipid molecule. DMPC molecules were modeled using a united atom force field (Smondyrev and Berkowitz, 1999b). The water model used in our simulation was TIP3P (Jorgensen et al., 1983). All bond lengths were constrained using the SHAKE algorithm with a tolerance of 10⁴, allowing the use of an integration time step of 0.002 ps. Electrostatic interactions were calculated via the particle mesh Ewald (Essmann et al., 1995; Sagui and Darden, 1999) method using a tolerance of 10⁴. The real space part of the Ewald sum and van der Waals interactions were cut off at 10 Å. The temperature of the system was kept constant at T = 308 K using the Nose thermostat with a relaxation time of 0.2 ps. Initial structures of a DMPC bilayer in water (Smondyrev and Berkowitz, 1999a) were equilibrated for several nanoseconds.

	A CONSTITUTIVE RELATION FOR A DMPC MEMBRANE

TOP ABSTRACT INTRODUCTION METHODS A CONSTITUTIVE RELATION FOR... NONEQUILIBRIUM MOLECULAR... RESULTS CONCLUSIONS REFERENCES

In this section, the constitutive relation applicable to a bilayer is presented. At the continuum level, the membrane contains no explicit molecular structure, and instead, its material properties are contained within the constitutive relations. The constitutive relations relate membrane deformations to external stresses via area dilation (area changes), surface shear, and bending. In this membrane model, the membrane thickness and changes in thickness are assumed to be negligible. The elastic (conservative) and viscous-liquid (nonconservative) constitutive relations appropriate for thin materials such as bilayers has been discussed previously (Evans and Needham, 1987; Hallet et al., 1993; Needham and Nunn, 1990). Here we will only present the relevant results.

For a homogeneous membrane in the xy plane with thickness in the z direction, L_z, small compared to the other directions L_x and L_y (i.e., L_z/ ~ 0 where A = L_xL_y), the constitutive relation for stress and strain is expressed as

(2)

where the stress, , is the average of the plane stress in the x and y directions,  = (_x + _y)/2, and likewise, the strain is  = _x = _y. Lipid bilayers in the liquid crystal phase exhibit nonzero diffusion within the plane of the membrane, and are thus defined by a state of zero shear modulus. The bulk modulus is given by  and relates uniform area changes to stress.

In this formulation, the sign of the strain (i.e., whether  > 0, corresponding to an expansion, or  < 0 and the system is contracted) is not specified. In principle, an expansion modulus _exp can be defined for strains with  > 0, and likewise a contraction modulus _cont can be applied to situations where  < 0. Depending on the material and the initial conditions (for example, zero stress, or a state with positive surface tension), these two moduli may or may not be equal. We note that, in almost all experimental situations (Kwok and Evans, 1981; Schneider et al., 1984; Evans and Needham, 1987; Mui et al., 1993; Koenig et al., 1997; Olbrich et al., 2000; Rawicz et al., 2000; Needham and Nunn, 1990), either the expansion or contraction modulus is determined, and, thus, to compare with experiments, it is necessary to separate the bulk modulus in this way.

In atomistic level molecular dynamics simulation, the thickness of the membrane is modeled in detail. In bilayer computer simulations, periodic boundary conditions in the x and y directions are used to model a macroscopic membrane. The z direction is also under periodic boundaries resulting in a periodic lamellar structure. However, if the solvent barrier between periodic bilayers is sufficiently large, then the correlated motions of the adjacent layers should be small, and the system will reasonably approximate a single bilayer. Regardless, the inclusion of periodic boundaries in the x and y directions ensures that a membrane constitutive model as proposed in Eq. 2 is valid.

The finite thickness of the membrane at the microscopic level requires that the response of the membrane in the normal direction to dilations within the plane of the membrane is specified. For example, under dilations in the xy plane, the thickness of membrane may slightly contract. The optimal boundary condition is that the membrane thickness can freely respond to area dilations, corresponding to a state of constant zero stress in the z direction. The details of how various boundary conditions are implemented at the microscopic level will be discussed in the next section.

	NONEQUILIBRIUM MOLECULAR DYNAMICS

TOP ABSTRACT INTRODUCTION METHODS A CONSTITUTIVE RELATION FOR... NONEQUILIBRIUM MOLECULAR... RESULTS CONCLUSIONS REFERENCES

A large number of mechanical thermodynamic observables are easily calculated by using equilibrium computer simulation techniques. For example, the internal energy U, hydrostatic pressure P, and diffusion coefficient D are all most readily evaluated at equilibrium. However, many quantities are more efficiently calculated by computer simulation when the system is out of equilibrium. For example, NEMD has been used successfully to calculate the viscosity of ethanol (Wheeler et al., 1997). With NEMD, an external field, force, or flux is imposed on an equilibrium system by directly altering the equations of motion (Evans and Holian, 1985; Hoover, 1985; Evans and Morriss, 1990). In fact, the present algorithms used to maintain constant temperature and pressure (Hoover et al., 1980a; Hoover, 1985; Allen and Tildesley, 1987; Evans and Morriss, 1990), as well as more membrane-specific algorithms used to maintain constant surface area (Feller and Pastor, 1996, 1999), use the same extended dynamics as NEMD. However, there are a number of subtle aspects concerning the statistical mechanics associated with NEMD. In essence, any viscous heat generated must be removed, and the mechanism by which heat is removed results in changes to the multidimensional phase space in which the system is represented as a single multidimensional point. Because we are interested in the zero strain-rate limit only, these aspects of the statistical mechanics of NEMD will be not be discussed. A complete description is given in detail elsewhere (Evans and Morriss, 1990; Mundy et al., 2000).

For the purpose of calculating material properties, NEMD is an efficient and elegant method where the quantities of interest are evaluated in the limit that the nonequilibrium perturbation goes to zero. The advantage of this approach over equilibrium techniques for the calculation of material properties (e.g., Green-Kubo time-correlation functions) is that NEMD does not usually require either long or large simulations. In particular, to calculate the bulk modulus, a form of NEMD known as cyclic compression will be used, and is presented in the next subsections.

Cyclic compression NEMD

Cyclic compression NEMD has been used to calculate the bulk viscosity of fluids, and has proven successful compared to traditional methods (e.g., Chapman-Enskog (Hoover et al., 1980a,b)). This method can be applied to both elastic and viscous materials, and has the advantage that it does not assume the form of the stress response. The idea is to introduce an artificial strain-rate in the form of a deterministic oscillating area change and then correspondingly monitor the resultant stress response (or time derivative of the stress response). The stress response is found by imposing Newton's First Law where a stress component is defined to be the negative of the pressure tensor component,  = P. The modulus is then calculated in the limit that the imposed strain-rate approaches zero.

With cyclic compression NEMD a strain rate is imposed as specified by either a sine or cosine form for the area oscillation. In the case of a sinusoidal oscillation, the strain rate is given by

(3)

where  is the dimensionless amplitude of the oscillation,  = 2/, and  is the corresponding wavelength. A strain rate with a cosine formulation is similarly written; however, the resultant oscillations in area (A = L_xL_y) are different. With a strain rate given by Eq. 3, the area oscillations will not be symmetric about the initial area A. Rather, for  > 0, the oscillations will be biased to only areas greater than the initial area and will have a (1  cos(t)) form. Likewise, negative amplitudes ( < 0) will result in compressions less than or equal to the initial area. This particular strain-rate formulation is useful in simulations where the expansion or contraction modulus (_exp or _cont, respectively) are to be calculated.

A strain rate with a cosine form,  =  cos(t) will result in area oscillations that are symmetric about the initial area, and is the form that was used earlier (Ayton et al., 2001b). For this work, we will use this formulation only for an initial examination of the behavior of the stress response (that is, to determine whether the system can be modeled by an elastic constitutive relation). To compare with experimental estimates of the expansion and contraction modulus for DMPC, the sinusoid formulation of cyclic compression NEMD is then implemented directly into the equations of motion in a spirit similar to traditional NPT algorithms. Also, in conjunction with the artificial cyclic-compression strain rate, a constant stress state of zero in the z direction (the normal of the membrane) is imposed. The state of zero stress implies that the membrane is free to contract or expand under the cyclic compression area dilations and contractions. As the membrane is expanded in the xy plane, it is allowed to freely contract in the z-direction. Likewise, under compression in the xy plane, the membrane can freely expand in the normal direction.

This formulation of cyclic compression NEMD is implemented directly into the equations of motion of both the particles, and the periodic cell vectors. For the simulations of membranes, the equations of motion of each atom are supplemented with terms that generate the cyclic compression in the xy direction, maintain zero stress in the z direction, and also maintain constant kinetic temperature by allowing heat exchange with the heat bath. For a membrane system composed of both membrane and solvent such that, in total, there are N atoms each of mass m_i with positions r = r₁, r₂, r₃, ... , r_N and conjugate momenta p = p₁, p₂, p₃, ... , p_N in a volume V = L_xL_yL_z, the equations of motion required to implement cyclic compression NEMD are

(4)

(5)

(6)

(7)

Here, the force on particle i is F_i = U, where U is the potential of the system. The equation of motion for , as given by Eq. 7, ensures that the stress in the z direction is given by _z = 0. The value of Q_p must be chosen such that the system responds to changes in stress at a faster rate than as given by .

The term p_i is a thermostatting term where the equation of motion for , the thermostat multiplier, is given by Nose-Hoover feedback (Evans and Morriss, 1990),

(8)

where Q is an arbitrary constant, T is the input thermodynamic temperature, and T_k(t) is the instantaneous kinetic temperature given by

(9)

and d is the number of degrees of freedom. It should be noted that, without the strain-rate NEMD terms in the equations of motion, these dynamics are equivalent to the NP_NAT constant normal stress ensemble at equilibrium (Zhang et al., 1995; Feller and Pastor, 1999). The above equations of motion result in a dilation of the system volume given by

(10)

which, as   0, becomes

(11)

As previously mentioned, cyclic compression NEMD does not assume the form of the constitutive relation. The system is subjected to an artificially imposed strain rate (area oscillation), and, essentially, the corresponding stress response can be of any form: elastic, viscous, viscoelastic, or something altogether different. In the previous section, an elastic constitutive model for a membrane was proposed. For this model to be valid, under NEMD cyclic compression, a sinusoidal strain rate should result in a stress response with a cosine form. Likewise, a strain rate described by a cosine will result in a sinusoidal stress response.

A viscous response (as would be the case in a Newtonian fluid) will have a stress response with the same form as the imposed strain rate. In fact, it is this "phase lag" that is responsible for viscous heating (Evans and Morriss, 1990). The experimental consensus is that biological membranes exhibit an elastic bulk response (Evans and Needham, 1987), with a constitutive model as given by Eq. 2, and, as will be further shown in the results, we find that this behavior is also observed in simulation. With the elastic constitutive model confirmed, there are multiple methods in which to extract either _exp or _cont.

In either case (_exp or _cont) the bulk elastic modulus  is defined as the derivative of stress versus strain evaluated at  = 0, and can be written as

(12)

where the time dependence of  and  is explicitly shown. With NEMD cyclic compression the time dependence of  occurs through the imposed strain-rate , and in the case of a sinusoidal strain rate, the stress will oscillate from  = 0 to  = 2 at a rate of . The factor of two arises in the sinusoidal formulation from the asymmetric oscillation about the initial area. Furthermore, if the material exhibits a linear elastic response about the initial stress (in this case  = 0), then  is simply found from the slope of  versus  about  = 0. Because the strain is time dependent under cyclic compression, a continuous spectrum of stress versus strain ranging from  = 0 to  = 2 is found, and it can be shown that a least-squares fit of the stress response from  = 0 to  = 2 is equivalent to averaging the slope along all strains. In the linear region about  = 0, () = ₀, and  = ₀/2. The coefficient, ₀, can be found from a least squares minimization from  = 0 to  = 2 by defining the residual R,

(13)

and by the minimization R/₀ = 0. The linear coefficient of the stress is given by

(14)

In this way, stress versus strain plots obtained from cyclic compression can be used to calculate the bulk modulus. This particular method is different, although complimentary to that used in our earlier work (Ayton et al., 2001a,b), and it is more in the spirit of vesicle experiments (Evans and Needham, 1987; Hallet et al., 1993). Using stress versus strain can only be used for systems that exhibit an elastic response, and furthermore can only be applied to systems that have a linear response about the initial strain.

Eq. 14 applies to only one compression cycle. To obtain a statistically meaningful result, the slope must be evaluated for several NEMD trajectories. In complex systems such as biological membranes, vastly different conformations over large time scales can occur. To this end, a methodology was adopted to reasonably sample the conformational space. An initial NPT equilibrium trajectory for a DMPC membrane in water was performed over 3 ns, where the box-lengths of the system were allowed to vary independently such that the diagonal components of the pressure tensor were each zero. It should be noted that the choice of the initial starting ensemble is not limited to an NPT ensemble. Equilibrium initial configurations could also be selected from constant surface-area ensembles (Feller et al., 1995; Zhang et al., 1995; Feller and Pastor, 1996, 1999; Venable et al., 2000), and indeed future studies will explore how the nonequilibrium stress response depends on the initial choice of the equilibrium ensemble. Presumably, with larger systems, the dependence on the choice of initial ensemble will decrease.

Configurations at regular intervals from 1.0 to 3.0 ns were selected as starting configurations for NEMD cyclic compressions. This is a standard method for initializing NEMD trajectories where the initial ensemble is of importance (Ayton and Evans, 1999). Six to ten configurations from this equilibrium trajectory were then used.

By averaging over multiple initial configurations occurring over a reasonably long simulation time (in this case the trajectories were selected over a range of 2 ns), the conformational space of the system can be better sampled. Ideally, the configurational sampling should occur over as long an equilibrium trajectory as possible. Each starting configuration is then used for a series of NEMD runs with varying strain rates, where the strain rate is systematically decreased by reducing . The frequency of the cyclic compression, , in the stress-strain form, can be thought of as the quantity that determines the rate at which the system is expanded or contracted. By averaging over multiple initial equilibrium configurations, the bulk modulus for the system is now given by

(15)

where ··· means an average over multiple configurations selected from an equilibrium trajectory.

For the NEMD simulations, the selection of amplitudes and frequencies used for the modulus extrapolation must be done carefully. For this work, we chose strains not greater than those expected to cause membrane lysis (Needham and Nunn, 1990), and frequencies on the order of the presently accepted times for reasonable conformational sampling (on the order of 1 ns). The combination of reasonably long cyclic compressions, along with the six different starting trajectories, will give a reasonable sampling of the conformational space. The NEMD cyclic compression runs are not designed to explore large conformational changes; this is accomplished by selecting different starting configurations along the equilibrium NPT trajectory. In a sense, the NEMD cyclic compression gives the modulus representative of the starting configuration, with configurations "similar" in a conformational sense. By averaging the modulus over various disparate equilibrium starting configurations, an average over different conformational spaces can be found.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS A CONSTITUTIVE RELATION FOR... NONEQUILIBRIUM MOLECULAR... RESULTS CONCLUSIONS REFERENCES

The results will be divided here into first a qualitative description of the nature of the nonequilibrium stress response for a DMPC bilayer and then a detailed examination of the bulk modulus for a pure DMPC bilayer. The first section will be used to demonstrate how NEMD cyclic compression can be used to determine the exact nature of the stress response, i.e., whether the response is purely elastic, viscous, or something else. With the nature of the response predetermined, and the appropriate constitutive relation established, techniques as described in the Cyclic Compression NEMD subsection will be used to calculate both the expansion and contraction bulk modulus.

Cyclic compression at finite frequency

As stated earlier, cyclic compression does not assume the form of the constitutive relation, and thus can be used to determine whether a material responds to external strains with an elastic, viscous, or viscoelastic response. To verify whether a proposed constitutive model is applicable, the stress response, without assuming an elastic or viscous constitutive relation, must first be verified. In the finite frequency regime, the response of the system is determined by complex conformational changes of the lipids due to the imposed strain rate. At very high strain-rate frequencies, the system will be unable to respond to the strain, and a highly elastic response is predicted. In the low-frequency limit, the lipids are able to respond to the cyclic compressions.

In the finite frequency case, the stress response to an imposed cosine strain rate is given in Fig. 1. The cosine, rather than sine form, of the strain rate results in a symmetric oscillation about the initial area, and this is used as a preliminary probe to determine the qualitative nature of the stress response. For this system, the frequency of the oscillation is 0.0157 ps¹, corresponding to a cyclic compression wavelength of 0.4 ns. Qualitatively, the stress response is sinusoidal, indicating an elastic response. However, closer inspection shows that the curve is not symmetric about zero stress. The stress response is much smaller in cases when the membrane is expanded from its initial area. In the situation where the system is under contraction, the stress response is greater. In this formulation, where the cyclic compression is symmetric about the initial area, the expansion modulus is defined from the stress response from those times where the system is expanding from zero stress (i.e., 0-500 ps, and 700-900 ps), and likewise the contraction modulus is determined from those regions where the system is contracting from zero stress (500-700 ps, 900-1100 ps).

View larger version (33K): [in this window] [in a new window]   FIGURE 1   The stress response  × 107 Nm2 for  = 0.0157 ps1 plotted versus time t.

From this result, the origin of the asymmetric response is not clear. Clearly, in the nonequilibrium regime, the system exhibits a more intense stress response on compressions than on expansions. Qualitatively, however, this effect can be explained by considering the origin of the stress response in terms of molecular interactions. Upon compression, the repulsive cores of the individual atoms result in a steeply rising pressure, whereas, on expansions, the longer-ranged attractive interactions dominate. Also, the frequency of the imposed strain rate will affect the intensity of the response: for frequencies that are too great, the system, upon either dilation or contraction, cannot adjust.

For systems such as biological membranes, a determination of the appropriate magnitude of strain rates is not trivial. Too large of a strain rate will result in essentially an artificially large elastic response, as the molecules cannot explore the conformational space within the NEMD cyclic compressions. Likewise, too small of a strain rate will result in a poor signal-to-noise ratio. Because the magnitude of the strain rate is given by the product of the dimensionless amplitude  and frequency , the choice of the optimum NEMD conditions can be further complicated. The choice of  is bounded by (in terms of expansions) strains on the order of those that result in lysis (Evans and Needham, 1987; Needham and Nunn, 1990; Hallet et al., 1993). The boundaries for allowed frequencies are less clear, but a conservative choice is that the allowed frequencies should correspond to oscillation wavelengths at least on the order of the time required to obtain reasonable estimates for equilibrium properties. Of course, "reasonable," at this time, is bounded by computational tractability, and is on the order of a nanosecond.

Heat production and elastic response

For this membrane system, we have proposed an elastic form for the constitutive relation as given by Eq. 2. However, the possibility of a viscous response component must still be considered, and, if a significant viscous component is found, then the constitutive relation as proposed in Eq. 2 must be modified. There are two methods that can determine whether a system exhibits a viscous component in the response to the strain rate: (1) Direct calculation of the stress response versus as described in the previous subsection. With cyclic compression, as previously described in the Cyclic Compression NEMD subsection, the xy components of the pressure tensor oscillate in cycle with . A viscous response could result in a small phase shift in the stress response; if no phase shift is observed from the direct response, then the viscous response (if any) is small. A more robust test, however, is to examine the heat (or entropy) production over the course of expansion or compression cycles. A purely elastic response will result in no irreversible heat being produced, whereas a viscous response will. This can be shown by using thermodynamic arguments as follows.

The internal energy change for this system is given by dU = PdV + dq, where dq is the heat added to the system. This expression can also be written in terms of rates as

(16)

and for a system of N atoms with equations of motion given by Eqs. 4 and 5, the instantaneous rate of heat removal is given by q/t = 3dk_BT_k(t)(t), where (t) is the instantaneous value of the thermostat multiplier appearing in Eq. 5, and whose equation of motion is given by Eq. 8. Under isothermal conditions U/t = 0, thus,

(17)

where t_cycle is the time of a compression cycle. Elastic conditions imply that no viscous heat is produced. Thus, for Eq. 17 to hold, (t) = 0 over compression cycles, expressed as

(18)

Thus, an examination of the heat produced over compression cycles can give a direct determination whether the system can be modeled with a purely elastic model. If (t)  0, within simulation error, then a viscous component to the response is indicated.

In all systems studied, whether using a sinusoid or cosine form of cyclic compression, or whether undergoing expansion or contraction, no heat was produced, further supporting the notion that the response of the membrane is elastic. The only possible source of a heat-generating response could arise from the water layers above and below the membrane. However, as will be discussed in more detail later, it will be shown that, under the zero-stress normal boundary condition, and with low enough frequencies, the water will readjust to the changing area without generating viscous heat.

Expansion modulus for the DMPC membrane

In the giant bilayer vesicle experiments of Evans (Kwok and Evans, 1981; Evans and Needham, 1987; Mui et al., 1993; Olbrich et al., 2000; Rawicz et al., 2000), a lipid vesicle is expanded from a zero lateral stress state, and the area expansion modulus K_A (Eq. 1) is calculated from essentially a stress versus tension experiment. In the present cosine formulation of cyclic compression NEMD used in the subsection Cyclic Compression at Finite Frequency, the imposed strain rate induces volume oscillations that are symmetric about the initial area. So that it is possible to directly mirror experimental conditions where only the expansion regime originating from the initial zero-stress state is considered, the sinusoidal form of cyclic compression as given in the Cyclic Compression NEMD subsection and Eq. 3 will be used.

To obtain a benchmark measurement of the bulk modulus for a DMPC bilayer using NEMD, a detailed set of NEMD calculations were performed as described in the Cyclic Compression NEMD subsection. The first step in the modulus calculation is an examination of the stress-strain behavior at finite frequency. In Fig. 2 is shown the stress-strain plot for a frequency of  = 0.0157 ps¹ and amplitude of  = 0.015, where this frequency corresponds to an expansion cycle that occurs over the course of 0.2 ns. The resultant stress versus strain, as previously mentioned, is an average over six separate starting configurations. Qualitatively, the stress increases linearly with strain, and a linear least-squares fit to the data gives a slope of (24.4 ± 3) × 10⁷ Nm². The nonequilibrium modulus can be calculated from Eq. 15, giving a value of (12.2 ± 2) × 10⁷ Nm². A lower frequency system with  = 0.0105 ps¹, corresponding to an expansion over of 0.3 ns also gives a linear stress versus strain relation, but with a smaller least-squares slope of 19.2 ± 3 × 10⁷ Nm² with a modulus of 9.6 ± 3 × 10⁷ Nm².

View larger version (15K): [in this window] [in a new window]   FIGURE 2   Expansion stress versus strain for a frequency of  = 0.0157 ps1 and amplitude of  = 0.015.

An important point is as follows. With cyclic compression NEMD, the system is expanded or contracted in a smooth analytic form, in contrast to the equilibrium method used to calculate K_A (Feller and Pastor, 1999). The smooth expansion from an initial strain of zero to a maximum strain of 2 implies that the rate of stress response along the entire expansion cycle, rather than just the initial and final states, is considered. In other words, with NEMD, one obtains a complete spectrum from  = 0 to  = 2 via the strain rate , rather than just the initial and final points. This essentially gives information on the stress-strain behavior at all strains less than 2, and, in the case of a linear stress-strain relationship, results in a complete stress-strain spectrum.

The nonequilibrium moduli found from a set of frequency-dependent stress versus strain NEMD runs can be collated to find the zero-frequency estimate for the modulus. The frequency dependence can be thought of as the rate that the system goes from a state of zero strain to 2, thus, as   0, the system is being expanded at a slower and slower rate. At the lowest examined frequency, the system undergoes expansion from  = 0 to  = 2 over the course of 1.0 ns. In Fig. 3 the extrapolated expansion modulus for the DMPC membrane is shown. The size of the symbols are roughly the size of the error bars obtained from the least-squares slope of the corresponding averaged stress versus strain plot. A linear least-squares fit to this data gives the zero-frequency bulk expansion modulus as 5.4 × 10⁷ Nm² with an asymptotic error of 4%.

View larger version (15K): [in this window] [in a new window]   FIGURE 3   The frequency dependence of the bulk modulus. The zero-frequency extrapolated value for the expansion modulus exp is found to be 5.4 ± 0.3 × 107 Nm2 with an asymptotic standard error of 4%.

Experimental estimates for the expansion modulus of lipids are on the order of 10⁷ to 10⁸ Nm² (Evans and Needham, 1987; Needham and Nunn, 1990; Hallet et al., 1993), so the present result is in good agreement with these values. An estimate for the bulk modulus can be found from the area compressibility modulus K_A (Hallet et al., 1993). The experimental value of K_A for a DMPC membrane at 30°C is 0.145 Nm¹ (Lipowsky and Sackmann, 1995). Using an estimated DMPC membrane thickness of 34 Å (Smondyrev and Berkowitz, 1999a), an experimental estimate for the bulk modulus is ~2.1 × 10⁷ Nm². The reasonable agreement between _exp from NEMD and the experimental estimate from K_A under area expansion suggests that, at least under expansion, changes in surface tension and changes in stress under increasing strain are similar.

In the vesicle experiment (Rawicz et al., 2000), it was pointed out that, under experimental conditions, a membrane under zero stress contains subvisible thermal bending undulations. The membrane must either be prestressed, or greater strains must be used to remove these undulations to calculate K_A. However, with the small system size used in NEMD simulations, the thermal bending modes are dampened out by definition.

The difference between changes in surface tension  and stress  due to strain can be seen from the definition of surface tension (Zhang et al., 1995; Feller and Pastor, 1999). If we consider the definition of the surface tension, along with Newton's first law  = P, then we can write

(19)

From the definition of surface tension (Eq. 19), it is clear that nonzero contributions to the integral arise only from those regions where   _z. Thus, K_A, as previously defined, is in contrast to the bulk modulus 2 = /, which does not consider differences in  and _z.

Compression modulus for DMPC

The compression modulus _cont can be calculated using NEMD cyclic compression in a similar manner to that in the previous subsection by using a negative value of . The resulting NEMD dynamics generate area oscillations such that the system is compressed from its initial starting state. The six initial equilibrium configurations originating from a zero-stress NPT trajectory used in the previous expansion modulus calculation were used, except that now the NEMD compresses, rather than expands, the system. The stress in the z direction normal to the plane of the membrane was maintained at zero. Frequencies ranging from 0.03142 ps¹ (corresponding to contractions occurring over 0.1 ns) down to 0.00314 ps¹ (corresponding to contractions occurring over 1.0 ns) were used. The dimensionless amplitude was set at  = 0.015, resulting in a total compressive strain of 0.03.

The stress versus strain for a typical compression cycle is shown in Fig. 4. For this system, the frequency is  = 0.0157 ps¹, corresponding to a compression from a zero-strain state to  = 0.03 over the course of 0.2 ns. For clarity, we report the absolute value of the strain. The negative stress is consistent with a positive pressure increase upon compression, and the stress versus strain relationship over the course of the compression is qualitatively linear. The finite frequency compression modulus for  = 0.0157 ps¹ is found to be 18.61 ± 6 × 10⁷ Nm².

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Compression stress versus strain for a frequency of  = 0.0157 ps1 and amplitude of  = 0.15. The negative stress upon compression corresponds to a positive pressure increase. The absolute value of the strain  is shown.

A similar frequency scan as in the expansion modulus case gives the zero-frequency extrapolation for the compression modulus. In Fig. 5, the compression bulk modulus extrapolation is shown along with the previous expansion modulus extrapolation. For all frequencies considered, the frequency-dependent compression bulk modulus is greater than the expansion bulk modulus. Again, the lowest frequency used,  = 0.00314 ps¹, corresponds to a compression of  = 0 to  = 0.03 over the course of a nanosecond. The zero-frequency compression modulus is found to be 16.5 × 10⁷ Nm² with an asymptotic standard error of 12%.

View larger version (17K): [in this window] [in a new window]   FIGURE 5   Compression (open squares) and expansion (solid squares) frequency extrapolations for the bulk modulus. The value of exp is found to be 5.4 ± 0.3 × 107 Nm2 with and asymptotic standard error of 4%, whereas cont is found to be 16.5 ± 2 × 107 Nm2 with an asymptotic standard error of 12%.

Comparison between expansion and contraction bulk moduli

It is useful to directly compare stress versus strain for the same frequency under expansions and contractions. In Fig. 6, the stress versus strain spectrum for  = 0.00628 ps¹ clearly shows that the stress response is greater upon compressions than expansions. This trend is observed for all frequencies. Furthermore, at least qualitatively, _cont > _exp holds for stresses less than || = 0.03. The origin of the difference between the finite frequency expansion and contraction moduli depends not only on the starting zero-stress configuration, but on the time-dependent nature of the expansion. In the expansion case, the system is being pulled from a zero stress configuration, resulting in a complex rearrangement of both lipids and water as they adjust to the larger area. Keeping in mind that P = , we can use the pressure response as a guide to understand the stress response. As the strain increases, the pressure becomes more and more negative (from the original equilibrium trajectory originated from zero stress). The negative pressure response originates from terms in the molecular pressure tensor arising from Coulombic interactions, intramolecular lipid bending and torsion interactions, and short-range van der Waals attractions. As the system is "pulled apart," these attractive interactions become dominant as the short-ranged repulsive Lennard-Jones interactions become less. In a sense, the lipids adjust to the slowly expanding area not so much by "getting out of the way" of other lipids, but by "filling the empty space."

View larger version (16K): [in this window] [in a new window]   FIGURE 6   Compression (squares) and expansion (triangles) stress versus strain spectrum for  = 0.00628 ps1. The dashed line designates () = 0.

The scenario upon compression, especially from the zero-stress state, is quite different. The zero-stress state, from a molecular perspective, represents a balance between short-ranged molecular core repulsions, and complex inter- and intramolecular attractions. Given that the density of the membrane is already quite high in the LC phase, any compression from the original state will result in pressure increases due to strong short-ranged repulsions. The lipids must respond to the compression by adjusting conformations. Clearly, for fast compression frequencies, the lipids cannot adjust to the new area, and the pressure response will be given by a sharp rise in pressure as the short-ranged repulsions dominate. For slower frequencies, the system can begin to adjust and relax to the compression. At some critical frequency, the system can relax to the compression faster than the imposed strain rate.

For the lowest frequencies ( = 0.00314 ps¹), it is unlikely that the difference in the moduli is a result of the nonequilibrium dynamics and the artificial strain rate. The zero-stress state represents a balance point between short-ranged strong repulsions and longer-ranged attractive interactions. Compressions result in sharper increases in pressure (or corresponding decreases in stress) than do expansions for the same strain. To understand what is happening under cyclic compression, three simulation snapshots under the compressed NEMD state (1), the initial zero-stress state (2), and the expanded NEMD state (3) are shown in Fig. 7. These images only apply to one of the chosen zero-stress starting states, and thus can only give qualitative information; however, some insightful trends can be observed. First, it is clear that the NEMD cyclic compression results in a compressed (1) and expanded (3) state of the initial zero-stress configuration. Upon compression, the thickness of the entire system increases such that the total volume of the system essentially remains unchanged. In fact, during the course of a compression, the thickness of the membrane and solvent combined, given by L_z, will increase on the order of 6%, whereas the volume change of the system is on average less than 0.5%. It is interesting that, qualitatively, the thickness of the membrane remains basically unchanged during both expansion and contractions (from the snapshots, it is hard to draw a more precise conclusion). The change in system thickness is mostly due to the reorganization of water above and below the membrane, and this is consistent with no viscous heat being produced.

View larger version (109K): [in this window] [in a new window]   FIGURE 7   Simulation snapshots of (1) the compressed NEMD state, (2) the initial zero-stress state, and (3) the expanded NEMD state for  = 0.00628 ps1. The images on the left are the bilayer structure with water not shown, and the images on the right are for the same system with the water shown.

A more detailed examination of how the bilayer thickness is affected by NEMD expansion can be found from phosphate density profiles obtained during the course of NEMD expansions and contractions. Density profiles obtained over the expansion can be compared to corresponding profiles obtaining during the compression, thus giving an indication of the changes in the bilayer thickness. In Fig. 8, the two NEMD density profiles are superimposed, and show a small change in the bilayer thickness depending on whether the system is undergoing expansion or contraction. Upon NEMD expansion, the bilayer thickness is slightly reduced to maintain zero stress in the normal direction. Conversely, upon NEMD compression, the bilayer thickness increases slightly in response to the decreased area. In both cases, the NEMD perturbation results in a coupling between the membrane area and thickness such that a state of zero normal stress is maintained.

View larger version (22K): [in this window] [in a new window]   FIGURE 8   Distribution of phosphate atom positions along the bilayer normal during NEMD expansion. The solid line corresponds to the expansion ( = 0, ... , 0.03), and the dashed line is for the compression ( = 0.0, ... , 0.03). The frequency of the expansion is  = 0.01571 ps1.

It is possible that a shear stress, due to the small size of the system combined with lipid packing commensuration effects, could account for the difference between the expansion and compression moduli. If the expansion or compression rate was too large, it is possible for a shear stress to manifest, especially upon contraction. In Fig. 9, the off-diagonal component of the pressure tensor, P_xy is shown for both the NEMD expansion and contraction with a frequency of  = 0.01571 ps¹. In both cases, within the fluctuations, P_xy does not vary significantly with the strain, whether  > 0 or  < 0. This result suggests that the lipids have time to adjust to the altered unit cell geometry, and that the magnitude of the imposed strain is not so large as to impose commensuration effects. Furthermore, because we are interested in the extrapolated zero-frequency value of , the absence of a system-size shear stress due to commensuration effects indicates that the NEMD strain is not too large.

View larger version (22K): [in this window] [in a new window]   FIGURE 9   The off-diagonal component of the pressure tensor, Pxy as a function of  for contractions (solid symbols) and expansions (hatched symbols). During the course of both the expansion and contraction, no significant shear stress is observed. The frequency of the expansion is  = 0.01571 ps