This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.045716v1 87/5/3242    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by McWhirter, J. L. Articles by Voth, G. A. Search for Related Content PubMed PubMed Citation Articles by McWhirter, J. L. Articles by Voth, G. A.

Coupling Field Theory with Mesoscopic Dynamical Simulations of Multicomponent Lipid Bilayers

Department of Chemistry and Henry Eyring Center for Theoretical Chemistry, University of Utah, Salt Lake City, Utah

Correspondence: Address reprint requests to Professor Gregory A. Voth, Dept. of Chemistry, University of Utah, 315 S. 1400 E. Rm 2020, Salt Lake City, UT 84112-0850. Tel.: 801-581-7272; E-mail: voth{at}chem.utah.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY SIMULATION METHODOLOGY PARAMETERS RESULTS DISCUSSION AND CONCLUSIONS APPENDIX A: FUNCTIONAL... APPENDIX B: EM MODEL... ACKNOWLEDGEMENTS REFERENCES

 
A method for simulating a two-component lipid bilayer membrane in the mesoscopic regime is presented. The membrane is modeled as an elastic network of bonded points; the spring constants of these bonds are parameterized by the microscopic bulk modulus estimated from earlier atomistic nonequilibrium molecular dynamics simulations for several bilayer mixtures of DMPC and cholesterol. The modulus depends on the composition of a point in the elastic membrane model. The dynamics of the composition field is governed by the Cahn-Hilliard equation where a free energy functional models the coupling between the composition and curvature fields. The strength of the bonds in the elastic network are then modulated noting local changes in the composition and using a fit to the nonequilibrium molecular dynamics simulation data. Estimates for the magnitude and sign of the coupling parameter in the free energy model are made treating the bending modulus as a function of composition. A procedure for assigning the remaining parameters in the free energy model is also outlined. It is found that the square of the mean curvature averaged over the entire simulation box is enhanced if the strength of the bonds in the elastic network are modulated in response to local changes in the composition field. We suggest that this simulation method could also be used to determine if phase coexistence affects the stress response of the membrane to uniform dilations in area. This response, measured in the mesoscopic regime, is already known to be conditioned or renormalized by thermal undulations.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY SIMULATION METHODOLOGY PARAMETERS RESULTS DISCUSSION AND CONCLUSIONS APPENDIX A: FUNCTIONAL... APPENDIX B: EM MODEL... ACKNOWLEDGEMENTS REFERENCES

 
The coupling between curvature and internal degrees of freedom is responsible for many shape transformations in biological membranes (Lipowsky, 1991; Sackmann, 1994). One example of such a coupling is echinocytosis, where a human red blood cell undergoes a transition from a biconcave shape to a crenated shape induced by an asymmetric adsorption of drugs into the two membrane monolayers (Sheetz et al., 1976). Another example, which is perhaps more relevant to our study, is giant unilamellar vesicle (GUV) budding induced by lateral phase separations within a closed lipid bilayer membrane composed of several molecular species (Jülicher and Lipowsky, 1993; Sackmann and Feder, 1995; Lipowsky and Sackmann, 1995; Baumgart et al., 2003). During budding, the interface between the two phases is located in a narrow neck interconnecting the mother and daughter cells. Provided the energy required to form a region of curvature at one phase is not large in comparison to the energy required to have an interface between two phase domains, the overall free energy of the membrane can be reduced by having a smaller, daughter vesicle of one phase sprout from the larger, mother vesicle of the other phase. Budding is sometimes the precursor to membrane fission; as a result, this phenomenon has received much attention. Experimentally, lipid domains in vesicles have been directly visualized using electron, two-photon excitation fluorescence, and confocal fluorescence microscopy techniques (Parasassi and Gratton, 1995; Bagatolli and Gratton, 1999, 2000a,b). Furthermore, the phases themselves have been characterized by fluorescence correlation spectroscopy, which monitors the translational diffusion of a fluorescent probe, typically 3,3'-dioadecylindocarboncynine (Korlach et al., 1999).

One motivation for studying GUVs is to understand the origins of the various shape changes in a biological vesicle subjected to variations in external conditions such as temperature, and internal conditions such as composition, at the different regions in a membrane (Baumgart et al., 2003). Composition is related to the concentrations of the different constituent molecules in the bilayer. A real, biological membrane has far more structure than an artificially synthesized GUV; for example, a biological membrane is supported by a cytoskeleton. Nonetheless, because of their simple structure, GUVs have been examined to simplify the identification of important physical mechanisms that might contribute to shape changes observed in real cells. Lateral phase separation is believed to be one of these mechanisms. Recent theoretical studies performed by Taniguchi (1995) and Jiang et al. (2000) have focused on determining the dominant, equilibrium states of a two-component membrane in the macroscopic regime by minimizing a Helfrich free energy model (Helfrich, 1973) extended to include composition as a degree of freedom. Here the membrane is treated as a continuum, that is, as a two-dimensional surface embedded in a three-dimensional space. Jiang et al. (2000) found these minima by solving analytically an equation where the gradient of the chemical potential associated with the free energy model is set to zero at all points on the membrane. Alternatively, Taniguchi (1995) found these minima by numerically integrating the Cahn-Hilliard (CH) equation (Cahn and Hilliard, 1958). Here the dynamics of the composition field (the composition at all points on the membrane) is driven by a diffusive flux, which is proportional to the gradient of the chemical potential (Langer, 1971; Castellano and Glotzer, 1995; Taniguchi, 1995; Kumar and Rao, 1998; Groves et al., 2000; Jiang et al., 2000). Some mixtures exhibit phase coexistence across a range of composition values, so the free energy model is made sufficiently general to allow for the possibility of phase separation. These studies have reproduced some of the vesicle structures seen in recent experiments on three-component GUVs (Baumgart et al., 2003).

Phase diagrams deduced from calorimetric, ESR, and SANS experiments on two-component membranes composed of DMPC and cholesterol suggest that these mixtures can phase-separate into coexisting liquid crystalline and fluid phases across a wide range of cholesterol mole fractions (Recktenwald and McConnell, 1981; Knoll et al., 1985; Vist and Davis, 1990; Almeida et al., 1992). However, although lateral phase separation has been observed directly using fluorescence microscopy for binary cholesterol lipid monolayers (Subramanian and McConnell, 1987; Keller and McConnell, 1999), phase domains have not been observed for binary cholesterol lipid bilayers. Ternary bilayer mixtures composed of cholesterol along with saturated and unsaturated lipids have exhibited binary phase separation (Dietrich et al., 2001; Baumgart et al., 2003) where it has been suggested that a cholesterol molecule becomes strongly correlated or binds with the unsaturated lipids to form a complex (Radhakrashnan et al., 2000). A simple theoretical model might suppose that each cholesterol molecule irreversibly dimerizes with a single unsaturated lipid so that the membrane behaves effectively as a two-component bilayer mixture where one of the components is now a complex (Komura et al., 2004). The material properties of these bilayer mixtures in GUVs have been measured using micropipette suction or pressurization (Kwok and Evans, 1981; Evans and Needham, 1987; Needham and Evans, 1988; Mui et al., 1993). These experiments show that even small amounts of cholesterol can dramatically change the stretching and bending elasticities of the membrane (Needham and Nunn, 1990). Equilibrium and nonequilibrium molecular dynamics (EMD and NEMD) simulations of atomistic model bilayers have also shown that the microscopic bulk modulus, , is affected by the addition of cholesterol to a membrane composed of primarily DMPC (Smondyrev and Berkowitz, 2001; Ayton et al., 2002b). In particular, increases with an increase in the concentration of cholesterol. The microscopic bulk modulus is a measure of the stress response of a small microscopic volume of bilayer with dimensions 3–10 nm to a local, isotropic in-plane strain: the stress, , is related to the strain, , via the constitutive relation, = 2 (Evans and Needham, 1987). Across longer length scales, however, the stress versus strain curve of a GUV determined experimentally is nonlinear at low strains or surface area dilations (Rawicz et al., 2000). This nonlinear behavior is thought to be related to the existence of subvisible (not directly visible by microscopy) thermal ripples or bending modes (Evans and Rawicz, 1990, 1997). At low strains, the apparent, projected area is less than the actual, nonprojected area because of thermal undulations of the membrane's surface. The stress response of a mesoscopic volume of bilayer with dimensions 20–100 nm to a change in apparent area contains a contribution from the intrinsic, in-plane local area changes and a contribution from the force required to either smooth or buckle the membrane's surface. The stress response becomes linear once these soft bending modes are removed and the apparent area approaches the actual area. Therefore, the mesoscopic bulk modulus, _meso, characterizing this mesoscopic response, is not the same as the microscopic bulk modulus, : the stress response is conditioned by the existence of thermal undulations. The theoretical studies performed by Taniguchi (1995) and Jiang et al. (2000) do not treat thermal undulations; these studies do not examine how the material properties of a membrane are conditioned by thermal undulations and the coexistence of phase domains whose spatial dimensions are mesoscopic in size.

To theoretically verify whether these undulations do indeed affect the material properties of a membrane, various simulation strategies can be employed. Large-scale atomistic molecular dynamics (MD) simulations with length scales on the order of 20 nm have been performed (Lindahl and Edholm, 2000a,b). System sizes of this magnitude are required to remove finite size effects arising from periodic boundary conditions. With simulation box lengths on the order of the membrane thickness, 3–5 nm, these conditions inhibit the production of any long-wavelength, thermal bending undulations. Presently, to perform atomistic MD simulations with length scales on the order of 100 nm, electrostatic interactions are truncated (Marrink and Mark, 2001). However, it is known that transport properties such as the shear viscosity can be sensitive to these truncations (Wheeler et al., 1997). Alternatively, mesoscopic simulation methods such as dissipative particle dynamics (DPD) can examine longer length- and timescales (Espanol and Warren, 1995; Groot and Warren, 1997). However, the soft particle interactions employed in DPD allow particles to pass through one another; this might yield a large viscous response to uniform changes in area in comparison to the elastic response. Therefore, to bridge the gap between the microscopic regime, with length- and timescales on the order of 3–10 nm and several nanoseconds, and the mesoscopic regime, with length- and timescales on the order of 100 nm and 100 ns, we require a simulation method that accurately accounts for the effects of electrostatic particle interactions and accurately depicts the material behavior of the membrane subjected to various surface deformations.

One approach to overcoming these two simulation problems has been to abandon the molecular-level detail and instead adopt a coarse-grained description. The material properties of the membrane (for example, ) are measured using microscopic atomistic EMD and NEMD simulations, and subsequently employed to parameterize effective forces between points in a mesoscopic, elastic membrane (EM) model (Ayton and Voth, 2002). Within the EM model, a membrane is decomposed into pseudoparticles that do not represent coarse-grained molecular models (Shelley et al., 2001; Marrink and Mark, 2003), but rather are more like continuum-level volumes that interact via a predetermined constitutive relation. The EM model is essentially a discretized constitutive model where interactions between neighboring points reproduces the microscopic stress response, , of a small, microscopic volume of bilayer to a local, in-plane isotropic strain, . In turn, the mesoscopic stress response of a patch of membrane with dimensions on the order of 20–100 nm can be determined by evaluating a virial expression for the pressure that incorporates the effective, constitutive forces.

In this article, we are interested in examining the dynamics of a composition field integrated in time over a thermally undulating surface. The EM model represents a computationally simple means for generating these undulations. For a membrane composed of many molecular species, an EM particle (point) can carry information about the composition at a particular region in space. Composition within the model is treated as an extended degree of freedom with its own governing equation of motion, the CH equation (Cahn and Hilliard, 1958). The strength of the bonds in the EM model can then be modulated given the dependence of the microscopic bulk modulus on composition (Ayton et al., 2002a,b). Again our main objective is to demonstrate how the dynamics of a composition field can be coupled to mesoscopic membrane surface dynamics that exhibits significant thermal undulations. The EM model serves as a tool for demonstrating our strategy; the EM model can be extended or replaced by another mesoscopic membrane model provided this new model and its accompanying dynamics generates accurate thermal undulations. For simplicity and the purpose of demonstration, we treat our two-component bilayer as if it could exhibit lateral phase separation; however, we stress that this phase separation has been exhibited for cholesterol lipid bilayers with at least three components. As a result, to obtain predictive results in the future our method should be generalized to accommodate more than just two components.

The remainder of this article is organized as follows. In Theory, we present the CH equation and its application to a membrane with a quasiplanar geometry. In Simulation Methodology, we describe the simulation methodology where the curvature and composition fields are coupled within the existing framework of the EM model; this model is reviewed in Appendix B. In Parameters, we give a procedure for systematically assigning approximate/reasonable values to all parameters in the free energy model. The present model omits the effects of the solvent on the motion of the thermal undulations; we describe the possible consequences of this omission and how the values of some model parameters can be adjusted to compensate. Finally, in Results, we present our simulation results, which should highlight some notable features of the dynamics. For a simulation of a critical quench producing a phase separation, we observe the effects of the dynamical feedback between the curvature and composition fields on the overall curvature of the membrane (the square of the mean curvature averaged over all EM points) and the local curvature at one of the two coexisting phases (the square of the mean curvature averaged over all EM points within a given phase). Finally a table of symbols, Table 1, has been included for the reader's convenience.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY SIMULATION METHODOLOGY PARAMETERS RESULTS DISCUSSION AND CONCLUSIONS APPENDIX A: FUNCTIONAL... APPENDIX B: EM MODEL... ACKNOWLEDGEMENTS REFERENCES

The concentration of species (either DMPC or cholesterol), C, satisfies the following balance equation in a local coordinate frame on the membrane's surface:

(1)

Here _m is the two-dimensional differential gradient operator tangent to the surface, is a local orthonormal basis set on the surface where is the normal unit vector. In addition, is the flow or velocity of species tangent to the membrane's surface. We emphasize that a fixed point on the membrane's surface where may not correspond to a fixed point in the lab frame if the membrane is undergoing thermal undulations. As a result, with the concentration of species may appear to change with time in the lab frame even though C/t|_m = 0 in a local frame on the membrane. Given Eq. 1 and provided the total concentration is a constant, the balance equation for the mole fraction, of species is

(2)

The baricentric velocity, is zero (deGroot and Mazur, 1962). The composition of the membrane at a point on its surface is then defined as

(3)

where A is DMPC and B is cholesterol. A membrane pure in DMPC or cholesterol has respectively = 1 or = – 1; furthermore, the DMPC-rich, ld phase occurs at > 0.8 (<10% cholesterol), whereas the cholesterol-rich, lo phase occurs at < 0.4 (>30% cholesterol) (Recktenwald and McConnell, 1981; Knoll et al., 1985; Almeida et al., 1992). If the composition averaged over all points in the membrane is between 0.4 and 0.8, then our simulation method is constructed to generate a phase separation starting with a uniform composition and ending with a mixture of ld and lo phase domains. Therefore, gives a measure of not only the composition of the membrane, but also some information on the structure of the resulting phase domains. Given Eq. 2, the balance equation for the composition is

(4)

where, within the diffusive flux, J, we have a velocity If there is a net flow of material in the membrane where V_m 0, then a convective term must be added to Eq. 4 (this term is described later in Undulating Membrane; see below). Note that is the composition averaged over both monolayers or leaflets of our membrane. In some situations, for example a membrane composed of molecules with either conical or anticonical shapes (Lipowsky and Sackmann, 1995), cannot uniquely describe the state of the membrane at a point; as a result, here it is necessary to describe the composition of the membrane using two variables, ⁽¹⁾ and ⁽²⁾, which are the compositions of each individual monolayer.

Rather than including the species velocities in the dynamics, we write the flux in terms of a chemical potential, µ, via the constitutive relation (deGroot and Mazur, 1962)

(5)

Here µ = µ_A – µ_B is the chemical potential or free energy required to change the composition, specifically the energy required to remove a number of molecules of B = cholesterol and replace them with an equal number of molecules of A = DMPC. The constitutive coefficient, M > 0, is the mobility that can be expressed in terms of the self-diffusion constant of cholesterol at infinite dilution in a solvent of DMPC (see Parameters, below). The chemical potential in turn is written as the functional derivative with respect to of a free energy model, so that Eq. 4 becomes (Cahn and Hilliard, 1958; Langer, 1971; Metiu et al., 1976)

(6)

In principle, F[,H] is an effective, mesoscopic Hamiltonian obtained by coarse-graining the microscopic phase space integrated-over within the partition function, Z; the microscopic degrees of freedom spanning this space are weighted by a Boltzmann factor whose microscopic Hamiltonian represents the kinetic energies and potential interactions of the molecules in the system (Gunton et al., 1983). The partition function associated with F[,H] is

(7)

The curvature field, H, describes the texture of the membrane's surface, and is a field integral over all possible composition fields. The dominant contributions to this field integral come from those fields where F/ = µ*. In practice, F[, H] is set to a physically reasonable phenomenological model. In turn, Eq. 6 is used to model the phase separation dynamics of a membrane suddenly quenched from a high temperature to a low temperature, T. Initially the membrane composition, , will be homogeneous but then, provided lies between the coexistence curves of the ld and lo phases, a phase separation occurs that is coupled to the membrane's curvature. We assume that after being suddenly quenched the membrane reaches thermal equilibrium much faster than phase equilibrium; in other words, we assume that the phase separation is driven by gradients in the chemical potential, µ, not the temperature. The CH dynamics drives µ toward a target chemical potential, µ*. If we are primarily interested in the movement of the phase domains in response to changes in the curvature field, H, and how this movement in turn affects H, then, at late times after the quench, Eq. 6 should give an accurate picture of this movement. The parameters in F[, H], which depend on T, are assigned values commensurate with what is known about the near-equilibrium phase domain structure obtained from atomistic MD simulations and experiment.

Equation 6 is the Cahn-Hilliard (CH) equation (Cahn and Hilliard, 1958) where F[, H] = F + F_H. The contribution to F from there being a binary mixture is

(8)

whereas the contribution from the coupling between the composition and curvature fields is (Lipowsky and Sackmann, 1995)

(9)

Here _r is a reference composition defined in Parameters, below (see the discussion/text surrounding Eq. 29). The first term on the right-hand side of Eq. 8 models the energy cost for having a nonuniform composition field. The second term is a double-well potential,

(10)

which should have minima at = 0.8 (ld phase) and at = 0.4 (lo phase), and should be zero for pure DMPC, U_bin( = 1) = 0. A differential element of area on the surface of the membrane is dA_m = g(x, y)dxdy where g is the metric (strictly speaking, the square-root of the metric). For an almost flat membrane parallel to the x,y plane with no overlaps, the shape of the membrane can be conveniently described in the Monge representation, [x, y, h(x, y)], where h(x, y) is the height of the membrane surface at (x, y) measured parallel to the z axis (Safran, 1994). In this representation, the metric is g(x, y) = [1 + (_xh)² + (_yh)²]^1/2, and the mean curvature, H, is approximately

(11)

To be concise in our notation, we have rewritten h/x as _xh and ²h/x² as The cubic coupling given by Eq. 9 is only suitable if the solvent environment above the membrane is identical to that below. However, if these environments are not the same so that there is a distinction between up and down (or for a closed membrane vesicle, a difference between in and out), then Eq. 9 should be replaced with a bilinear coupling (Jülicher and Lipowsky, 1993; Taniguchi, 1995; Kumar and Rao, 1998; Jiang et al., 2000),

(12)

Since the atomistic NEMD simulations of different DMPC/cholesterol mixtures were performed under symmetric, periodic boundary conditions, the appropriate coupling model is Eq. 9. Knowing the interfacial width and line tension between the two phase domains as well as the coexistence curve, we can then assign values to all the parameters in Eq. 8; however, we treat the coupling, , as an adjustable parameter. Provided the deformations in the membrane from being flat are not appreciable, we have (see Appendix A)

(13)

where

(14)

B₁₁(x, y), B₂₂(x, y), and B₁₂(x, y) are given in Appendix A. To neglect the secondary coupling to curvature, we simply ignore the spatial derivatives of g arising from the operators and in Eq. 13 acting on g(x, y)I.

The CH equation is not the only possibility available for generating the composition dynamics. We could alternatively use the time-dependent Landau-Ginzburg (LG) equation (Metiu et al., 1976),

(15)

where > 0. Both equations result in the free energy, approaching a local minimum; however, only the CH equation will conserve the mean composition, keeping a constant. As a result, the CH equation appears at first glance more suitable than the LG equation for studying a closed system such as a membrane vesicle, but the LG equation can be extended by adding a composition-stat (analogous to a thermostat), that is, a Nosé-Hoover feedback term that approximately conserves the mean composition (Evans and Holian, 1985; Evans and Morriss, 1990).

	SIMULATION METHODOLOGY

TOP ABSTRACT INTRODUCTION THEORY SIMULATION METHODOLOGY PARAMETERS RESULTS DISCUSSION AND CONCLUSIONS APPENDIX A: FUNCTIONAL... APPENDIX B: EM MODEL... ACKNOWLEDGEMENTS REFERENCES

 
To couple the composition and the curvature fields, the key quantity to be resolved instantaneously is the membrane's height function, h(x, y). Details in the construction of the EM model are highlighted in Appendix B and explained fully in Ayton and Voth (2002). The EM simulation dynamics generates a fluctuating curvature field whose Fourier modes have an average amplitude satisfying equipartition (Helfrich, 1973; Lindahl and Edholm, 2000a; Brown, 2003). These simulations use periodic boundary conditions where the membrane's surface is approximately parallel to the x,y plane; thus, h(x, y) can be written as a Fourier series expansion,

(16)

Here and L_{x, y} are the dimensions of our central simulation box. The symbol in the sum over n_y signifies that when n_x = 0, terms where n_y < 0 are excluded from the sum. Given the positions of the EM particles (points) projected onto the x,y plane, r_xy(j), and their heights above this plane, z_j, the coefficients, c_l, in this expansion can be found by minimizing a residual (Travis et al., 1995),

(17)

This minimization (R/c_l = 0) yields a system of linear equations for the coefficients; solving this system during a simulation at every timestep is time-consuming and hence not feasible. Alternatively, if the number density (per unit area) or distribution of the EM points on the x,y plane is sufficiently dense, then the following approximation can be used:

(18)

Here the l = 0 coefficient is the mean height, h.

The composition dynamics is coupled to the curvature by using Eq. 16 to evaluate H, then entering this H in F[, H]. To couple the curvature dynamics to the composition, the spring constants of the bonds in the EM model are set, given the compositions _i of the points within this network. Specifically, h' and in (Eq. 45) are modulated based on fits to the thicknesses and microscopic bulk moduli of three DMPC/cholesterol mixtures estimated from atomistic EMD and NEMD simulations (Ayton et al., 2002b). This data was presented in an earlier article; here the data is presented again, except in units better reflecting the length- and timescales pertinent to our mesoscopic simulations (Table 2). The data for h and are well fit to f() = a( – 1) exp [b( – 1)] + c with constants a, b, and c. Increasing the concentration of cholesterol increases the conformational ordering in the DMPC hydrocarbon tails (McMullen et al., 1994; Smondyrev and Berkowitz, 2001; Ayton et al., 2002b); as a result, the membrane's thickness, h, increases and the membrane becomes more resistant to local, in-plane changes in area (Table 2).

View this table: [in this window] [in a new window]   TABLE 2  Various atomistic NEMD estimates for the membrane thickness, h', and microscopic bulk modulus,

	PARAMETERS

TOP ABSTRACT INTRODUCTION THEORY SIMULATION METHODOLOGY PARAMETERS RESULTS DISCUSSION AND CONCLUSIONS APPENDIX A: FUNCTIONAL... APPENDIX B: EM MODEL... ACKNOWLEDGEMENTS REFERENCES

d

(19)

To determine A and ², we minimize _L for a perfectly flat membrane.

The line tension, _L, is the difference in the free energy per unit interfacial length (measured along a line tangent to the interface between the two phases) and the free energy per unit interfacial length that the system would have if the properties of the bulk phases on both sides of the interface were constant up to a plane located so as to satisfy the following constraint (Rowlinson and Widom, 1982),

(20)

Performing this minimization for U_bin() symmetric at = _c, we find given H = 0 that

(21)

where the composition profile is and the interfacial width is

(22)

In deriving Eq. 21, we considered an interface/phase boundary parallel to the y direction, where are positions normal to the interface (x = 0) sufficiently far away from x = 0 that and Substituting (x) into Eq. 21, we find

(23)

Therefore, if we know _L and x, then we can determine A and ², and in turn , , ß, and . The remaining parameter, , can be set using the condition, U_bin( = 1) = 0 (Table 3), but since the composition dynamics is independent of , can be ignored. Crude experimental and theoretical estimates of the line tension yield (Lipowsky, 2002; Lipowsky and Dimova, 2003) _L 2 x 10^–12N to 10^–11N; we chose _L = 6 x 10^–12N or _L 3.595 amu · nm/ps². For the interfacial width, we chose x 5.6nm. The true interfacial width might be on the order of a few molecular diameters, so this choice is probably an overestimation. Choosing a smaller x will result in larger gradients in the composition field that will require more numerical effort to calculate accurately (that is, a smaller x requires a denser spatial grid to accurately calculate the spatial derivatives in the CH equation whose parameters are chosen using x as a reference). Since our main objective is simply to demonstrate that we have a mesoscopic simulation method that couples the dynamics of the curvature and composition fields, we chose for convenience a larger value of x than might be realistic.

(24)

provided is small and such that f + g = 4 are negligible relative to and Here D_eff() is an effective diffusion constant: D_eff() = [3² + 2 – ß]M. Remembering = (C_A – C_B)/C then noting the diffusion equation (deGroot and Mazur, 1962), valid only near equilibrium, we find from Eq. 24

(25)

(26)

If the total concentration C_A + C_B = C is constant, then further simplifications can be made; however, for the purpose of estimating M these simplifications are not necessary. In the infinite dilution limit where B (cholesterol) is the solute and A (DMPC) is the solvent, = 1, Eq. 25 shows that D_BA | ₌₁ = 0. Summing together Eqs. 25 and 26, then setting to 1, we obtain an expression for the mobility of

(27)

Here D_BB | ₌₁ is the lateral self-diffusion constant of cholesterol in DMPC. From experiment and simulation (Lipowsky and Sackmann, 1995; Smondyrev and Berkowitz, 2001), D_BB | ₌₁ 10^–5 nm²/ps; Eq. 27 then yields M 2 x 10^–8 nm² · ps/amu. Short-range diffusion methods such as quasielastic neutron scattering give lateral diffusion constants that are two orders-of-magnitude larger than long-range methods such as fluorescence (Korlach et al., 1999). One possible reason for this discrepancy is that thermal undulations make the actual distances traveled by molecules much larger than the apparent, projected distances; as a result, the diffusion constant is conditioned by these undulations with increases in length scale, or equivalently decreases in the spatial resolution of the measurements. The quantity D_BB|₌₁ entered into Eq. 27 is the short-range diffusion constant.

Another estimate of the mobility can be made by treating M more generally as a function of and then comparing the approximated CH equation, where such that f + g = 3 are ignored, to the diffusion equation, C_B/t _xy · [D_BB()_xyC_B]. This comparison yields M() ²U_bin/² = D_BB(). Replacing ²U_bin/² with its average from = 0.4 to 0.8 and D_BB() with D_BB( = 1), we have the crude estimate, M 9 x 10^–8 nm²ps/amu.

The EM model of a two-component membrane is composed of two parts: the Hamiltonian, F[, H], and the elastic energy model, 2E^l (Eq. 43), which is generalized to include cross-links between the two membrane leaflets (Ayton and Voth, 2002). Although (Eq. 45) is derived considering a perfectly flat membrane, Eq. 43 with this is applied to a nonflat membrane. Provided the peristaltic modes are ignored, the two-leaflet, EM model parameterized by can be equated to a single surface, continuum model parameterized by a bending modulus, _b, and a surface tension, _S (Safran, 1994). The Hamiltonian for this continuum model is the Helfrich free energy (Helfrich, 1973; Lipowsky and Sackmann, 1995),

(28)

where H₀ is the spontaneous curvature. To identify the cubic and bilinear coupling parameters in F[, H], we treat _b as a function of , expanding _b() at some reference composition, _r,

(29)

If the spontaneous curvature, H₀, also depends on , then the expression for ' should be modified. To relate _b to , we imagine a small, local bend in a bilayer membrane with two equal, finite principle radii of curvature where the local expansion in area of one monolayer is equal in magnitude to the local compression in area of the opposite monolayer (Safran, 1994; Sackmann, 1994; Lipowsky and Sackmann, 1995). Treating this asymmetric stretch as a bend is physically reasonable only across length scales comparable to the thickness of the bilayer, 3–5 nm. Therefore, the _b so derived is a microscopic bending modulus, not conditioned or renormalized by mesoscopic thermal undulations (Peliti and Leibler, 1985). Taking the membrane thickness, h', to be fixed, then equating the elastic free energy per unit area given by Eq. 44 to the Helfrich free energy of bending per unit area, we find _b C(h')³; depending on how the contribution to change in the elastic energy is distributed across the membrane, the prefactor C ranges from 1/12 to 1. Noting the fits to the NEMD and EMD data for and h', we then have – 2800–33,700, 1400–16,600, and 700–8100 amu · nm²/ps² for at _r = _I = 0.4, _r = _c = 0.6, and _r = _II = 0.8, respectively. The derivation above is crude, so the most we can say with some confidence is that _b depends linearly on . Curiously, the sign of appears unique. Nonetheless, in our simulations is treated as an adjustable parameter so that we can examine the effects of the dynamics coupling and H. As such, is assigned positive and negative values. The free energy model does not have to be restricted to a quadratic or cubic coupling between and H; the coupling model can made completely nonlinear by setting F_H to

The values in our simulation assigned to the parameters in the free energy model, F[,H], are summarized in Table 2. Note that the values for are 1–2 orders-of-magnitude lower than our estimates discussed above. Preliminary simulations showed that for very large || the composition of some EM points could be driven above = 1 (pure DMPC) for < 0. A more modest assignment for prevents this unphysical drift. Alternatively, we could have added to U_bin() the following term representing the Flory-Huggins entropy of mixing for a two-component mixture consisting of a cholesterol-lipid dimer and a monomeric lipid (a simple model where each cholesterol irreversibly dimerizes with a single lipid) (Komura et al., 2004):

(30)

Here A_l is the cross-sectional area of a lipid. The derivative of this term becomes infinite at = 0 and = 1, so this term acts as a caging potential. However, employing such a term necessitates abandoning our present, simple procedure for assigning , , ß, , and . In addition, M in our simulations is set four orders-of-magnitude larger than our estimate of 2 x 10^–8 nm²ps/amu (Table 2). As we will discuss in Bare EM Dynamics, below, the EM model in the absence of a solvent moves 10–100 times faster than predicted by hydrodynamic theory and observed in experiments. To compensate for the absence of viscous dampening effects to the undulating motion of the membrane's surface, the mobility, M, should be set 10–100 times larger than its estimate. Once phase separation has occurred, the phase domains are expected to move in response to the curvature in the membrane. The normal front speed of a phase boundary should be roughly proportional to the local curvature, specifically (see the end of Undulating Membranes, below). Since || is set 10–100 times smaller than our estimate, M should also be multiplied by another factor of 10–100. Therefore, to balance both the effects of the absence of a solvent and a low setting for ||, we set M in our simulations to 10,000 times its estimate, that is, to 2 x 10^–4 nm² ps/amu.

In the microscopic regime, finite-size effects arising from periodic boundary conditions (Allen and Tildesley, 1987) may influence the phase separation observed in an atomistic simulation of a binary mixture (L_x 3–5 nm L_y), specifically the fluid structure at the interface. Similarly, these conditions might even affect the values of the mesoscopic bulk moduli, _meso, estimated using EM simulations if the dimensions of the central simulation box (L_x 80 nm L_y) are not chosen sufficiently large to contain several phase domains. Given the values listed in Table 2, it is therefore prudent to check before a simulation whether the Fourier modes characterizing the phase separation of the composition field fit within this box. According to a linear stability analysis of the CH equation for a perfectly flat membrane (Binder, 1986), immediately after a critical quench to T = 308 Kelvin at = _c = 0.6, the wavenumber of the fastest growing Fourier mode in the composition field is k_fast = [U''_bin(_c)/2²]^1/2. Here U'' denotes the second derivative of U_bin with respect to . The corresponding wavelength is _fast 2/k_fast 18 nm, which is approximately four times smaller than the length of our simulation box.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY SIMULATION METHODOLOGY PARAMETERS RESULTS DISCUSSION AND CONCLUSIONS APPENDIX A: FUNCTIONAL... APPENDIX B: EM MODEL... ACKNOWLEDGEMENTS REFERENCES

 
The results are presented in three parts. In Bare EM Dynamics, we examine the dynamics of the bare EM model where the composition field is static and uniform. Next we consider the composition dynamics integrated in time over a configurationally frozen EM snapshot. Finally in Undulating Membranes, we combine the two dynamical schemes propagating in time both the curvature field and the phase-separating composition field.

Bare EM dynamics
Including the dissipative and random forces in the EM equations of motion will generate a canonical distribution and thus correct thermodynamic averages (Espanol and Warren, 1995; Groot and Warren, 1997). However, whereas the effects of the membrane's hydration layer on the stress response are embedded in the microscopic bulk modulus, , an explicit, mesoscopic model for the solvent is still required to reproduce the dampening effects of the bulk solvent on the membrane's motion. We therefore suspect that the thermal undulations of the EM model, being in a vacuum, move faster than those of a membrane immersed in a bulk solvent. To check this motion, we performed an EM simulation with a homogeneous composition field, = _c = 0.6, where the area, A = L_xL_y, was fixed. The bond strengths were set according to the fits to the EMD and NEMD data for _m, h', and :_m = 722.835 amu/nm³, h' = 3.74 nm, and = 58.09 amu/(nm · ps²). The simulation was started from the final configuration with dimensions L_x 80.30 nm L_y taken from a zero surface tension (_S = 0) simulation of a pure DMPC EM model ( = 1). The system was allowed to equilibrate for 500,000 timesteps of size t = 0.04 ps, followed by a production run of 1,000,000 t. The time autocorrelation functions of the undulation modes were averaged over 20 blocks of 100,000 t. The autocorrelation functions of those modes whose wavenumber, k_l, lay within the same interval, (p–1)k to pk, were then averaged:

(31)

Here p is an integer running from p = 2 to p = 10, k = k_max/10, and ·_T denotes a time average. The maximum value of k_l is where n_max = 4 is the maximum value chosen for n_x and n_y in the Fourier series expansion, h(x, y, t), Eq. 16.

Standard hydrodynamic coupling theory (Kramer, 1971; Crilly and Earnshaw, 1983; Hirn et al., 1999; Lindahl and Edholm, 2000a) predicts an oscillating regime at low wavenumbers, k, characterized by damped, oscillatory motion of the undulating modes, and an overdamped regime at high k characterized by overdamped, nonpropagated motion. This theory treats only the bending modes in the membrane, neglecting any possible coupling to the internal peristaltic modes that describe fluctuations in the membrane's thickness. The crossover between the two regimes occurs at a critical wavenumber, k_crit, which satisfies where _{S, eff} = _S + _bk². Here _w and _w are, respectively, the mass density and shear viscosity of bulk water, _b is the bending modulus of the membrane, and _S is the surface tension of the membrane. Assuming _b 4 x 10^–20 J, _w 0.891 x 10^–3 kg · m^–1 s^–1, and _S 0.1 x 10^–3 to 1.0 x 10^–3 Nm^–1, we find k_crit 10³ to 10⁴ cm^–1. The wavenumbers, k_l, that we sample in our simulation range from 10⁵ to 10⁶ cm^–1. Since the initial EM configuration was selected from a zero-surface-tension, mesoscopic simulation (Ayton and Voth, 2002), the surface tension in our simulation is small, lying within 0.1 x 10^–3 and 1.0 x 10^–3 Nm^–1. Therefore, according to this estimate of k_crit, our simulation should be within the overdamped regime provided we include a mesoscopic model dynamics for bulk water. Fig. 1 shows the normalized autocorrelation function for several p. The function for p 10 undergoes oscillations. In the oscillating regime, the period of the oscillation in C_k(t) should be T₀ = 2/₀ where ₀ = [_{S, eff}k³/(2_w)]^1/2. As shown in Fig. 1, the first peak in C_pk(t) for p = 8, 9, and 10 match the predicted times T₀ 700 ps, 500 ps, and 400 ps. For p > 10, C_pk(t) becomes a monotonically decreasing function of time. Therefore, our simulation actually lies on the crossover between the two undulating regimes. Despite this discrepancy, our simulation results are still consistent with the hydrodynamic theory. According to this theory, k_crit should increase as _w, and thereby the dampening effects of the bulk water are decreased.

View larger version (19K): [in this window] [in a new window]   FIGURE 1  Time autocorrelation functions of the coefficients in the Fourier series fit to the membrane's height, h(x,y,t). The composition of each point in the EM network is set to = c = 0.6. Here kmax 4.44 x 106 cm–1 and k = kmax/10. The correlation function labeled k* corresponds to a Fourier mode with a wavenumber 6.3 x 106 cm–1. As the wavenumbers increase above kmax, the undulating modes go from being damped to overdamped. Persistent oscillations (not shown) occur for p 5. These oscillations are particularly apparent if a single (no averaging over modes) Fourier mode is inspected.

pk

pk

Static membrane
We simulate the phase separation of a critically quenched binary mixture by integrating the CH equation of motion for the composition field, . Since the membrane is static, so that the curvature field, H, does not vary in time, we can isolate the dynamics of . In addition, we can determine the timescales required for to reach a "steady state," a phase-separated mesoscopic state where the free energy,