INTRODUCTION

TOP ABSTRACT INTRODUCTION A CONTINUUM-LEVEL MODEL FOR... THE ATOMISTIC MODEL RESULTS CONCLUSIONS APPENDIX A APPENDIX B APPENDIX C REFERENCES

In the giant unilamellar vesicle (GUV) experiments of Evans and co-workers (Kwok and Evans, 1981; Rawicz et al., 2000; Olbrich et al., 2000), pure lipid bilayers form spherical vesicles of radius on the order of 20 µm. This experimental model for a cell membrane has proven to be extremely useful in understanding the mechanical properties of lipid bilayers. Micropipette suction of GUVs allows the measurement of various material properties; for example, the elastic area modulus, bulk modulus, and water permeability. GUV experiments have also been studied in the context of lipid domains (Bagatolli and Gratton, 1999, 2000a,b; Bagatolli et al., 2000), and as a template in which to examine lipid rafts (Dietrich et al., 2001). In both cases, it was found that the GUV may contain domains of different phases (liquid crystal/gel) (Bagatolli and Gratton, 1999, 2000b), different compositions (Bagatolli and Gratton, 2000a), and lipid rafts (Dietrich et al., 2001). The complex nature of the GUV (in terms of the various compositions and phases of different domains) is directly related to the bilayer's composite material properties. Moreover, in Gheber and Edidin (1999) a model was proposed that related membrane domains to lipid lateral diffusion.

The length-scale associated with a GUV is considered to be macroscopic relative to molecular dimensions. Typically, the thickness of a single bilayer is on the order of 30-40 Å; thus the length-scale ratio of bilayer thickness to its area is very small. Molecular dynamics (MD) simulations are not capable of modeling biological assemblies with length-scales on the order of a GUV, which is a significant drawback given the wealth of atomistic-level information which is, in principle, accessible from MD. To be able to directly correlate macro-scale membrane deformations and material properties to microscopic interactions would offer a new level of understanding. Furthermore, not only are the accessible length-scales of MD orders of magnitude below that required to model a GUV, but the time-scales of MD are presently computationally restricted to tens of nanoseconds.

Thus, an alternative approach to model a GUV and other large-scale biological assemblies is in order. Continuum-level simulation, traditionally used in engineering fields, can, in principle, be adapted to model biological structures. In continuum-level simulations the material or system of interest is represented as a homogeneous continuum with no molecular-level detail. The time- and length-scales are no longer restricted to microscopic domains. However, the cost of removing molecular-level detail is that the dynamics must be described by a constitutive relation which, at a minimum, relates stress to strain (in the case of an elastic material), or stress to strain-rate (for viscous systems).

The constitutive relation relevant to biological membranes was proposed in Evans and Needham (1987) and is expressed as

(1)

where  is the in-plane stress of the membrane,  is the bulk modulus, and  is the corresponding plane strain. Thus, to perform a continuum-level simulation of a membrane, the first task is to determine the value of the bulk modulus. In many continuum-level problems (for example, modeling an elastic band), the modulus can easily be determined experimentally and, as such, is often considered as an input constant. Experimental measurements of the bulk modulus for biological membranes have indeed been performed, and the experimental value of  might then be used. However, this "top-down" approach to molecular modeling requires macroscopic input to study macroscopic phenomena. It can be far more insightful if the bulk modulus is determined from atomistic models, and then propagated up to the macro scale. In this way, a bridging between micro-scale processes and macroscopic properties could be constructed.

In two recent papers we proposed a method whereby the bulk modulus is calculated from an atomistic model, and this density-dependent quantity is then interfaced with a continuum-level simulation (Ayton et al., 2001a,b). In Ayton et al. (2002) the bulk modulus was calculated for a real membrane (dimyristoylphosphatidylcholine (DMPC)) using nonequilibrium molecular dynamics (NEMD). In all cases it was found that the bulk modulus is dependent on density, phase, temperature, and the underlying molecular interaction model. An obvious extension to the work in Ayton et al. (2002) is to design an appropriate continuum-level membrane model that can bridge the microscopic and macroscopic domains.

However, continuum-level models are not as easily implemented as atomistic-level MD. The calculation of strains, and strain-rates, about discrete points in the material demands some method for evaluating what are essentially global properties (flow fields). An incorrect evaluation of strain, via the constitutive relation, can result in poorly sampled stresses. The situation is further complicated in a continuum-level membrane formulation, where the stress and strain must be decomposed into components parallel and perpendicular to the membrane.

A general continuum-level membrane formulation using the material point method (MPM) has been proposed and tested (York et al., 1999, 2000); however, for it to be applied to a biological membrane, special modifications must be made. First, biological membranes exhibit the property of a bulk modulus and a shear viscosity. This effect arises from the fact that the membrane behaves as a two-dimensional fluid, and thus the constitutive relation must be carefully chosen. Second, in the context of a spherical GUV, significant simplifications to the continuum-level model can be made. Almost all experiments of GUVs involve non-flaccid vesicles; that is, vesicles under a pre-stressed state. Under low osmotic stress, the vesicle (if its material properties are uniform) will adopt a spherical shape to minimize the surface area. However, as discussed in Olbrich et al. (2000), subvisible thermal oscillations will persist at low strain, and the vesicle must be pre-stressed to remove these thermal oscillations. In this pre-stressed regime, a continuum model that is 1) restricted to spherical shapes and 2) can respond to changing osmotic stresses, concentration gradients, and temperature can model a variety of experimental scenarios. Furthermore, if the model is constructed from a "bottom-up" approach, then atomistically determined information can be used to predict macroscopic behavior. The resulting macro-scale behavior can then be traced to effects propagating from length- and time-scales originating from microscopic domains.

In this paper we present such a continuum-level model for a GUV under osmotic pressure gradients. The focus at this point is on simplicity, rather than complete generality. We will show how osmotic pressure gradients can be modeled by coupling atomistically determined information to a continuum-level model. More extended continuum-level GUV models can also be developed from the one presented here. In the next section the continuum-level model is described where the relationships between osmotic pressure and plane stress are shown. In The Atomistic Model the microscopic systems and simulation methods are presented. The bulk modulus for various DMPC/cholesterol mixtures are calculated, and the methods used are briefly discussed. The Results section is divided into two subsections where the microscopic modulus results are presented, and then the micro-macro bridge to the continuum-level simulation is demonstrated.

	A CONTINUUM-LEVEL MODEL FOR GIANT UNILAMELLAR VESICLES

TOP ABSTRACT INTRODUCTION A CONTINUUM-LEVEL MODEL FOR... THE ATOMISTIC MODEL RESULTS CONCLUSIONS APPENDIX A APPENDIX B APPENDIX C REFERENCES

In a continuum-level simulation the system of interest can be discretized into a finite number of connected regions. The motion of any point is given by a continuous function defined throughout the computational domain, and the numerical solution is given by the displacements of the discretized points. The discretization of the material results in certain similarities to microscopic-level simulations (molecular dynamics); however, there are no independent "particles" and likewise the kinetic energy at some region is not to be directly associated with temperature. This latter property is a direct consequence of the macroscopic time- and length-scale associated with the continuum model. To apply continuum-level models, care must be taken to ensure that the chosen length and time-scales are beyond the fundamental molecular "coarse-grained" regime. In the case of lipid bilayers this regime is attained in a GUV, as the length-scales associated with thermal oscillations are subvisible (Rawicz et al., 2000).

There are a variety of continuum-level simulation techniques. Here we use the continuum-level simulation technique known as the material point method (MPM), which is a well-known numerical technique for solving large-deformation problems in continuum mechanics (Sulsky et al., 1994, 1995; Sulsky and Schreyer, 1996). Specific details relevant to thin membranes can be found in Ayton et al. (2001a,b) and York et al. (1999, 2000), and we will only briefly describe them here. In general, the equation of motion for a specific region in the material is given by

(2)

where is the stress tensor, F is an external force,  is the mass density, n is the number density, and a =  is the acceleration. For elastic systems, a constitutive relation that relates the stress to the strain is required to solve the dynamics. In this case, the incremental constitutive relation can be given by the relation

(3)

where the strain rate is given by the symmetric part of the velocity (or flow) gradient

(4)

Here, T is the modulus tensor and v is the flow, or velocity, at that point.

The MPM algorithm involves "mapping" material point velocities and masses to a regular grid v  v_g, m  m_g at each time step. This involves a transformation between a material point and a grid node representation. The specific tranformation functions used in this work are presented in Appendix B. With the material point velocities and masses evaluated in this new representation, new lab fixed momenta are found, p_g = m_gv_g, and then are integrated in time,

(5)

Once the grid momenta are found, they are "mapped" back to the material points, and new masses and positions of the material points are updated. Compared to finite-element methods, the mapping of properties such as mass and momenta from the material points to the lab fixed grid enables complete spatial freedom for the material points.

In the next subsections some details relevant to the problem of a GUV under osmotic stress will be discussed. This will include a method that can incorporate osmotic stress within the MPM framework and a representation of the GUV in an applicable coordinate frame. In this particular application of MPM, using a spherical-polar coordinate frame rather than a usual Cartesian system results in significant simplifications. However, certain subtleties regarding plane stress, strain, and the coupling of these to normal pressure gradients must be analyzed first. To construct a continuum-level simulation of a spherical GUV, three core parameters must be analytically expressed: 1) the excess pressure due to the osmotic stress; 2) the opposing normal force generated by the plane stress of the GUV, resulting from its curvature; and 3) a constitutive relation between the in-plane stress and strain.

Osmotic stress gradients

Under the presence of an excess internal pressure, a GUV will go from the flaccid (zero stress) state to a spherical (on continuum length-scales) pressurized structure. Experimentally, a pressurized GUV can be created by an osmotic pressure gradient (Hallet et al., 1993) where the concentration of the solution surrounding the GUV is decreased relative to the internal concentration. Initially, the solute concentration inside and outside the vesicle is given by C₀. The concentration of the external medium is then gradually reduced to C^*_sol, where C^*_sol = C₀,  < 1. The decreased concentration in the surrounding environment results in osmotic swelling of the vesicle, and the resulting expansion reduces the concentration inside the vesicle from C₀ to C* with

(6)

where r₀ is the initial radius of the vesicle, and r* is the radius at C^*_sol. The osmotic pressure across the vesicle wall is then given by

(7)

where, with C^*_sol = C₀, we can write

(8)

For a spherical GUV, this equation can also apply to a local region, and if a spherical-polar coordinate system is chosen, then the radius at a particular point r(, ) is specified by the two angles  and , as shown in Fig. 1. The osmotic pressure at (, ) is given by P(, ) = RTC₀ [(r₀/r*(, ))³  ]. In situations where large deviations from sphericity occur, this method of calculating P will break down. The computational details associated with the implementation of osmotic stress within the MPM framework is discussed in detail in Appendix A.

View larger version (10K): [in this window] [in a new window]   FIGURE 1   Spherical polar coordinate system used to model a GUV.

The normal component of plane stress

As previously mentioned, a convenient coordinate system applicable to a GUV is a spherical polar description defined by the radius r and the two angles  and . Cartesian coordinates are also possible; however, as will be seen, a number of simplifications to the solution of the dynamics occur in spherical polar coordinates. Consider a small area element of the membrane at (r, , ), given by A(r, , ) = r² sin dd. Under a net internal excess pressure, A(r, , ) will be displaced along r to A(r + r, , ). For small area elements the normal force acting on A(r, , ) is given by F_P = PA, where is the unit vector in the r direction.

However, the curvature of the membrane results in a small component of stress, , along the normal direction , and this scenario is shown in Fig. 2. This normal force F is in the opposite direction to the normal force generated by the excess internal pressure P. If the vesicle is subject to an excess internal pressure generated by an osmotic pressure gradient, then it will swell until for each A(r, , ), F = F_P. The normal component arising from the plane stress can be found from geometrical arguments and is given by

(9)

where h is the thickness of the vesicle, r is the radius, and  is the plane stress. This result can also be found from the more familiar relation  = nF, where n is the number density and is the gradient in the normal direction. The derivation is given in detail in Appendix C.

View larger version (13K): [in this window] [in a new window]   FIGURE 2   The normal stress component n due to the membrane curvature.

Stress and strain in spherical polar coordinates

The components of strain in spherical polar coordinates, and , can be written in terms of the radial component, r, and the initial radial component r₀. In this representation, {r, , }, we consider small area elements A = r² sin dd where the corresponding strain in each of the components is given by the normalized changes in length in the angular directions,

(10)

which is simply (r  r₀)/r₀ for both and . The constitutive relation given by Eq. 1 can be used in this representation for small A. The strain-rate along the  coordinate is given by as

(11)

and likewise for ,  = /r₀.

With the constitutive relation for a bilayer given by Eq. 1 (Hallet et al., 1993; Ayton et al., 2001a,b), where  = 1/2( + ) and /t = 0, we can use the expression for the strain-rates and to give an equation of motion for the plane stress in terms of the velocity of the radial component , i.e.,

(12)

By coupling the velocity in the radial direction to the plane stress derivative via the constitutive relation, forces in the normal direction due to osmotic swelling and normal components of stress can resolve the dynamics.

	THE ATOMISTIC MODEL

TOP ABSTRACT INTRODUCTION A CONTINUUM-LEVEL MODEL FOR... THE ATOMISTIC MODEL RESULTS CONCLUSIONS APPENDIX A APPENDIX B APPENDIX C REFERENCES

To solve the continuum-level equation of motion, Eq. 2, for a GUV, the bulk modulus , as appearing in Eq. 1, must be specified. In terms of atomistic quantities, the actual value of the modulus is dependent not only on the specific molecular interactions of the molecules that compose the bilayer, but also on the interactions of the lipids with the surrounding polar solvent. In the spirit of the "bottom up" approach to multi-scale modeling of a GUV, we have used NEMD to calculate  for a specific model system. The microscopically determined  includes the complex hydrophobic/hydrophilic electrostatic interactions by virtue of explicit molecular-level models. Once the value has been determined for a specific system, it is then used to bridge the microscopic and macroscopic domains. In essence, the detailed micro-scale interactions have been collapsed down and averaged to a single key material property in this case. In Ayton et al. (2001a,b), a feedback mechanism was developed whereby continuum-level information was transferred back to the micro-scale, resulting in a feedback loop. For this study we have only constructed the micro to macro bridge; however, a similar feedback scenario will be constructed in the future.

The first subsection below will present the equilibrium atomistic simulation method and will describe the three different systems of interest: a pure DMPC bilayer, a 2:1 DMPC/cholesterol mixture, and a 1:1 DMPC/cholesterol mixture. The next subsection will then present the NEMD methodology used to calculate the bulk modulus. The NEMD methodology summarized here is discussed in detail in Ayton et al. (2002).

Equilibrium simulation methods

The MD simulations of DMPC and DMPC/cholesterol membranes 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. The DMPC/cholesterol 1:1 system contained 32 DMPC and 32 cholesterol molecules, which were uniformly distributed in the plane of the membrane. The DMPC/cholesterol membrane at the 2:1 ratio was prepared using the same procedure as in the simulations of DPPC and DMPC bilayers with cholesterol (Smondyrev and Berkowitz, 1999b, 2001). Cholesterol molecules were once again distributed uniformly in the bilayer leaflets. We have chosen a large number of lipid molecules in the membrane to accommodate such an arrangement, such that the DMPC/cholesterol membrane at the 2:1 ratio consisted of 48 and 24 molecules, respectively.

Initial equilibration, which utilized the same protocol as in recent simulations of a membrane with cholesterol (Smondyrev and Berkowitz, 1999b, 2001), was followed by a 2-ns equilibrium molecular dynamics simulation at constant pressure and temperature. After ~1 ns the structure converged and configurations from the last 1 ns of this trajectory were used as starting points of the NEMD simulations. DMPC and DMPC/cholesterol 1:1 membranes were hydrated by 1312 water molecules, while the DMPC/cholesterol 2:1 system was hydrated by 1476 waters. This corresponds to 20.5 waters per lipid molecule, which is close to a full hydration limit of ~23 waters/lipid (Nagle and Weiner, 1988). The DMPC molecules were modeled using a united atom force field (Smondyrev and Berkowitz, 1999c). Cholesterol molecules were modeled using a modified AMBER force field, which was used in recent MD simulations of DPPC and DMPC membranes with cholesterol (Smondyrev and Berkowitz, 1999b, 2001). The water model used in the simulation was TIP3P (Jorgensen et al., 1983), and was chosen to maintain consistency with previous work (Smondyrev and Berkowitz, 2001), as well as to maintain the correct parametrization of the force field in Smondyrev and Berkowitz (1999c). 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 (PME) (Sagui and Darden, 1999; Essmann et al., 1995) 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.

NEMD methods

The bulk modulus was calculated using NEMD, and the method used was described in detail in Ayton et al. (2001a,b, 2002) and Hoover et al. (1980a,b). Essentially, a set of configurations sampled from an equilibrium NPT trajectory are subjected to an artificial strain-rate in the form of a deterministic oscillating area change given by

(13)

where  is the dimensionless amplitude of the oscillation and  is the corresponding frequency.

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. The implementation of this technique involves altering the equations of motion of the system in a manner similar to the methods traditionally used to maintain constant temperature and pressure in equilibrium simulation.

The modulus for materials that exhibit a linear elastic response about an initial strain can be found from the slope of  versus  about  = 0. By averaging over multiple initial equilibrium configurations, the bulk modulus for the system can be expressed as

(14)

where  ...  means an average over multiple configurations selected from an equilibrium trajectory. The NEMD strains used were chosen to be less than those experimentally determined to result in lysis, and the frequencies used had wavelengths on the order of 1 ns. The bulk modulus for pure DMPC was previously calculated in Ayton et al. (2002) and the value agrees well with experimental estimates. Thus, the pure lipid NEMD estimate for  can be viewed as a benchmark calculation, and the subsequent DMPC/cholesterol mixtures can be examined relative to the pure case.

	RESULTS

TOP ABSTRACT INTRODUCTION A CONTINUUM-LEVEL MODEL FOR... THE ATOMISTIC MODEL RESULTS CONCLUSIONS APPENDIX A APPENDIX B APPENDIX C REFERENCES

The results will be broken down here into three subsections. First, some important equilibrium results for the DMPC/cholesterol 2:1 mixture will be presented, as this particular system has not been previously reported in detail. Second, the results of the NEMD modulus calculation for DMPC/cholesterol mixtures will be presented and discussed. Finally, these results will be used in a continuum-level vesicle simulation where an osmotic pressure gradient due to an imposed concentration gradient is imposed. From the continuum-level MPM simulation, changes in vesicle radius, surface tension, and strain due to the osmotic gradient will be examined for the different DMPC/cholesterol mixtures.

The model is further extended to handle mixtures of varying DMPC/cholesterol concentrations within a single GUV. Randomly generated regions of pure, 2:1, and 1:1 DMPC/cholesterol domains will be generated on the surface of a GUV. The combination of different moduli, densities, and membrane thicknesses results in a variety of responses to osmotic pressure gradients.

DMPC/cholesterol 2:1 mixture

DMPC bilayers containing cholesterol at varying concentrations have been previously reported (Robinson et al., 1995; Smondyrev and Berkowitz, 2001; Tu et al., 1998; Nielsen et al., 1999; Pasenkiewicz-Gierula et al., 2000; Scott, 1991), and pure DMPC results have been reported in Smondyrev and Berkowitz (1999a). The present DMPC/cholesterol 2:1 mixture represents an intermediate concentration between two limiting cases: low sterol concentration membranes corresponding to a predicted mixed phase and sterol-rich membranes in a liquid-ordered phase (Lipowsky and Sackmann, 1995). In the context of the bulk modulus, experimentally it was found that the compositional dependence of the area compressibility modulus, with K_A ~ h, indicated ideal mixing. The present 2:1 system was examined to probe this effect via simulation.

The average area per DMPC headgroup can be used as a convenient measure of molecular packing in a DMPC/cholesterol bilayer. We can use a simple estimate to extract this value assuming that the effective area per cholesterol molecule is A_Chol ~ 32 Å² and is concentration-independent. Although it is likely that the average cholesterol area varies with sterol concentration, this estimate will provide a first measure of the molecular packing in the membrane. In Fig. 3 the time evolution for A_DMPC(t) shows that the area per DMPC molecule begins to stabilize after the initial 1.5-ns equilibration. The average values for A_DMPC calculated for the last 1.4 ns are 50.6 Å² and 45.5 Å² for membranes with 33 mol % (2:1 DMPC/cholesterol) and 50 mol % (1:1 ratio) of cholesterol, respectively. As the membrane area expands at lower cholesterol concentration, its thickness decreases correspondingly. The membrane thickness can be estimated as the average distance between DMPC phosphorus atoms in the opposing bilayer leaflets. The P-P distances averaged over the last 1.0 ns of simulations are 39.8 and 42.8 Å for membranes with 2:1 and 1:1 DMPC/cholesterol ratios, respectively. As a result, it would be expected that the DMPC hydrocarbon tails become more disordered.

View larger version (25K): [in this window] [in a new window]   FIGURE 3   Time evolution of the area per DMPC molecule in the 2:1 DMPC/cholesterol mixture. The upper spectrum refers to the 2:1 system, while the lower corresponds to the 1:1 case. The average values for ADMPC calculated for the last 1.0 ns are 50.6 Å2 and 45.5 Å2 for the two respective systems.

Ordering of the lipid hydrocarbon tails can be examined by S_CD, the deuterium order parameter, and it can be calculated in computer simulation via the following expression (Egberts and Berendsen, 1988)

(15)

where S_ij = 3/2 cos _i cos _j  1/2_ij and _i is the angle between the z-axis and the ith particle's molecular axis. Fig. 4 shows S_CD averages over the Sn-1 and Sn-2 chains in membranes with 30 and 50 mol % cholesterol. When compared to the pure DMPC membrane (with typical S_CD values of 0.2 in the plateau region) (Smondyrev and Berkowitz, 1999b) and membranes with low (0.11 mol %) cholesterol content (Smondyrev and Berkowitz, 1999b), both bilayers with 30 and 50 mol % cholesterol exhibit a much more pronounced ordering of DMPC hydrocarbon tails. However, the cholesterol condensing effect is weaker in membranes with a 2:1 DMPC/cholesterol ratio. The difference between the results for S_CD is somewhat more pronounced than the observed differences in membrane geometry described above. These results indicate that the membrane with a 2:1 DMPC/cholesterol should be "softer" than the one with 50 mol % cholesterol. In the following section we show that this is in fact true by calculating the bulk modulus directly.

View larger version (10K): [in this window] [in a new window]   FIGURE 4   Deuterium order parameter. SCD in the DMPC Sn-2 tails for the DMPC/cholesterol 2:1 membrane (lower curve) and for the DMPC/cholesterol 1:1 membrane (upper curve).

NEMD: the bulk modulus

The zero-frequency estimate for the bulk modulus  for the three membranes (DMPC, DMPC/cholesterol 1:1, DMPC/cholesterol 2:1 was calculated with NEMD as described above by using Eq. 14. The systems used in the previous equilibrium study were used as initial configurations for the NEMD.

The amplitude was set at  = 0.015 and the lowest frequency examined was  = 0.00314 ps¹. The frequency dependence can be thought of as the rate that the system goes from a state of zero strain to one of 2; thus, as   0, the system is being expanded at slower and slower rates. At the lowest examined frequency, the system undergoes expansion from  = 0 to  = 2 over the course of 1.0 ns. In Fig. 5 the extrapolated expansion moduli for the various systems are shown. The sizes of the symbols are roughly the size of the error bars obtained from the least-squares slope analysis of the corresponding averaged stress versus strain plot. The actual values of the moduli can be found in Table 1, and here we first comment on the relative values of the moduli. Clearly, as the concentration of cholesterol increases, so does the modulus. (However, without more detailed interim concentrations it is difficult to conclude whether the concentration dependence of the modulus is linear.) The microscopic origin of the increased "stiffness" of the membrane due to cholesterol effects is discussed in Smondyrev and Berkowitz (1999b). With NEMD we can further probe the cause of the substantially increased modulus.

View larger version (16K): [in this window] [in a new window]   FIGURE 5   NEMD extrapolations for the three different DMPC/cholesterol systems. The pure DMPC system has a zero frequency modulus of 5.4 ± 0.3 × 107 Nm2 (bottom line). The DMPC/cholesterol 2:1 mixture was found to have  = 17.8 ± 0.8 × 107 Nm2 (middle line), while the DMPC/cholesterol 1:1 system had  = 42.5 ± 1.4 × 107 Nm2.

View this table: [in this window] [in a new window]   TABLE 1   Various NEMD estimates for the expansion bulk modulus

We note that the modulus for the 1:1 system is approximately seven times that for the pure DMPC membrane. In Evans and Needham (1987) and Lipowsky and Sackmann (1995) it was experimentally found to be on the order of four to five times the size. A direct comparison is difficult, however, because different temperatures were used. Nevertheless, it appears that the NEMD estimate for  for the DMPC/cholesterol 1:1 system has somewhat overestimated the experimental value.

Given that the NEMD estimate for  in pure DMPC agrees well with experiment, there are two obvious possibilities for the large value found for the DMPC/cholesterol mixture. First, as the system is a mixture of two components, the possible spatial configurations of the two species are large; that is, for a given number of lipid and cholesterol molecules, there are a number of possible configurational starting points, ranging from a perfectly alternating intermixed array to two separated domains. In the small finite systems used in NEMD, it is almost impossible to sample all possible mixture configurations. For this work an alternating array was used, and thus the calculated modulus is representative of this and similar configurations. A "domained" system may have a different  value, and in the experimental case, by virtue of exploring large systems, essentially all possible plaquettes are sampled. Second, the spatial structure of the DMPC/cholesterol mixture contains strong solid-like correlations, included a tilt angle. To respond to changes in area via NEMD expansion, it is possible that the more "rigid" cholesterol molecules have less conformational degrees of freedom in which to respond. In other words, the pure DMPC system is more fluid-like, and the response to strain-rates is not constrained to global molecular motions. We note that under NEMD expansion no shear-stress component was observed, indicating that commensuration effects are small.

However, we feel that a calculated factor of 7 versus an experimental factor of 5 in the modulus difference is still quite reasonable in quantitative terms, so that the result obtained here is valid for the molecular models and the chosen geometry. If a more exact NEMD estimate is desirable, a series of mixture simulations with different starting conditions would need to be performed and included in Eq. 14. This will be the subject of future research.

MPM: osmotic pressure effects

The continuum-level MPM simulation parameters are summarized in Table 2. The size of the simulation in terms of the number of material points is not crucial in the case that the modulus and density are constant throughout the membrane surface. In fact, in this case all material points should be equally radially displaced under the osmotic pressure gradient. This criterion was used to test the MPM algorithm. We have chosen to use a fairly large system (in computational terms) for two reasons: first, to demonstrate the computational efficiency of the algorithm. For a simulation of this size, 100,000 time steps take ~10 min on a slow workstation (SGI 02, R5000, 200 mHz) without any parallel decomposition. Second, and more importantly, this model can be extended to vesicles where different regions of the vesicle can have different material properties, as each material point is, in essence, independent of the others. This particular point will be addressed in more detail later in the results.

Each MPM simulation was initialized from a GUV with a radius of 10 µm. At time zero the concentration ratio was gradually reduced to  = C^*_sol/C₀ for various values of  ranging from  = 0.75 ... 0.95. To mirror the experiments of Hallet et al. (1993) the concentration of the external solution was decreased from C₀ gradually to the final value C^*_sol over the course of 2 µs. We note that this is faster than is attainable in an experiment, but for this study we are only interested in the long-time steady state. In reality, a GUV responds to the altered concentration by expanding to a new constant radius (provided the strain is below the lysis point). In terms of the MPM algorithm, to obtain a non-oscillating solution, a viscous dampening term was included within the MPM time-integration, and has the effect of dampening out oscillatory motion. Details of this procedure are given in Appendix A. We note that the final steady-state GUV radius is independent of the strength of the viscous dampening so long as the dampening frequency is faster than the characteristic oscillatory solution.

In Fig. 6 the time dependence of the vesicle radius r(t) is shown for the three different DMPC/cholesterol systems with  = 0.90, corresponding to C^*_sol = 0.0036 mol/l. In all cases the vesicle responds to the altered external concentration by initially expanding until a steady state is reached. The pure DMPC system expands the most, with a final radius of 10.1 µm, while the DMPC/cholesterol 1:1 expands the least for the same concentration gradient. Interestingly, the final radius for the 2:1 system is closer to the 1:1 system than the pure DMPC membrane. The 1:1 DMPC/cholesterol mixture expands the least. The increased stiffness associated with the cholesterol concentration results in a decreased expansion for the same vesicle size and concentration difference.

View larger version (12K): [in this window] [in a new window]   FIGURE 6   Time-dependent relaxation to the steady-state radius for the three different DMPC/cholesterol lipid mixtures. The solid line is for pure DMPC, the dashed line for DMPC/cholesterol 2:1, and the dotted line is for DMPC/cholesterol 1:1. In all cases, C*sol/C0 = 0.9.

Within the MPM framework it is possible to examine the strain and surface tension that results from the osmotic pressure gradient. The final steady state for the vesicle will occur when the normal component of stress, due to the curvature of the vesicle, balances the opposing normal pressure. The surface tension predicted by considering a hemisphere of a thin spherical membrane (Hallet et al., 1993) is given by  = Pr/2, where P is given by Eq. 7 and r is the radius of the vesicle. For  = 0.90, C₀ = 0.004,  = 3.66 mN/m, while from the MPM simulation with the same parameters (Fig. 7)  = 3.66 ± 0.002 mN/m, indicating that the simulation is giving the correct solution. It was found that the point where the osmotic force was included in the MPM algorithm was crucial, and the correct solution was found when all normal forces were transformed to the grid node representation.

View larger version (11K): [in this window] [in a new window]   FIGURE 7   Time-dependent relaxation to the steady-state surface tension for the three different DMPC/cholesterol lipid mixtures. The solid line is for pure DMPC, the dashed line for DMPC/cholesterol 2:1, and the dotted line is for DMPC/cholesterol 1:1. In all cases, C*sol/C0 = 0.9.

In Fig. 8, the strain  corresponding to the radial displacement in Fig. 6 shows the effect of the increased modulus going from pure DMPC to a 1:1 mixture of DMPC and cholesterol. Although the resulting surface tension as shown in Fig. 7 is similar in magnitude ( ~ 4-5) for the various systems, the corresponding strain decreases as the cholesterol concentration increases. This is due to the increased resistance to area dilation resulting from the cholesterol concentration. For a given stress the DMPC/cholesterol 1:1 system does not expand as much as the pure system for the same osmotic conditions (C^*_sol/C₀ = 0.9 after 2 µs).

View larger version (12K): [in this window] [in a new window]   FIGURE 8   Time-dependent relaxation to the steady-state strain for the three different DMPC/cholesterol lipid mixtures. The solid line is for pure DMPC, the dashed line for DMPC/cholesterol 2:1, and the dotted line is for DMPC/cholesterol 1:1. In all cases, C*sol/C0 = 0.9.

The relative trends continue as C^*_sol/C₀ is decreased. For DMPC vesicles, lysis occurs at ~2-3 mN/m (Evans and Needham, 1987), which occurs for C^*_sol/C₀ ~ 0.9. Without further microscopic information bridged to the continuum model, the MPM model cannot predict lysis. A method for this inclusion was presented in Ayton et al. (2001a), and involves an equilibrium MD study of phase stability. From this model the resulting surface tension arising from various concentrations C^*_sol, along with the corresponding strain, can be examined. As can be seen from Fig. 9, for C^*_sol/C₀ ~ 0.9 the resulting tension is similar for all three systems. However, as C^*_sol/C₀ decreases, the surface tension of the mixtures increases to  ~ 10 for C^*_sol/C₀ ~ 0.8, while the pure system has  ~ 7.5. The resulting strain for the mixtures is also less (Fig. 10) in accordance with the smaller increase in vesicle volume due to the increased modulus.

View larger version (10K): [in this window] [in a new window]   FIGURE 9   The variation of the surface tension  as a function of C*sol/C0 for three different DMPC/cholesterol lipid mixtures. The open squares correspond to pure DMPC, the open triangles for DMPC/cholesterol 2:1, and the open circles are DMPC/cholesterol 1:1.

View larger version (10K): [in this window] [in a new window]   FIGURE 10   The variation of the strain  as a function of C*sol/C0 for three different DMPC/cholesterol lipid mixtures. The open squares correspond to pure DMPC, the open triangles to DMPC/cholesterol 2:1, and the open circles are DMPC/cholesterol 1:1.

Mixtures of domains

Each material point in the MPM simulation responds to local stresses via local strains. The result of the balance between plane stress  and the osmotic pressure P at some point r, , is a displacement of a small area element A(r, , ). The cumulative effect of all area elements on the surface of the GUV results in the total volume V of the vesicle expanding from V₀ = 4r to V_f = 4(r + r)³, where

(16)

and N_mp is the number of material points. In the case that the osmotic pressure is given by Eq. 8 and the material properties of all the material points that compose the surface of the GUV are identical (i.e., bulk modulus, density, thickness), then (r_i + r_i)³ = (r + r)³ for all i, and a uniform GUV expansion results.

In principle, for the MPM model to handle more complex scenarios, for example, large deformations due to vesicle collision, the decomposition of stress into normal and plane components would have to be performed. However, for small deformations due to small changes in material properties, the present model can still be used. (The approximation that the membrane normal at some r, , lies along the radial vector implies that no large deformations can take place.)

From the experimental work of Bagatolli and Gratton (1999, 2000a,b) and Bagatolli et al. (2000) it was observed that domains of varying composition and phase can occur within one GUV. GUVs composed of mixtures of DMPC with 30 mol % cholesterol, however, exhibited a homogenizing effect where no macro-scale gel-lipid domains coexisted (Bagatolli and Gratton, 1999) and no strong vesicle shape changes occurred. This particular scenario is possible within the framework of the present MPM continuum model. Mixtures of pure, 2:1, and/or 1:1 DMPC/cholesterol microdomains can be randomly generated on the surface of the MPM GUV, where the material properties of the domain are determined from the previously determined microscopic parameters. Situations where one system (pure, 2:1, or 1:1) is adjacent to another are handled by applying the appropriate constitutive relation as given in Eq. 1. The strain at that region is found from the MPM transformation between the grid nodes and material points (given in Appendix B). Here is an important point: the MPM grid node/material point transformation naturally "blends" the properties of the material points. For domains that are sufficiently small and randomly placed, the present MPM algorithm should therefore be able to resolve the average behavior.

To examine the mixture-MPM model, two MPM GUV mixture simulations were performed. One was constructed with a DMPC/cholesterol mix of 8.4:1, while the other was an intermediate mixture of 1.53:1; in each case they were created by generating random domains on the GUV surface. The simulation size was set at 2500 grid nodes and 39,200 material points, and was run for a total of 100 µs, 10 times longer than in the pure case. Also, a larger simulation with 6400 grid nodes and 101120 material points was also performed to examine the effect of grid spacing.

To generate the GUV mixture, small domains of pure, 2:1, or 1:1 DMPC/cholesterol were created by the following algorithm: 1) initially create a pure DMPC vesicle; 2) randomly select a region on the surface of the vesicle, and target a small surface area with radius r_target = 0.344 µm about this point; 3) select all material points within this region, and then randomly change the material properties of that point (modulus, thickness, density) to correspond to the pure, 2:1, or 1:1 case. It was found that 1000 randomly generated regions under this algorithm resulted in an 8.4:1 mixture, where the domains were randomly scattered on the vesicle surface. A snapshot of the GUV simulation is given in Fig. 11. For this mixture, the average modulus was found from _avg = /2, where  ...  is the average over all material points. For the 8.4:1 mix, _avg = 8.5 ± 3 × 10⁷ Nm² which, when compared to the values in Table 1, is qualitatively between the pure and 2:1 systems.

View larger version (184K): [in this window] [in a new window]   FIGURE 11   A snapshot of the DMPC/cholesterol 8.4:1 GUV. The inset shows the entire vesicle, and the boxed area is shown in more detail in the foreground. The dark material points correspond to a DMPC/cholesterol 1:1 system with material properties summarized in Table 1. Likewise, the white points are for the 2:1 system. The gray material points represent the pure system. For this system, the average radius is r = 10.065 ± 0.001 µm, with an average modulus of avg = 8.4 ± 3 × 107 Nm2. The fluctuations arise from the different regions with varying material properties (modulus, density, thickness).

Inspection of Fig. 11 shows that after 100 µs, small undulations in the membrane surface form. On average, the membrane radius has expanded from r₀ = 10 µm to r = 10.065 µm with a fluctuation of 0.001 µm, indicating that the perturbations resulting from the mixture extension are small. However, it is noticeable that the 1:1 regions, with a modulus seven times that of pure DMPC, clearly result in small "dints" in the vesicle surface. Of course, in a real GUV these domains have the ability to migrate, and this extension will be necessary in the future to fully model the system. Still, at this level, the qualitative behavior of cholesterol-rich domains on a GUV result in an intermediate average modulus, consistent with the 8.4:1 mixture, obtained without ever specifically calculating the modulus for an 8.4:1 system with NEMD.

Ten thousand random test regions on the surface of the GUV resulted in a 1.53:1 DMPC/cholesterol mix, and in this case _avg = 28 ± 7 × 10⁷ Nm². A snapshot of the final steady state for the 39200 material point system is shown in Fig. 12, where in contrast to the 8.4:1 DMPC/cholesterol system (Fig. 11) the GUV does not exhibit membrane surface irregularities. This is consistent with the slightly smaller steady-state radius of the 1.53:1 system (r = 10.02 µm) in conjunction with a smaller fluctuation ±0.0004 µm, indicating that the more dense and spatially uniform cholesterol-rich domains result in a behavior that is intermediate between the 2:1 and 1:1 system. Also, the lack of any undulations in the GUV surface is well within the present MPM model, as the membrane normal at each material point is very close to the radial normal.

View larger version (183K): [in this window] [in a new window]   FIGURE 12   A snapshot of the DMPC/cholesterol 1.529:1 GUV. The inset shows the entire vesicle, and the boxed area is shown in more detail in the foreground. The dark material points correspond to a DMPC/cholesterol 1:1 system with material properties summarized in Table 1. Likewise, the white points are for the 2:1 system. The gray material points represent the pure system. For this system, the average radius is r = 10.02 ± 0.0004 µm, with an average modulus of avg = 28.2 ± 6 × 107 Nm2. The fluctuations arise from the different regions with varying material properties.

It is instructive to compare the _avg as obtained from the MPM simulations to the results for different DMPC/cholesterol mixtures found from NEMD. In Fig. 13 the modulus as a function of cholesterol fraction for both the MPM and NEMD simulations are shown. The MPM GUV, constructed from essentially small domains of the systems examined with NEMD, exhibits an average modulus that appears to follow the trend sets from the detailed NEMD simulations. The large error bars from the MPM GUV simulation are a result of the small domains having different material properties. The differences in _avg for the two MPM system sizes are smaller than the symbols, indicating that the grid node spacing does not have a large effect.

View larger version (11K): [in this window] [in a new window]   FIGURE 13   A summary of the bulk modulus  versus cholesterol fraction. The open symbols correspond to microscopic NEMD calculations, while the solid symbols correspond to continuum-level MPM GUV simulations. In the MPM GUV simulation,  is an average of  and  over all material points, and the error bars represent the fluctuation arising from different domains on the surface of the GUV with different material properties.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION A CONTINUUM-LEVEL MODEL FOR... THE ATOMISTIC MODEL RESULTS CONCLUSIONS APPENDIX A APPENDIX B APPENDIX C REFERENCES

This paper has demonstrated how a GUV can be modeled at the continuum level by including atomistically determined information from NEMD simulations. Material properties obtained from detailed microscopic MD and NEMD simulations can therefore be bridged to continuum-level MPM simulations. In the pure case, the MPM algorithm can be used to examine the time-dependent response of a GUV to osmotic pressure gradients generated by external solvent concentration changes. The steady-state results are seen to agree with an analytic solution. The model can also be extended to a GUV composed of different domains. For this first study, a DMPC/cholesterol mixture was created by rando