Originally published as Biophys J. BioFAST on March 25, 2005.
doi:10.1529/biophysj.105.059436
Biophysical Journal 88:3855-3869 (2005)
© 2005 The Biophysical Society
Coupling Field Theory with Continuum Mechanics: A Simulation of Domain Formation in Giant Unilamellar Vesicles
Gary S. Ayton *,
J. Liam McWhirter *,
Patrick McMurtry 
and
Gregory A. Voth *
* Center for Biophysical Modeling and Simulation and Department of Chemistry, and Department of Mechanical Engineering, University of Utah, Salt Lake City, Utah
Correspondence: Address reprint requests to Professor Gregory A. Voth, Tel.: 801-581-7272; E-mail: voth{at}chem.utah.edu.
 |
ABSTRACT
|
|---|
Domain formation is modeled on the surface of giant unilamellar vesicles using a Landau field theory model for phase coexistence coupled to elastic deformation mechanics (e.g., membrane curvature). Smooth particle applied mechanics, a form of smoothed particle continuum mechanics, is used to solve either the time-dependent Landau-Ginzburg or Cahn-Hilliard free-energy models for the composition dynamics. At the same time, the underlying elastic membrane is modeled using smooth particle applied mechanics, resulting in a unified computational scheme capable of treating the response of the composition fields to arbitrary deformations of the vesicle and vice versa. The results indicate that curvature coupling, along with the field theory model for composition free energy, gives domain formations that are correlated with surface defects on the vesicle. In the case that external deformations are included, the domain structures are seen to respond to such deformations. The present simulation capability provides a significant step forward toward the simulation of realistic cellular membrane processes.
|
INTRODUCTION
|
|---|
Recent experimental evidence has shown the existence of fluid-fluid phase coexistence in the form of dramatic domain structures in giant unilamellar vesicles (GUVs) (Baumgart et al., 2004
; Veatch and Keller, 2002, 2003; Veatch et al., 2004). In the case where ternary mixtures are considered, composed of cholesterol, sphingomyelin, and dioleoylphosphatidylcholine (DOPC) (Baumgart et al., 2004), cholesterol, dilauroyl phosphatidylcholine (DLPC), and dipalmitoyl phosphatidylcholine (DPPC) (Korlach et al., 1999; Feigenson and Buboltz, 2001), or cholesterol/DOPC/DPPC (Veatch and Keller, 2002, 2003; Veatch et al., 2004), the resulting fluid-fluid domains exhibit clear structure. The correlation with global curvature seems to be weak, and it is believed that line tension is more likely the key player in determining the domain sizes and shapes (Baumgart et al., 2004). Furthermore, in some cases, the domains are accompanied by distinct deformations on the vesicle surface, often in the form of circular bulges (Baumgart et al., 2004; Veatch and Keller, 2003). On the other hand, gel-liquid crystal domain coexistence has been observed in DPPE/DPPC mixtures (Bagatolli and Gratton, 2000), where, qualitatively, the shape of the domain is not highly correlated to the shape, or curvature, of the vesicle. Rather, the gel DPPE domains appear to be more painted on the surface of the vesicle. DPPC/DLPC gel-liquid crystal domains also seem to exhibit this behavior (Feigenson and Buboltz, 2001; Korlach et al., 1999), where the inclusion of cholesterol eventually results in solubilization of DPPC into the fluid phase, until the entire vesicle consists of only one fluid phase.
Domain formation on vesicles has also been examined theoretically (Seifert, 1993; Julicher and Lipowsky, 1993; Taniguchi, 1996; Jiang et al., 2000). Depending on whether the theoretical model is based explicitly on line tension (Julicher and Lipowsky, 1993) or composition (Taniguchi, 1996; Jiang et al., 2000; McWhirter et al., 2004), different results are predicted. Models based on line tension require that the system be pre-phase-separated into well-defined domains. From there, this particular free-energy framework can predict whether or not bulging or budding will form. The budding, or at least bulging, deformation is found from minimizing the free-energy subject to bending energy, line tension, and the constraint that the area of the domain is constant. Conversely, in the case that the free-energy model is based on composition, for example a Landau model for phase coexistence (Jiang et al., 2000; Van, 2002), then the model has the capability of predicting whether or not domains will form subject to the geometrical constraints of the problem. In this framework, domain formation can be further extended to include composition-curvature coupling, where the local curvature can perturb domain formation, along with curvature-composition coupling, where the local composition can then alter the local membrane structure. In the case where the latter type of coupling is restricted to variations in the membrane's material properties (i.e., in terms of a composition-dependent bulk and bending modulus), the result is soft and stiff regions on the membrane surface, depending on the location of the domains. This scenario has been examined in mesoscopic regimes (McWhirter et al., 2004) where the effects of thermal undulations were also considered. The drawback of employing a Landau model for composition is that, without additional coupling terms, it cannot directly influence the structure of the underlying membrane. In other words, it can predict whether or not domains will form, but once the domains have annealed, it cannot predict whether or not budding, or bulging, will occur. This fact is a simple consequence that the Landau model for composition, as it stands, will not result in specific composition-dependent stresses acting on the membrane.
Interestingly, real systems seem to exhibit behavior that is spanned by both theoretical frameworks: sometimes the formation of a domain is accompanied by an obvious shape change of the vesicle (for example, a bulge), and sometimes it is not. Thus, at the moment, two apparently different theoretical frameworks are required to model domain formation on vesicles. First, one model is required to predict whether domains will form, and then a second model is required to predict whether or not the resultant domains will alter the vesicle shape.
Regardless of the specific theoretical framework, a key limitation of these types of theoretical approaches is the evaluation of the actual free-energy functional itself. In the case of a Landau model for composition, only in relatively simple geometries can analytic solutions for the free-energy minima be found (Jiang et al., 2000), which is problematic if the GUV undergoes deformations. As such, the predictive power of these theoretical frameworks can be severely limited by the imposed geometrical constraints. If some means of efficiently evaluating these free-energy functionals in a completely unrestricted geometry could be found, then it may be possible to examine how domains form and couple to the shape of the underlying vesicle.
A tempting route of action is to model the formation of domains on GUVs employing molecular dynamics (MD) simulation, as this would, in principle, be able to distinguish, and perhaps even validate, the form of the ideal free-energy model framework by which to describe domain formation. However, examining domain formation on GUVs from MD is presently impossible as the lengthscales of GUVs are typically in the range of 20 µm, with a rough estimate of 
109 lipids constituting the vesicle. As such, even with coarse-grained MD methods (Marrink et al., 2004; Shelley et al., 2001; Rudd and Broughton, 1998; Kumar and Rao, 1998; Kumar et al., 2001; Laradji and Kumar, 2004), systems such as these are far beyond attainable simulation system sizes. Furthermore, the timescales, on the order of up to seconds, simply cannot be reached with present computational power and algorithms.
Another option is to employ continuum-level mechanics, where at least the underlying membrane dynamics, in the absence of composition field dynamics, can be examined using continuum-level modeling methods such as the material point method (Ayton et al., 2002b; York et al., 1999). However, once again, geometrical constraints are still a problem, which arise because of the grid-based framework employed in many continuum-level algorithms. In the material point method, the grid is used as a computational scratch-pad to evaluate continuum-level strains and strain-rates, and can be visualized as a three-dimensional lattice of points spanning the accessible space of interest. When thin structures like membranes are embedded within the grid, the transformation of in-plane quantities like the plane stress and strain to the grid can break down. As such, this type of scheme primarily works when the grid is actually bound in the plane of the membrane (Ayton et al., 2002b). Still, even though some computational problems exist with modeling membranes at the continuum level, the governing constitutive relation for membranes has been well studied (Hallet et al., 1993; Needham and Nunn, 1990; Rawicz et al., 2000; Olbrich et al., 2000), and, in the case of small deformations, it employs an elastic bulk modulus and a viscous shear viscosity (Evans and Needham, 1987; Ayton et al., 2002a). This combination of an elastic (solidlike) material property and a viscous (fluidlike) component is crucial for the membrane to perform all of its key functions (e.g., ion transport, lipid diffusion, lysis, fusion). In the case of very simple deformations (for example, the swelling of a vesicle due to osmosis or small surface deformations), an elastic constitutive model can be used (Ayton et al., 2002b). However, with a simple continuum-level membrane model, incorporating domains is still very difficult. At best, different prespecified regions on the surface of the GUV can be given different material properties corresponding to different domains. These can include, for example, a composition-dependent bulk modulus, thickness, and density. However these domains are fixed on the surface of the vesicle and, as such, they cannot move, coalesce, or change shape.
Another significant problem is that, without additional effort, continuum-level models make no contact with the underlying molecular-level interactions that ultimately determine the details of both the shape and the properties associated with domain formation on GUVs. A multiscale framework is required where the atomistic spatial temporal regime is connected in some way to the continuum-level. At its most minimal level, the bridge can be accomplished when, for example, the material properties such as the bulk modulus vary depending on the local composition, where the exact dependence is determined from atomistic-level non-equilibrium molecular dynamics (NEMD) simulations (Evans and Morriss, 1990; Ayton et al., 2002a; Ayton and Voth, 2004). In the context of our previous multiscale work (Ayton et al., 2002a,b; Ayton and Voth, 2002, 2004) this approach offers a computational means of bridging microscopic representations of the system with other representations that operate in higher length- and timescales. In the case of vesicles under osmotic stress, the lack of high amplitude thermal undulations (Marrink and Mark, 2001) allows for a direct atomistic, to continuum, bridge (Ayton et al., 2002b).
A Landau model for composition (Taniguchi, 1996; Jiang et al., 2000; McWhirter et al., 2004) can be incorporated into a multiscale simulation scheme by coupling the composition dynamics to the underlying membrane dynamics. The basic idea is that the composition fields are coupled with local surface deformations, and then the resulting domains affect the underlying membrane structure. The degree of this coupling can be varied, but the result is that domains of varying composition are not fixed on the surface of the vesicle, but can instead move, reform, and change shape. In fact, the formation of buds (Lipowsky and Sackmann, 1995; Seifert, 1993; Julicher and Lipowsky, 1993) is a classic example of how domains can couple to the underlying membrane structure.
Formally, at the continuum-level, the domains are defined by regions of a specific composition denoted by a variable 
, whose dynamics, i.e., d/dt, can be governed by either the time-dependent Landau-Ginzburg (LG) (Metiu et al., 1976; Van, 2002) or the Cahn-Hilliard (CH) equation (Metiu et al., 1976; Langer, 1971; Cahn and Hilliard, 1958). To bridge the underlying continuum-level membrane dynamics with the composition dynamics, a careful examination of the couplings between the two must be carried out. In fact, this approach has already been developed in the case of a small patch of membrane in the x,y plane (McWhirter et al., 2004), where the lengthscales were such that a mesoscopic membrane model (Ayton and Voth, 2002) could be employed. The multiscale bridge was accomplished by employing key results as found in previous work (Ayton et al., 2002b), to parameterize the mesoscopic model (Ayton and Voth, 2002). The planar geometry of the system (with the addition of small thermal undulations) allowed for the composition dynamics to be resolved via a Fourier transform method, such that the domain dynamics could be essentially be projected onto the undulating membrane surface. However, beyond this non-overlap geometry (the so-called Monge representation; Lin and Brown, 2004; Ayton and Voth, 2002), this scheme becomes exceedingly difficult.
To summarize the preceding discussion, if domains are to be successfully modeled on GUVs and other complex membrane surfaces, a computationally tractable continuum-level scheme is required for both the underlying membrane dynamics and the composition dynamics, with full coupling between the two. Fortunately, there exist alternative grid-free continuum-level methodologies known as the smooth particle methods, in particular smoothed particle hydrodynamics (SPH) (Lucy, 1977; Monaghan, 1992; Ellero et al., 2002) and smooth particle applied mechanics (SPAM) (Kum et al., 1995; Hoover and Hoover, 2003; Hoover and Posch, 1996). These schemes avoid many difficulties associated with traditional continuum-level simulation methods. SPH was originally developed to examine large-scale continuum-dynamics problems (in fact, astrophysical scenarios) (Lucy, 1977) and efforts have been made with SPH to improve the accuracy and stability of the method (see, e.g., Bonet and Kulasegaram, 2002; Bonet et al., 2004). On the other hand, the development of SPAM has tended more toward smaller-scale hydrodynamics, and in fact makes an elegant link with molecular dynamics (Hoover and Hoover, 2003). The correction schemes that appear in the SPH method (Bonet and Kulasegaram, 2002; Bonet et al., 2004) (typically associated with improving the numerical accuracy) have generally not been incorporated.
With SPAM, a continuum object can be thought of as being partitioned into a number of separate subsystems, such that each subsystem can be described by a state of local thermodynamic equilibrium (Evans and Morriss, 1990). In the case of a GUV with a diameter of 20 µm, for example, the surface of the vesicle may be partitioned into a large number of smaller areas, where each area is still, in itself, a large system relative to atomistic scales. The concept underlying SPAM is to then formulate the required continuum-level conservation equations for this partitioned system so the result is essentially a set of interacting free-energy particles that possess not only mass, position, and velocity, but possibly various other properties, for example composition, gradients of composition, and chemical potential. In a Lagrangian scheme, the SPAM particles not only translate according to local momentum conservation, but they can exchange various properties between one another. In the case of composition dynamics (either LG or CH), the SPAM particles would exchange composition, according to the underlying free-energy functional that governs the system. As such, SPAM can be thought of as both a computationally efficient means to solve continuum-level equations and as a conceptually attractive framework in which to recast continuum-level problems. In fact, if carefully formulated, the entire problem of coupled domain and vesicle deformation dynamics can be recast in a SPAM representation, and this is the topic of the present article. In doing so, almost all of the issues associated with the evaluation of complex free-energy functionals and grid-based schemes are avoided. Of course, the SPAM particles do not represent molecules, or anything remotely close. In the case of a membrane, they are best thought of as thin discs of mass that can exchange properties based on an underlying free-energy functional. It should be noted that solvent is included implicitly in the present model (i.e., bilayer properties), but explicit solvents (e.g., hydrodynamic effects such as fluid flows) can be readily included since SPAM has its origins in SPH. The coupling of a SPAM vesicle to a viscous SPAM fluid is clearly possible.
This article will therefore define a methodology to examine domain formation on GUVs coupled to its elastic deformations. The computational methodology will employ either LG or CH dynamics to model the composition dynamics. A fairly generic Landau model for composition (Taniguchi, 1996; Jiang et al., 2000) will be utilized as the underlying choice for the free-energy functional, in order to examine the form of the domain structures that emerge when the free-energy model is allowed to explore and couple with the underlying membrane deformations. In this approach, the composition fields are allowed freely flow over the surface of the deforming membrane so that they can find a free-energy minimum commensurate with the evolving and coupled geometry of the system. The entire problem, both composition and membrane dynamics, is recast with SPAM, resulting in a unified continuum-level description of the GUV system that can also be linked, in a multiscale sense, to atomistic-level properties.
|
COMPOSITION IN BINARY AND TERNARY MIXTURES
|
|---|
For a system of two components labeled a and b, the composition variable is defined as
 | (1) |
where 
a is the mass density of component a, i.e., a = ma/
V, and ma is the mass of component a that is found in the small volume V. The case where V is constant will be considered here. Likewise, the mass density for component b is b, and the obvious condition that = a + b.
In the case where 
/t = 0, one can write
 | (2) |
where d/dt is the Lagrangian time-derivative of the composition and u is the flow.
In the case where the system is a ternary mixture of components a, b, and c, which phase-separates into an effective binary mixture, a similar definition of composition can be derived. Considering a situation similar to what was observed in Veatch et al. (2004), two phases, 
and ß, are defined where the mass density of the phase is given by 
and 
is the mass density of the a-component that exists in the -phase and likewise for the components b and c. Furthermore, the total mass density of component a in the two effective phases is given by 
and likewise for components b and c. In a similar manner, the mass density in the ß-phase is defined as 
The composition variable for this effective binary system is given by
 | (3) |
Again, the total mass density is = + ß, and in the case where /t = 0, one again arrives at Eq. 2.
|
A LANDAU MODEL FOR PHASE SEPARATION
|
|---|
A generic Landau model will be employed in the present work to describe phase separation in a binary mixture (Van, 2002; Taniguchi, 1996; Jiang et al., 2000; McWhirter et al., 2004),
 | (4) |
where FM corresponds to a Helfrich bending energy and local area dilation free-energy contribution (den Otter and Briels, 2003; Lin and Brown, 2004; Brown, 2003; Brannigan and Brown, 2004; McWhirter et al., 2004), given by
 | (5) |
where kc is the bending modulus, H is the mean curvature, h is the membrane thickness, 
is a local bulk modulus, and 
refers to the local plane strain. It should be noted that this functional is appropriate for liquid and possibly gel bilayer phases, but not for the solid phase. In Eq. 4, F is the standard Landau model for phase separation (Taniguchi, 1996; Jiang et al., 2000; Van, 2002; McWhirter et al., 2004)
 | (6) |
where 
2 gives the strength of the nonlocal gradient term. Since we will specifically be dealing with membranes, dr = dA where dA is a local area element on the membrane evaluated in the correct in-plane reference frame. In cases where the membrane surface has complex undulations, or is an enclosed surface, the evaluation of these integrals can be quite complex (Jiang et al., 2000; McWhirter et al., 2004). In fact, one of the limiting factors in applying such free-energy models for membranes is evaluating these integrals on complicated, and perhaps even time-dependent, surfaces (Jiang et al., 2000).
Returning to the free-energy model, here V() is a simple double-well potential given by
 | (7) |
where n > m (both are positive), and a and b are constants. This strictly local term drives the composition within some area element dA to one of the minima in the potential. Other, more complex, free-energy models can be employed (McWhirter et al., 2004).
The functional F in Eq. 6 is in fact a free energy even though it is sometimes referred to as an effective Hamiltonian (McWhirter et al., 2004). To appreciate this distinction, one can consider a partition function Z for a microscopic system given by
 | (8) |
where 
= {r1,r2,r3,...rN,p1,p2,p3,...pN}, ri is the position of atom or molecule i, pi = mvi is the corresponding momentum, and C equals (m/2
2ß)3N/2. The atomistic-level Hamiltonian is given by 
and ß = 1/(kBT) where kB is Boltzmann's constant and T is the thermodynamic temperature. We now imagine coarse-graining the spatial extent of the system into M cells centered at locations {R1, R2, R3,...RM}. Note that R1, the location of the cell 1, is different from r1, the location of particle 1.
The composition of one of the cells, for example cell i, as determined from the real microscopic system, is denoted by 
where only the molecules within cell i contribute the local coarse-grained value of the composition. One can integrate over all compositions in cell i such that
 | (9) |
where the notation i = (Ri) is employed here, and (Ri) is an integration variable. The partition function Z can be rewritten as
 | (10) |
where 
d1d2d3...dM. However, except for simple lattice and spin models, F cannot be easily evaluated (Mazenko, 2003). In practice, one employs a model F whose parameters are chosen to reproduce experimentally determined phase structures at a specified resolution of measurement, for example, the width of a interface between two phases. The perceived or measured width at the mesoscopic scale will be different from the true microscopic width due to thermally induced variations of the location of the true interface. The formal procedure above, although not easily performed for real systems, does illustrate that F is a scale-dependent free-energy functional. That is, the values of the parameters that enter into the model F will depend on the size of the cells centered about the positions Ri (i.e., the degree of coarse-graining). Importantly, the integral in the last line of Eq. 10 is performed over all composition fields; however, below a critical temperature, those composition fields that correspond to a phase-separated configuration will receive the largest Boltzmann weight, and so are dominant in the contribution to Z.
The nonlocal behavior of F occurs through the gradient-dependent composition term (Van, 2002), and, roughly speaking, drives the system to a uniform state of composition. As such, this term will be denoted by Fmix. Likewise, the term that contains V(), since it favors phase separation, will be denoted as Fdemix. Thus, in this Landau model the free-energy minimum is found from the balance between nonlocal mixing and local demixing contributions to the free energy.
The final term in Eq. 4, F, H, couples the composition to the curvature H via
 | (11) |
where, in contrast to previous theoretical studies (Taniguchi, 1996; Jiang et al., 2000), a quadratic curvature coupling is employed here which can be justified by including a linear composition dependence to the bending modulus, i.e., k() = kc + k, where kc is the usual one-component bending modulus that enters into FM (Sackmann, 1994; Marrink and Mark, 2001; Lindahl and Edholm, 2000; Ayton and Voth, 2002) and k = 
. No other assumptions were employed in determining the functional form of the curvature coupling. It should be noted that under spherical geometries the result is that the membrane has no incentive to bend in either direction (i.e., bulge in, or bulge out). This quadratic form is in contrast to the more tradition linear coupling (Taniguchi, 1996; Jiang et al., 2000), which can result in dents and bulges depending on the local value of the composition. With the linear coupling model, a region with a negative composition will favor dents (i.e., H > 0), whereas regions with positive composition will favor bulges (i.e., H < 0). Although this form of the coupling does indeed result in interesting phase behavior, the justification of the linear form must be traced back to atomistic-level phenomena (i.e., explaining in terms of lipid structure, asymmetry of the bilayer, etc., why the domain with negative composition prefers dents in the GUV surface rather than bulges). The same situation applies for the domain with positive composition. In the present model, curvature coupling simply arises from the fact that the domain with the smaller bending modulus will have a correspondingly smaller free-energy cost to supporting a certain square of the curvature. Thus, the domain with the larger bending modulus will favor regions of smaller curvature, regardless of the sign. This form of curvature coupling, when tied to the underlying elastic membrane dynamics, results in a positive feedback scenario where domains with a smaller bending modulus collate in regions of locally higher curvature. The material properties of the domain are now modulated, and in this case, the curvature can actually be enhanced due to the local softening of the membrane.
|
COMPOSITION DYNAMICS
|
|---|
In this section the time evolution of the composition dynamics will be discussed. Both the time-dependent Landau-Ginzburg (LG) equation, and Cahn-Hilliard dynamics (CH) will be employed.
Landau-Ginzburg dynamics
The LG dynamics (Van, 2002; Metiu et al., 1976; McWhirter et al., 2004) for the composition field, , is given by
 | (12) |
where F[, H] = F + F, H is the free-energy functional and the phenomenological coefficient, , is positive. Under LG dynamics, a system should eventually reach the free-energy minimum; however, the actual dynamics are best thought of as relaxational, and as such, rigorously determining for a particular system can be difficult. This equation of motion will drive the chemical potential of the composition field, µ = F/, to a target µ*, which is the chemical potential of the environment. Once µ = µ* a state of equilibrium is achieved, and the dynamics stops; however, the mean composition, 
, is not conserved under this dynamics.
To evaluate the functional derivative, F[, H]/, care must be taken when considering the gradient term in Fmix, since strictly speaking this gradient is constrained to be in the plane of the membrane. Thus, if the gradient were evaluated in an unconstrained fashion, the normal component could appear. When evaluated in the plane of the membrane, and with b = a, the free-energy functional derivative, F[H]/, is given by
 | (13) |
where the gradient, 
, is the required in-plane gradient regardless of the local orientation of the membrane. Of course, if this functional derivative were to be evaluated in a lab-reference frame, then much more work would be required.
In the case that is non-zero, the dynamics must be constrained to conserve composition. Thus, in a Lagrangian form the constrained LG equation of motion can be written as
 | (14) |
where constrains the total composition of the system to be constant. This composition-stat, , is related to µ* via (t) = – µ*. The means by which such a constraint can be implemented will be discussed later.
The next step in the dynamical evaluation is to simplify the number of free parameters (i.e., , 2, a, and ). One option is to define a new set of scaled parameters as
 | (15) |
where now all the strengths of the different components are expressed relative to the mixing term strength. This new parameter set is not unique, but it manages to factor at least one set of parameters (i.e., 2) into the LG prefactor, . In the case that a dynamical simulation methodology is employed (as will be described later), the prefactor * is essentially combined with the fundamental timestep of the simulation. A value of * = 2 µm2/µs, combined with the timestep t (as given in Table 1), resulted in stable dynamics in the present application. Much larger values of *, for the given value of t, resulted in the total composition not being conserved (even in the presence of the composition-stat). We note that the fundamental timestep t was selected based on the underlying membrane dynamics.
It should be noted that in McWhirter et al. (2004) a significant effort was made to parameterize the model to a known system (Ayton et al., 2002b). In the present study, such a detailed parameterization will not be made, as here we are interested in more generic phase behavior. However, in the spirit of a multiscale methodology, effort will be made to assign reasonable values to specific parameters (i.e., bulk moduli, curvature coupling strengths). Also, since the lengthscales for this system are on the order of micrometers, the coupling between thermal undulations and composition do not come into play, as any thermal undulation modes are subvisible, and thereby effectively renormalize the material properties.
, is used to define the continuum mass distribution. Here, the mass density at some position r, is given by
A = A – A0 and A0 is the initial area. In this expression, a local coordinate frame is defined where the membrane normal is along the local z axis, and thus the x and y axes are orthogonal and in the plane of the membrane. In the case of a perfectly spherical vesicle for example, the local z axis would lie along the radial vector of the vesicle. Likewise P corresponds to the negative of the plane stress. The bulk modulus is given by 
, 
is equal to 
TOP
ABSTRACT
INTRODUCTION