Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2929801

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Two-dimensional membrane
- 3. Governing equations
- 4. Numerical method
- 5. Summary
- 6. Numerical results
- 7. Conclusions and future work
- References

Authors

Related links

J Comput Phys. Author manuscript; available in PMC 2010 August 30.

Published in final edited form as:

J Comput Phys. 2010; 229(1): 119–144.

PMCID: PMC2929801

NIHMSID: NIHMS224579

Jin Sun Sohn,^{a} Yu-Hau Tseng,^{a,}^{b} Shuwang Li,^{c} Axel Voigt,^{d} and John S. Lowengrub^{a,}^{*}

ude.icu.htam@knusnij (J.S. Sohn), wt.ude.utcn.htam@gnesthy (Y.-H. Tseng), ude.tii.htam@ils (S. Li), ed.nedserd-ut@tgiov.lexa (A. Voigt).

See other articles in PMC that cite the published article.

We develop and investigate numerically a thermodynamically consistent model of two-dimensional multicomponent vesicles in an incompressible viscous fluid. The model is derived using an energy variation approach that accounts for different lipid surface phases, the excess energy (line energy) associated with surface phase domain boundaries, bending energy, spontaneous curvature, local inextensibility and fluid flow via the Stokes equations. The equations are high-order (fourth order) nonlinear and nonlocal due to incompressibil-ity of the fluid and the local inextensibility of the vesicle membrane. To solve the equations numerically, we develop a nonstiff, pseudo-spectral boundary integral method that relies on an analysis of the equations at small scales. The algorithm is closely related to that developed very recently by Veerapaneni et al. [81] for homogeneous vesicles although we use a different and more efficient time stepping algorithm and a reformulation of the inextensibility equation. We present simulations of multicomponent vesicles in an initially quiescent fluid and investigate the effect of varying the average surface concentration of an initially unstable mixture of lipid phases. The phases then redistribute and alter the morphology of the vesicle and its dynamics. When an applied shear is introduced, an initially elliptical vesicle tank-treads and attains a steady shape and surface phase distribution. A sufficiently elongated vesicle tumbles and the presence of different surface phases with different bending stiffnesses and spontaneous curvatures yields a complex evolution of the vesicle morphology as the vesicle bends in regions where the bending stiffness and spontaneous curvature are small.

Biological membranes play an active and critical role in cell functions such as solute transport, cell locomotion, adhesion, signal transduction, etc. (e.g. see [1]). In addition, cells frequently use small membrane-bound carriers (vesicles) to transport proteins and molecules. The basic structural component of all biological membranes are lipid bilayers. Multicomponent vesicles are hollow, closed biomembranes with a lipid bilayer membrane containing different types of lipids and cholesterol. Vesicles serve as important, but simplified models of more complex cell-membranes [67]. Vesicle membranes are liquid-like, resist bending and are inextensible, e.g. [51,67]. Recent experiments on giant unilamellar vesicles demonstrate that there exists a rich variety of behavior of multicomponent vesicles. Spinodal decomposition into distinct surface domains (e.g. liquid-ordered, liquid-disordered) known as rafts (or domains), raft coarsening, viscous fingering, vesicle budding, fission and fusion are all observed with concomitant changes in membrane morphology. See, for example, [79,80,8,14,9,77, 18,22,84,68,5]. The resulting rafts may play an important role in regulating protein activity because of the existence of raft-resident proteins such as hemagglutinin which is an envelope protein of influenza [29].

While there have been many theoretical and numerical studies of homogeneous vesicles using discrete and continuum models (e.g. see the reviews [59,51,53,62,61]), the recent papers [16,17,15,12,11,10,35,36,42,6,81,72,63] (and the references therein), there are far fewer studies of inhomogeneous systems although there has been an increasing focus on the inhomo-geneous vesicles in the past 10 years.

Discrete approaches, such as Monte Carlo methods, dissipative particle dynamics and molecular dynamics, have been used to simulate the dynamics of phase separation and domain formation, vesicle fission and fusion (e.g. see [46,48,69,82,60,34,49,21,57,25,70,20]). However, the length and time scales achieved in such simulations are significantly limited by computational cost.

Continuum methods provide a good modeling alternative to reach larger length and time scales. Continuum models are also easier to analyze and parametrize. The continuum approach is based on the generalized bending energy proposed by Helfrich [28] supplemented by a line energy associated with surface phase domain boundaries (e.g. Lipowsky [52], Seifert [66] and Jülicher and Lipowsky [40,39]). Until recently, studies of multicomponent vesicles have been limited to equilibrium investigations (e.g. see [4,43,76,66,26,39,23,24,27,7,68,71]) or dynamical simulations limited to small deformations or special shapes (e.g. see [75,37,58,55,19]). Very recently, phase-field models developed for single-component vesicles (e.g. see [16,17,15,12,11,10,35,36,63]) have been extended to the multicomponent case [13,83,56]. In the phase-field approach, the vesicle membrane is given a finite thickness and the governing Helfrich equations are modeled with diffuse interface counterparts. While this approach is capable of simulating complex dynamics including vesicle adhesion, fusion and fission, the Helfrich equations are not solved directly and the local inextensibility constraints are typically replaced with global constraints or simplified approximations.

To our knowledge, there are no sharp interface models where the dynamics of multicomponent inextensible membranes in a viscous fluid are simulated directly. While the sharp interface approach cannot simulate vesicle fusion and fission without ad hoc cut-and-connect techniques, this method is valuable nevertheless because it provides highly accurate solutions that can be used to validate other more general approaches (e.g. phase-field) and the solutions provided by the sharp interface method are of intrinsic interest. In fact, there have been few studies using sharp interface models of inextensible homogeneous vesicles in a viscous fluid [44,85,74,81] because of the difficulties involved in solving the high-order nonlinear systems of equations. The stiffness introduced by bending forces posed a formidable challenge that was only recently overcome in two dimensions by Veerapaneni et al. [81] who developed a nonstiff algorithm by extracting dominant contributions to the equations at small spatial scales, following Hou et al. [30], Kropinski [45], Tornberg and Shelley [78] and Hou and Shi [32].

In this paper, we develop and investigate numerically a thermodynamically consistent model of two-dimensional multi-component vesicles in an incompressible viscous fluid. The model is derived using an energy variation approach that accounts for different surface phases, the excess energy (line energy) associated with surface phase domain boundaries, bending energy, spontaneous curvature, local inextensibility and fluid flow via the Stokes equations. The equations are high-order (fourth order) nonlinear and nonlocal due to incompressibility of the fluid and the local inextensibility of the vesicle membrane. To solve the equations numerically, we use a boundary integral representation of the fluid velocity and develop a nonstiff, pseudo-spectral boundary integral method to solve the coupled system that relies on an analysis of the equations at small scales. The algorithm is closely related to that developed very recently by Veerapaneni et al. [81] for homogeneous vesicles although we use a different and more efficient time stepping algorithm and a reformulation of the inextensibility equation. We present simulations of multicomponent vesicles in an initially quiescent fluid and investigate the effect of varying the average surface concentration of an initially unstable mixture of lipid phases. The phases then redistribute and alter the morphology of the vesicle and its dynamics. When an applied shear is introduced, an initially elliptical vesicle tank-treads and attains a steady shape and surface phase distribution. A sufficiently elongated vesicle tumbles and the presence of different surface phases with different bending stiffnesses and spontaneous curvatures yields a complex evolution of the vesicle morphology as the vesicle bends in regions where the bending stiffness and spontaneous curvature are small.

The rest of the paper is organized as follows: in Section 2, we present the boundary integral formulation of two-dimensional Stokes flow. In Section 3, the equations governing the dynamics of the inextensible, multicomponent vesicle and surface phases are derived. In Section 4, the numerical method is presented. Numerical results are given in Section 6. Conclusions are drawn and future work is discussed in Section 7.

Let **X** = **X**(*s*, *t*) = (*x*(*s*, *t*), *y*(*s*, *t*)) denote the position of a membrane bounding a vesicle as a function of arclength *s* and time *t*. The membrane evolution is given by

$$\frac{d\mathbf{X}}{\mathit{dt}}=V\mathbf{n}+T\mathbf{s},$$

(1)

where *V* is the outward normal velocity, **n** is the outward normal vector, *T* is the tangential velocity and **s** is the tangent vector. The velocity is determined via the solution of the fluid flow equations. We introduce the following nondimensional variables:

$${\mathbf{X}}^{\prime}=\frac{\mathbf{X}}{l},\phantom{\rule{thickmathspace}{0ex}}{t}^{\prime}=\frac{t}{\tau},\phantom{\rule{thickmathspace}{0ex}}{\mathbf{u}}^{\prime}=\frac{\tau}{l}\mathbf{u},$$

where *l* is a characteristic length of the membrane and *τ* = *v*_{2}*l*^{3}/ is the characteristic time, where *v*_{2} is the viscosity of the matrix fluid exterior to the vesicle and is a characteristic bending stiffness. We also introduce another parameter *Δ* known as the reduced surface area: *Δ* = *L/R* − 2*π*, where *L* is the surface area (total length in 2D) and *R* is the radius of the equivalent circle (circle with the same enclosed area as the vesicle). Note that *Δ* = 0 for a circle. This parameter has been used to classify vesicle shapes and in particular the transition to tumbling under an applied shear flow [41]. Hereafter, we drop the prime notation and present only the nondimensional equations.

Define the stress tensor

$${\mathbf{P}}_{i}=-{p}_{i}\mathbf{I}+{\nu}_{i}{\mathbf{D}}_{i},\phantom{\rule{1em}{0ex}}{\mathbf{D}}_{i}=(\nabla {\mathbf{u}}_{i}+\nabla {\mathbf{u}}_{i}^{T}),$$

(2)

where *p*_{i} is the pressure, with *i* = 1, 2 denoting the interior and exterior fluids, *v*_{i} are the nondimensional viscosities (e.g. *v*_{2} = 1 and *v*_{1} is the viscosity ratio), **D**_{i} is the deformation tensor and **u**_{i} is the velocity.

We will assume that the length and velocity scales are small so that Stokes flow governs the motion of the fluid

$$\nabla \cdot {\mathbf{P}}_{i}=0,\phantom{\rule{1em}{0ex}}\nabla \cdot {\mathbf{u}}_{i}=0,$$

(3)

where a velocity **u = u**_{∞} may be imposed in the far-field. Across the membrane, a stress jump condition holds [**P** · **n**]_{Σ} (**P** · **n**)_{1} − (**P** · **n**)_{2} is given by

$${[\mathbf{P}\cdot \mathbf{n}]}_{\Sigma}=-\mathbf{F},$$

(4)

where the stress **F** is obtained from energy variation and thermodynamic consistency (see Section 3.4). The velocity is assumed to be continuous across the membrane:

$${\left[\mathbf{u}\right]}_{\Sigma}=0.$$

(5)

In the remainder of the paper, we assume the viscosities of the fluids interior and exterior to the vesicle are matched: *v*_{1} = *v*_{2} = 1. An analysis of the effect of viscosity contrast is currently under study. Defining the stress componentwise to be **F**(*s, t*) = (*f*1,*f*2), the velocity **u** = (*u*1,*u*_{2}) of the vesicle membrane is given by the boundary integral formulae (e.g. [61])

$${u}_{1}=\frac{1}{4\pi}{\int}_{\Sigma}\left(-\mathrm{log}\phantom{\rule{thinmathspace}{0ex}}{\mathit{rf}}_{1}\left({s}^{\prime}\right)+{f}_{1}\frac{{r}_{1}^{2}}{{r}^{2}}+{f}_{2}\frac{{r}_{1}{r}_{2}}{{r}^{2}}\right){\mathit{ds}}^{\prime}+{u}_{\infty ,1},$$

(6)

$${u}_{2}=\frac{1}{4\pi}{\int}_{\Sigma}\left(-\mathrm{log}\phantom{\rule{thinmathspace}{0ex}}{\mathit{rf}}_{2}\left({s}^{\prime}\right)+{f}_{2}\frac{{r}_{2}^{2}}{{r}^{2}}+{f}_{1}\frac{{r}_{1}{r}_{2}}{{r}^{2}}\right){\mathit{ds}}^{\prime}+{u}_{\infty ,2},$$

(7)

where *r* = |**X**(*s*) − **X**(*s*′)|, *r*_{1} = *x*(*s*) − *x*(*s*′), *r*_{2} = *y*(*s*) − *y*(*s*′) and **u**_{∞} = (*u*_{∞,1}, *u*_{∞,2}) is the far-field velocity. Next, we need to determine a constitutive relation for the stress **F**.

Let *α* denote a parametrization of the membrane that bounds the vesicle (≤ *α* < 2*π*) and **s** = (*x _{α}*,

$${\theta}_{t}=-{V}_{s}+\kappa T,\phantom{\rule{1em}{0ex}}{s}_{\alpha t}=({T}_{s}+\kappa V){s}_{\alpha}$$

(8)

and

$${\kappa}_{t}=-({\partial}_{s}^{2}+{\kappa}^{2})V+T{\kappa}_{s}.$$

(9)

A thermodynamically consistent constitutive relation for the stress **F** is determined by requiring that the dynamics of the membrane decreases the total membrane free energy *E*^{M} = *E*^{b} + *E*^{T}. Here, *E*^{b} is the normal bending energy and *E*^{T} is the excess (line) energy associated with the surface phase (e.g. ordered/disordered) domain boundaries. For isothermal systems (which we consider here) a decreasing free energy is equivalent to an increasing entropy since the free energy is the total energy minus the temperature times the entropy [47]. The energies are given as follows.

Let *u* be the concentration of one of the lipid phases (e.g. disordered phase). Then the normal bending energy is given by

$${E}^{b}=\frac{1}{2}{\int}_{\Sigma}b\left(u\right){(\kappa -\stackrel{-}{\kappa}\left(u\right))}^{2}d\Sigma ,$$

(10)

, where *b* is the nondimensional normal bending stiffness, *κ* is the total curvature and is the spontaneous curvature. The line energy *E*^{T} is

$${E}^{T}=\frac{a}{\u220a}{\int}_{\Sigma}\left(f\left(u\right)+\frac{{\u220a}^{2}}{2}{\left|{\nabla}_{\Sigma}u\right|}^{2}\right)d\Sigma ,$$

(11)

, where *f* is a double-well potential, e.g $f\left(u\right)=\frac{1}{4}{u}^{2}{(1-u)}^{2}$, a is a nondimensional parameter (line tension scaled by a characteristic bending stiffness ), _{Σ} = (**I** - **nn**) is the surface gradient, and is a nondimensional small scaling parameter where the scaling is chosen such that as → 0, *E ^{T}* converges to a finite constant times

Next, we consider the variational derivatives of each energy. We do this by taking the time derivative which is equivalent to varying *u* and the membrane *Σ* simultaneously. Taking the time derivative of Eqs. (11), (10), using Eqs. (8), (9), integrating by parts and using periodicity (since *Σ* is assumed to be a closed curve) we get:

$${\stackrel{.}{E}}^{T}=\frac{a}{\u220a}{\int}_{\Sigma}({u}_{t}-{\mathit{Tu}}_{s})({f}^{\prime}\left(u\right)-{\u220a}^{2}{u}_{\mathit{ss}})\mathit{ds}+\frac{a}{\u220a}{\int}_{\Sigma}V\kappa (f\left(u\right)-\frac{{\u220a}^{2}}{2}{u}_{s}^{2})\mathit{ds}$$

(12)

and

$${\stackrel{.}{E}}^{b}={\int}_{\Sigma}({u}_{t}-{\mathit{Tu}}_{s})\left(\frac{{b}^{\prime}\left(u\right)}{2}{(\kappa -\stackrel{-}{\kappa})}^{2}-b\left(u\right)(\kappa -\stackrel{-}{\kappa})\frac{\partial \stackrel{-}{\kappa}}{\partial u}(u,s)\right)\mathit{ds}+{\int}_{\Sigma}V\left(-{\left(b\left(u\right)(\kappa -\stackrel{-}{\kappa})\right)}_{\mathit{ss}}-\frac{b\left(u\right)}{2}\kappa ({\kappa}^{2}-{\stackrel{-}{\kappa}}^{2})\right)\mathit{ds},$$

where the overdots denote time derivatives and the primes denote derivatives with respect to *u*. Combining Eqs. (12) and (13), and defining the variational derivatives,

$$\frac{\delta {E}^{M}}{\delta u}=\frac{a}{\u220a}({f}^{\prime}\left(u\right)-{\u220a}^{2}{u}_{\mathit{ss}})+\frac{{b}^{\prime}\left(u\right)}{2}{(\kappa -\stackrel{-}{\kappa})}^{2}-b\left(u\right)(\kappa -\stackrel{-}{\kappa})\frac{\partial \stackrel{-}{\kappa}}{\partial u},$$

(13)

$$\frac{\delta {E}^{M}}{\delta {\Sigma}^{n}}=-\left({\left(b\left(u\right)(\kappa -\stackrel{-}{\kappa})\right)}_{\mathit{ss}}+\frac{b\left(u\right)}{2}\kappa ({\kappa}^{2}-{\stackrel{-}{\kappa}}^{2})-\frac{a}{\u220a}\kappa \left(f\left(u\right)-\frac{{\u220a}^{2}}{2}{u}_{s}^{2}\right)\right),$$

(14)

we get

$${\stackrel{.}{E}}^{M}={\int}_{\Sigma}{u}_{t}\frac{\delta {E}^{M}}{\delta u}\mathit{ds}+{\int}_{\Sigma}V\frac{\delta {E}^{M}}{\delta {\Sigma}^{n}}\mathit{ds}-{\int}_{\Sigma}{\mathit{Tu}}_{s}\frac{\delta {E}^{M}}{\delta u}\mathit{ds}.$$

(15)

To make further progress, we need to consider mass conservation of the surface phases.

In the absence of reactions or adsorption/desorption, the total amount of lipid phases should remain unchanged. This means that the mass

$$M={\int}_{\Sigma}u\phantom{\rule{thinmathspace}{0ex}}\mathit{ds}$$

(16)

should be invariant in time, i.e. = 0. Defining *M*_{Γ} = *∫*_{Γ}*u* ds then * _{Γ}* =

$${\stackrel{.}{M}}_{\Gamma}={\int}_{\Gamma}({u}_{t}+u({T}_{s}+\kappa V))\mathit{ds}.$$

(17)

Thus, local mass conservation implies that

$${u}_{t}+u({T}_{s}+\kappa V)={J}_{s},$$

(18)

where the mass flux *J* is determined below to ensure energy dissipation (thermodynamic consistency). Using Eq. (18) in Eq. (15) gives

$${\stackrel{.}{E}}^{M}={\int}_{\Sigma}({J}_{s}-u({T}_{s}+\kappa V))\frac{\delta {E}^{M}}{\delta u}\mathit{ds}+{\int}_{\Sigma}V\frac{\delta {E}^{M}}{\delta {\Sigma}^{n}}\mathit{ds}-{\int}_{\Sigma}{\mathit{Tu}}_{s}\frac{\delta {E}^{M}}{\delta u}\mathit{ds}.$$

(19)

.

Now we account for the local inextensibility of the membrane. That is, the membrane is not allowed to stretch. Consequently, the arclength parametrization should be time independent, *s*_{α}(*α, t*) = *s*_{α}(*α*,0), which implies

$${T}_{s}+\kappa V=0.$$

(20)

Note that an alternate form of Eq. (20) is **s** · u_{s} = 0. Following previous work, we pose this constraint via a Lagrange multiplier Λ^{LL}(*α, t*). Hence, we append the membrane energy considered earlier with a null Lagrangian *E ^{LL}*, where

$${E}^{\mathit{LL}}={\int}_{\Sigma}{\Lambda}^{\mathit{LL}}({s}_{\alpha}(\alpha ,t)-{s}_{\alpha}(\alpha ,0))d\alpha .$$

(21)

Note that *Λ ^{LL}* can also be interpreted as a tensile stress introduced to enforce inextensibility. Taking the time derivative, we get

$${\stackrel{.}{E}}^{\mathit{LL}}={\int}_{\Sigma}{\Lambda}_{t}^{\mathit{LL}}({s}_{\alpha}-{s}_{\alpha}(\alpha ,0))d\alpha +{\int}_{\Sigma}{\Lambda}^{\mathit{LL}}({T}_{s}+\kappa V)\mathit{ds}.$$

(22)

Putting everything together, and assuming that *Λ ^{LL}* is chosen such that

$${\stackrel{.}{E}}^{M}={\stackrel{.}{E}}^{T}+{\stackrel{.}{E}}^{b}-{\stackrel{.}{E}}^{\mathit{LL}}={\int}_{\Sigma}{J}_{s}\frac{\delta {E}^{M}}{\delta u}\mathit{ds}+{\int}_{\Sigma}V\left(\frac{\delta {E}^{M}}{\delta {\Sigma}^{n}}-{\Lambda}^{\mathit{LL}}\kappa \right)\mathit{ds}+{\int}_{\Sigma}T\left({\Lambda}_{s}^{\mathit{LL}}-{u}_{s}\frac{\delta {E}^{M}}{\delta u}\right)\mathit{ds}.$$

(23)

Note that local inextensibility makes Eq. (18) become

$${u}_{t}={J}_{s}.$$

(24)

.

We are now in a position to pose thermodynamically consistent constitutive relations for the flux *J* and the stress **F**. Assuming that each term in Eq. (23) is dissipative, we pose the following consitutive relations for the *J* and **F**:

$$J=\frac{1}{\mathit{Pe}}{\left(\frac{\delta {E}^{M}}{\delta u}\right)}_{s},$$

(25)

where *Pe* is a nondimensional Peclet number, and

$$\mathbf{F}=-\left(\frac{\delta {E}^{M}}{\delta {\Sigma}^{n}}-{\Lambda}^{\mathit{LL}}\kappa \right)\mathbf{n}-\left({\Lambda}_{s}^{\mathit{LL}}-{u}_{s}\frac{\delta {E}^{M}}{\delta u}\right)\mathbf{s}.$$

(26)

A calculation shows that the stress **F** can be written as

$$\mathbf{F}={\left(b\left(u\right)(\kappa -\stackrel{-}{\kappa})\mathbf{n}\right)}_{ss}-\frac{1}{2}{\left(b\left(u\right)(\kappa -\stackrel{-}{\kappa})(3\kappa +\stackrel{-}{\kappa})\mathbf{s}\right)}_{s}+{\left(\frac{a}{\u220a}\left(f\left(u\right)-\frac{{\u220a}^{2}}{2}{u}_{u}^{2}\right)\mathbf{s}\right)}_{s}-{\left({\Lambda}^{\mathit{LL}}\mathbf{s}\right)}_{s}.$$

(27)

This gives the jump in normal stress:

$${\left[\mathbf{Pn}\right]}_{\Sigma}=-{\left(b\left(u\right)(\kappa -\stackrel{-}{\kappa})\mathbf{n}\right)}_{\mathit{ss}}+\frac{1}{2}{\left(b\left(u\right)(\kappa -\stackrel{-}{\kappa})(3\kappa +\stackrel{-}{\kappa})\mathbf{s}\right)}_{s}-{\left(\frac{a}{\u220a}\left(f\left(u\right)-\frac{{\u220a}^{2}}{2}{u}_{s}^{2}\right)\mathbf{s}\right)}_{s}+{\left({\Lambda}^{\mathit{LL}}\mathbf{s}\right)}_{s},$$

(28)

which shows that the bending and surface phase evolution induce zero total force and torque on the system. With these choices, we obtain the dissipation formula:

$${\stackrel{.}{E}}^{M}=-\frac{1}{\mathit{Pe}}{\int}_{\Sigma}{\left({\left(\frac{\delta {E}^{M}}{\delta u}\right)}_{s}\right)}^{2}\mathit{ds}-\frac{1}{2}{\int}_{\Sigma}\mathbf{D}:\mathbf{D}\mathit{ds},$$

(29)

where $J=\frac{1}{\mathit{Pe}}{\left(\frac{\delta {E}^{M}}{\delta u}\right)}_{s}$ and **D** = (**u** + **u**^{T}).

To close the system of equations governing the membrane evolution, the local Lagrange multiplier *Λ ^{LL}* is chosen such that

$${(\mathbf{u}\cdot \mathbf{s})}_{s}=-\kappa (\mathbf{u}\cdot \mathbf{n}),$$

(30)

, where **u** · **s** and **u** · **n** are the tangential and normal components of the Stokes velocity on the membrane and thus these velocity components are nonlocal linear functionals of *Λ ^{LL}* through the stress

On the membrane, we decompose the Stokes velocity **u** as

$$\mathbf{u}=\mathbf{v}+\stackrel{~}{\mathbf{u}},$$

(31)

where **v** is a length conserving velocity field which is derived using a global Lagrange multiplier *Λ ^{GL}* that preserves the overall length of the membrane. The velocity

To begin, we decompose the velocity **v** as

$$\mathbf{v}={\mathbf{v}}^{u}+{\Lambda}^{\mathit{GL}}\mathbf{w},$$

(32)

, where $\mathbf{v}={\mathbf{v}}_{1}^{u}{\mid}_{\Sigma}={\mathbf{v}}_{2}^{u}{\mid}_{\Sigma}$ is the solution of the unconstrained Stokes equation:

$$\begin{array}{cc}\hfill & \nabla \cdot {\mathbf{P}}_{i}^{u}=0,\phantom{\rule{1em}{0ex}}{\mathbf{P}}_{i}^{u}=-{p}_{i}^{u}\mathbf{I}+{\mathbf{D}}_{i}^{u},\phantom{\rule{1em}{0ex}}{\mathbf{D}}_{i}^{u}=\left(\nabla {\mathbf{v}}_{i}^{u}+\nabla {{\mathbf{v}}_{i}^{u}}^{T}\right),\hfill \\ \hfill & \nabla \cdot {\mathbf{v}}_{i}^{u}=0,\hfill \\ \hfill & {\left[{\mathbf{P}}^{u}\mathbf{n}\right]}_{\Sigma}=\frac{\delta {E}^{M}}{\delta {\Sigma}^{n}}\mathbf{n}-{u}_{s}\frac{\delta {E}^{M}}{\delta u}\mathbf{s},\phantom{\rule{1em}{0ex}}{\left[{\mathbf{v}}^{u}\right]}_{\Sigma}=0\hfill \end{array}$$

(33)

and **w** = **w**_{1}|_{Σ} = **w**_{2}|_{Σ} satisfies

$$\begin{array}{cc}\hfill & \nabla \cdot {\mathbf{P}}_{i}^{w}=0,\phantom{\rule{1em}{0ex}}{\mathbf{P}}_{i}^{w}=-{p}_{i}^{w}\mathbf{I}+{\mathbf{D}}_{i}^{w},\phantom{\rule{1em}{0ex}}{\mathbf{D}}_{i}^{w}=(\nabla {\mathbf{w}}_{i}+\nabla {\mathbf{w}}_{i}^{T}),\hfill \\ \hfill & \nabla \cdot {\mathbf{w}}_{i}=0,\hfill \\ \hfill & {\left[{\mathbf{P}}^{w}\mathbf{n}\right]}_{\Sigma}=\kappa \mathbf{n},\phantom{\rule{1em}{0ex}}{\left[\mathbf{w}\right]}_{\Sigma}=0.\hfill \end{array}$$

(34)

Further, *Λ ^{GL}* is chosen such that

$${\int}_{\Sigma}\kappa \mathbf{v}\cdot \mathbf{n}\mathit{ds}=0,$$

(35)

so that the velocity field **v** is length conserving. Rewriting Eq. (35) using the decomposition of **v** from Eq. (32), we get

$${\int}_{\Sigma}\kappa ({\mathbf{v}}^{u}+{\Lambda}^{\mathit{GL}}\mathbf{w})\cdot \mathbf{n}\mathit{ds}=0,$$

(36)

which implies that

$${\Lambda}^{\mathit{GL}}=-\frac{{\int}_{\Sigma}\kappa {\mathbf{v}}^{u}\cdot \mathbf{n}\mathit{ds}}{{\int}_{\Sigma}\kappa \mathbf{w}\cdot \mathbf{n}\mathit{ds}}.$$

(37)

To ensure local inextensibility, we use the correction **ū** which satisfies:

$$\begin{array}{cc}\hfill & \nabla \cdot {\stackrel{~}{\mathbf{P}}}_{i}=0,\phantom{\rule{1em}{0ex}}{\stackrel{~}{\mathbf{P}}}_{i}=-{\stackrel{~}{p}}_{i}\mathbf{I}+{\stackrel{~}{\mathbf{D}}}_{i},\phantom{\rule{1em}{0ex}}{\stackrel{~}{\mathbf{D}}}_{i}=(\nabla {\stackrel{~}{\mathbf{u}}}_{i}+\nabla {\stackrel{~}{\mathbf{u}}}_{i}^{T}),\hfill \\ \hfill & \nabla \cdot {\stackrel{~}{\mathbf{u}}}_{i}=0,\hfill \end{array}$$

(38)

$${\left[\stackrel{~}{\mathbf{P}}\mathbf{n}\right]}_{\Sigma}={\left(\stackrel{~}{\Lambda}\mathbf{s}\right)}_{s},{\left[\stackrel{~}{\mathbf{u}}\right]}_{\Sigma}=0,$$

(39)

where (*α, t*) = *Λ ^{LL}*(

$${(\stackrel{~}{\mathbf{u}}\cdot \mathbf{s})}_{s}+\kappa \stackrel{~}{\mathbf{u}}\cdot \mathbf{n}=-({(\mathbf{v}\cdot \mathbf{s})}_{s}+\kappa \mathbf{v}\cdot \mathbf{n}).$$

(40)

Thus, Eqs. (38)-(40) comprise the linear nonlocal system for (*α, t*) that now needs to be solved (e.g. by iteration). This is done using the small scale decomposition described in the next section, following [81]. While this system has a similar form as the original system, the reformulated system has the advantage that the only external forcing is due to the local stretching present in the length conserving velocity **v** through the right hand side of Eq. (40). This tends to reduce the number of iterations required to solve the system.

We remark that there exists at least one other inextensible velocity field. This is seen as follows. Since the mean of *κ***v** · **n** is equal to zero (by choice of *Λ ^{GL}*), it is possible to solve (

$$\mathbf{v}={\mathbf{v}}^{u}+{\Lambda}^{\mathit{GL}}\mathbf{W}+\zeta \mathbf{s},$$

(41)

where *ζ* is chosen so that the velocity field **v** is actually locally inextensible:

$${\zeta}_{s}=-{\left(\left({\mathbf{v}}^{u}+{\Lambda}^{\mathit{GL}}\mathbf{w}\right)\cdot \mathbf{s}\right)}_{s}-\kappa \left({\mathbf{v}}^{u}+{\Lambda}^{\mathit{GL}}\mathbf{w}\right)\cdot \mathbf{n}.$$

(42)

Integrating this equation in *s* gives *ζ*. Although **v** is an inextensible velocity field, it is not generated by a single Lagrange multiplier. This is seen as follows. Given the velocity *ζ***s**, one may solve Eqs. (6) and (7) for the corresponding stress ${\mathbf{f}}^{\zeta}=({f}_{1}^{\zeta},{f}_{2}^{\zeta})$; this is a system of first kind Fredholm equations.^{1} Doing this, one observes that **f**^{ζ} ≠ (*β***s**)_{s} for some *β* which would be required in order to have a single Lagrange multiplier description. Nevertheless the velocity in Eq. (41) may be used to reformulate the original system for *Λ ^{LL}* instead of the velocity in Eq. (32) but ultimately both approaches require solving Eq. (40).

The small scale decomposition(SSD) extracts the dominant part of the equations at small spatial scales [30,31]. In the context of Stokes flows, this was first performed in [45] and later in [32] for immersed interface method and in [81] in the context of the dynamics of inextensible, homogeneous vesicles. Ref. [81] is most relevant to our work here although unlike [81] we follow [30] and use the SSD to develop an explicit, nonstiff time integration algorithm.

In Eqs. (6) and (7), the only singularity in the integrand comes from the logarithmic kernel. This can be analyzed as follows:

$$\mathrm{log}\phantom{\rule{thinmathspace}{0ex}}\left|\mathbf{r}\right|=\mathrm{log}\phantom{\rule{thinmathspace}{0ex}}2\left|\mathrm{sin}\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|+\mathrm{log}\frac{\left|\mathbf{r}\right|}{2\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\left|\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|},$$

(43)

where the second term in Eq. (43) is smooth. For *α* ≈ *α*′, we have,

$$\left|\mathbf{r}\right|=\sqrt{(\mathit{xs}\left(\alpha \right)-x{\left(s\left({\alpha}^{\prime}\right)\right)}^{2}+y\left(s\left(\alpha \right)\right)-y{\left(s\left({\alpha}^{\prime}\right)\right)}^{2}}={s}_{\alpha}|\alpha -{\alpha}^{\prime}|\sqrt{1+O(\alpha -{\alpha}^{\prime})}={s}_{\alpha}|\alpha -{\alpha}^{\prime}|(1+O(\alpha -{\alpha}^{\prime}))$$

(44)

where *O*(*α* – *α*′) denotes a smooth function that vanishes as *α*′ → *α* [30]. Thus, $\left|\mathbf{r}\right|\u2215\left(2\left|\mathrm{sin}\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|\right)$ is a smooth function.

Analyzing Eqs. (6) and (7) at small spatial scales, we obtain [45,32,81]:

$$\mathbf{u}(\alpha ,t)~\frac{{\stackrel{-}{s}}_{\alpha}}{4\pi}{\int}_{0}^{2\pi}\mathrm{log}\left(2\left|\mathrm{sin}\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|\right)\mathbf{F}({\alpha}^{\prime},t)d{\alpha}^{\prime},$$

(45)

where the notation *f* ~ *g* means that *f* – *g* is a smoothing operator, e.g. an integral operator with a smooth kernel. In addition, _{α} = *s _{α}* is constant in space and time because of inextensibility.

Let g(α) be a smooth function and define the integral operator

$$\mathcal{L}\left[g\right]\left(\alpha \right)=\frac{1}{\pi}{\int}_{0}^{2\pi}\mathit{log}\left(2\left|\mathit{sin}\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|\right)\mathit{g}\left({\alpha}^{\prime}\right)d{\alpha}^{\prime}.$$

(46)

Then the symbol of is = −1/|k|, where k is the Fourier wavenumber. Therefore, [g](k) = −ĝ(k)/|k|, for k ≠ 0. For k = 0, the symbol is equal to 0. The proof of this lemma can be found in [38,81*], for example.*

Next, consider the small scale decomposition of the normal *V* = **u**(*α, t*) · **n**(*α, t*) tangential *T* = **u**(*α, t*) · **s**(*α, t*) velocities that arise from the Stokes flow:

$$V(\alpha ,t)~\frac{{\stackrel{-}{s}}_{\alpha}}{4\pi}{\int}_{0}^{2\pi}\mathrm{log}\left(2\left|\mathrm{sin}\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|\right)\mathbf{F}({\alpha}^{\prime},t)\cdot \mathbf{n}({\alpha}^{\prime},t)d{\alpha}^{\prime},$$

(47)

$$T(\alpha ,t)~\frac{{\stackrel{-}{s}}_{\alpha}}{4\pi}{\int}_{0}^{2\pi}\mathrm{log}\left(2\left|\mathrm{sin}\left(\frac{\alpha -{\alpha}^{\prime}}{2}\right)\right|\right)\mathbf{F}({\alpha}^{\prime},t)\cdot \mathbf{s}({\alpha}^{\prime},t)d{\alpha}^{\prime}.$$

(48)

Further, from Eq. (27)

$$\mathbf{F}\cdot \mathbf{n}~b{\kappa}_{\mathit{ss}}+{\Lambda}^{\mathit{LL}}\kappa =\frac{b}{{\stackrel{-}{s}}_{\alpha}^{3}}{\theta}_{\alpha \alpha \alpha}+\frac{{\Lambda}^{\mathit{LL}}}{{\stackrel{-}{s}}_{\alpha}}{\theta}_{\alpha},$$

(49)

$$\mathbf{F}\cdot \mathbf{s}~-{\Lambda}_{s}^{\mathit{LL}}=\frac{1}{{\stackrel{-}{s}}_{\alpha}}+{\Lambda}_{\alpha}^{\mathit{LL}}.$$

(50)

Plugging Eqs. (49) and (50) into Eqs. (47) and (48), we obtain

$$V(\alpha ,t)~-\frac{b}{4{\stackrel{-}{s}}_{\alpha}^{2}}\mathcal{L}\left[{\theta}_{\alpha \alpha \alpha}\right](\alpha ,t)-\frac{{\theta}_{\alpha}}{4}\mathcal{L}\left[{\Lambda}^{\mathit{LL}}\right],$$

(51)

$$T(\alpha ,t)~\frac{1}{4}\mathcal{L}\left[{\Lambda}_{\alpha}^{\mathit{LL}}\right].$$

(52)

Next, we use the SSD to design an efficient, nonstiff explicit time integration scheme and an efficient preconditioner for the local Lagrange multiplier.

Following [30], we use the SSD to rewrite Eq. (8) as

$${\theta}_{t}=\frac{\stackrel{-}{b}}{4{\stackrel{-}{s}}_{\alpha}^{3}}{\partial}_{\alpha}\mathcal{L}\left[{\theta}_{\alpha \alpha \alpha}\right](\alpha ,t)+{N}_{1}(\alpha ,t),$$

(53)

where we have absorbed the *Λ ^{LL}* contributions into

$${N}_{1}(\alpha ,t)=-\frac{{V}_{\alpha}}{{\stackrel{-}{s}}_{\alpha}}+\kappa T(\alpha ,t)-\frac{\stackrel{-}{b}}{4{\stackrel{-}{s}}_{\alpha}^{3}}{\partial}_{\alpha}\mathcal{L}\left[{\theta}_{\alpha \alpha \alpha}\right](\alpha ,t).$$

(54)

In the equations above, is constant and we find that taking = max *b*(*u*) works well in practice. Taking the Fourier transform of Eq. (53), we get

$${\widehat{\theta}}_{t}(k,t)=-\frac{\stackrel{-}{b}}{4{\stackrel{-}{s}}_{\alpha}^{3}}{\left|k\right|}^{3}\widehat{\theta}(k,t)+{\widehat{N}}_{1}(k,t).$$

(55)

Following [30], we discretize Eq. (55) using the second-order accurate linear propagator method:

$${\widehat{\theta}}^{n+1}\left(k\right)={e}_{k}({t}_{n},{t}_{n+1}){\widehat{\theta}}^{n}\left(k\right)+\frac{\mathrm{\Delta}t}{2}\left(3{e}_{k}({t}_{n},{t}_{n+1}){\widehat{N}}_{1}^{n}\left(k\right)-{e}_{k}({t}_{n-1},{t}_{n+1}){\widehat{N}}_{1}^{n-1}\right),$$

(56)

where

$${e}_{k}({t}_{1},{t}_{2})=\mathrm{exp}\left(-\frac{\stackrel{-}{b}}{4}{\left|k\right|}^{3}{\int}_{{t}_{1}}^{{t}_{2}}\frac{\mathit{dt}\prime}{{\stackrel{-}{s}}_{\alpha}^{3}\left({t}^{\prime}\right)}\right).$$

(57)

We use the trapezoidal rule to approximate

$${\int}_{{t}_{n}}^{{t}_{n+1}}\frac{{\mathit{dt}}^{\prime}}{{\stackrel{-}{s}}_{\alpha}{\left({t}^{\prime}\right)}^{3}}\approx \frac{\mathrm{\Delta}t}{2}\left(\frac{1}{{\left({\stackrel{-}{s}}_{\alpha}^{n}\right)}^{3}}+\frac{1}{{\left({\stackrel{-}{s}}_{\alpha}^{n+1}\right)}^{3}}\right),$$

(58)

$${\int}_{{t}_{n-1}}^{{t}_{n+1}}\frac{{\mathit{dt}}^{\prime}}{{\stackrel{-}{s}}_{\alpha}{\left({t}^{\prime}\right)}^{3}}\approx \mathrm{\Delta}t\left(\frac{1}{2{\left({\stackrel{-}{s}}_{\alpha}^{n-1}\right)}^{3}}+\frac{1}{{\left({\stackrel{-}{s}}_{\alpha}^{n}\right)}^{3}}+\frac{1}{2{\left({\stackrel{-}{s}}_{\alpha}^{n+1}\right)}^{3}}\right).$$

(59)

The arclength * _{α}*(

$${\stackrel{-}{s}}_{\alpha}^{n+1}={\stackrel{-}{s}}_{\alpha}^{n}+\frac{\mathrm{\Delta}t}{2}\left(3{M}^{n}-{M}^{n-1}\right),$$

(60)

where

$$M=\frac{1}{2\pi}{\int}_{0}^{2\pi}{\theta}_{{\alpha}^{\prime}}V({\alpha}^{\prime},t)d{\alpha}^{\prime}.$$

(61)

To reconstruct the membrane interface (*x*(*α*, *t*_{n+1}), *y*(*α*, *t*_{n+1})) from the updated *θ*^{n+1}(*α*) and ${\stackrel{-}{S}}_{\alpha}^{n+1}$, we first update a reference point (*x*(0, *t*_{n+1}), *y*(0, *t*_{n+1})). This is done by using a second-order explicit Adams–Bashforth method. Once we update the reference point, we obtain the configuration of the interface from the *θ*^{n+1}(*α*) and ${\stackrel{-}{S}}_{\alpha}^{n+1}$ by [30]

$$x(\alpha ,{t}_{n+1})=x(0,{t}_{n+1})+{\stackrel{-}{s}}_{\alpha}^{n+1}\left({\int}_{0}^{\alpha}\mathrm{cos}\left({\theta}^{n+1}\left({\alpha}^{\prime}\right)\right)d{\alpha}^{\prime}-\frac{\alpha}{2\pi}{\int}_{0}^{2\pi}\mathrm{cos}\left({\theta}^{n+1}\left({\alpha}^{\prime}\right)\right)d{\alpha}^{\prime}\right),$$

$$y(\alpha ,{t}_{n+1})=y(0,{t}_{n+1})+{\stackrel{-}{s}}_{\alpha}^{n+1}\left({\int}_{0}^{\alpha}\mathrm{sin}\left({\theta}^{n+1}\left({\alpha}^{\prime}\right)\right)d{\alpha}^{\prime}-\frac{\alpha}{2\pi}{\int}_{0}^{2\pi}\mathrm{sin}\left({\theta}^{n+1}\left({\alpha}^{\prime}\right)\right)d{\alpha}^{\prime}\right),$$

where the integration is performed using the discrete Fourier transform.

Note that the linear propagator and Adams–Bashforth methods require two previous time steps. At the first time step to compute *θ*^{1}, we replace the second-order linear propagator with a first-order linear propagator of a similar form. For ${\stackrel{-}{S}}_{\alpha}^{1}$, we use the explicit Euler method.

To calculate the Stokes boundary integrals in Eqs. (6) and (7) we use the trapezoidal rule for the integrals with smooth integrands. To calculate the integral with the logarithmic integrand, we follow the approach in [38] and use the SSD and the lemma to obtain a spectrally accurate discretization using the discrete Fourier transform. Finally, following [30], all derivatives in the membrane, concentration and local Lagrange multiplier equations are performed using the discrete Fourier transform. A high-order (25th order) Fourier smoothing is used to control aliasing errors.

We solve the integro-differential inextensibility Eq. (40) for using GMRES with a preconditioner derived from the SSD (see also [81]). From the SSD, we obtain

$${(\stackrel{-}{\mathbf{u}}\cdot \mathbf{s})}_{s}+\kappa \stackrel{-}{\mathbf{u}}\cdot \mathbf{n}~\frac{1}{4{\stackrel{-}{s}}_{\alpha}}{\partial}_{\alpha}\mathcal{L}\left[{\widehat{\Lambda}}_{\alpha}\right]-\frac{{\theta}_{\alpha}^{2}}{4{\stackrel{-}{s}}_{\alpha}}\mathcal{L}\left[\widehat{\Lambda}\right].$$

(62)

Replacing *θ _{α}* by a characteristic value

$$A\left[\stackrel{~}{\Lambda}\right]=\frac{1}{4{\stackrel{-}{s}}_{\alpha}}{\partial}_{\alpha}\mathcal{L}\left[{\stackrel{~}{\Lambda}}_{\alpha}\right]-\frac{{\stackrel{-}{\theta}}_{\alpha}^{2}}{4{\stackrel{-}{s}}_{\alpha}}\mathcal{L}\left[\stackrel{~}{\Lambda}\right],$$

(63)

which can be inverted in Fourier space using the symbol:

$${\widehat{\mathcal{A}}}^{-1}=-\frac{1}{\left|k\right|}{\left(\frac{1}{4{\stackrel{-}{s}}_{\alpha}}-\frac{{\stackrel{-}{\theta}}_{\alpha}^{2}}{4{\stackrel{-}{s}}_{\alpha}}\frac{1}{{\left|k\right|}^{2}}\right)}^{-1}.$$

(64)

Interestingly, using *A*^{−1} as the preconditioner, we found that the best results are obtained by taking *θ _{α}* = 0.

Using the small scale decomposition for the diffusion flux *J*

$$J~-\frac{a\u220a}{\mathit{Pe}{\stackrel{-}{s}}_{\alpha}^{3}}{u}_{\alpha \alpha \alpha},$$

(65)

we may rewrite Eq. (24) as:

$${u}_{t}=-\frac{a\u220a}{\mathit{Pe}{\stackrel{-}{s}}_{\alpha}^{4}}{u}_{\alpha \alpha \alpha \alpha}+{N}_{2}(\alpha ,t),$$

(66)

where

$${N}_{2}(\alpha ,t)=\frac{1}{\mathit{Pe}}\left({\left(\frac{\delta {E}^{M}}{\delta u}\right)}_{\mathit{ss}}+\frac{a\u220a}{{\stackrel{-}{s}}_{\alpha}^{4}}{u}_{\alpha \alpha \alpha \alpha}\right).$$

(67)

Taking the Fourier transform of Eq. (66) to get

$${\widehat{u}}_{t}(k,t)=-\frac{a\u220a}{\mathit{Pe}{\stackrel{-}{s}}_{\alpha}^{4}}{\left|k\right|}^{4}\widehat{u}(k,t)+{\widehat{N}}_{2}(k,t),$$

(68)

we can use an analogous linear propagator method to perform the time discretization:

$${\widehat{u}}^{n+1}\left(k\right)={\eta}_{k}({t}_{n},{t}_{n+1}){\widehat{u}}^{n}\left(k\right)+\frac{\mathrm{\Delta}t}{2}\left(3{\eta}_{k}({t}_{n},{t}_{n+1}){\widehat{N}}_{2}^{n}\left(k\right)-{\eta}_{k}({t}_{n-1},{t}_{n+1}){\widehat{N}}_{2}^{n-1}\right),$$

(69)

where

$${\eta}_{k}({t}_{1},{t}_{2})=\mathrm{exp}\left(-\frac{a\u220a}{\mathit{Pe}}{\left|k\right|}^{4}{\int}_{{t}_{1}}^{{t}_{2}}\frac{{\mathit{dt}}^{\prime}}{{\stackrel{-}{s}}_{\alpha}^{4}\left({t}^{\prime}\right)}\right).$$

(70)

We also use the trapezoidal rule to approximate ${\int}_{{t}_{1}}^{{t}_{2}}\frac{{\mathit{dt}}^{\prime}}{{\stackrel{-}{S}}_{\alpha}^{4}\left({t}^{\prime}\right)}$ as in Eqs. (58) and (59). The first time step *u*^{1} is obtained using a first-order linear propagator method.

In this section, we summarize the governing equations and numerical procedure. For the fluid, the equations are:

$$\mathrm{\nabla}\cdot {\mathbf{P}}_{i}=0,\phantom{\rule{1em}{0ex}}\mathrm{\nabla}\cdot {\mathbf{u}}_{i}=0,$$

(71)

$${[\mathbf{P}\cdot \mathbf{n}]}_{\Sigma}=\left(\frac{\delta {E}^{b}}{\delta {\Sigma}^{n}}-{\Lambda}^{\mathit{LL}}\kappa \right)\mathbf{n}+\left({\Lambda}_{s}^{\mathit{LL}}-{u}_{s}\frac{\delta {E}^{M}}{\delta u}\right)\mathbf{s},$$

(72)

$${\left[\mathbf{u}\right]}_{\Sigma}=0,$$

(73)

$$\mathbf{u}\to {\mathbf{u}}_{\infty}\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}\text{in the far-field},$$

(74)

where ${\mathbf{P}}_{i}=-{p}_{i}\mathbf{I}+\nu (\nabla {\mathbf{u}}_{i}+\nabla {\mathbf{u}}_{i}^{T})$ and *i* = 1, 2 denotes the interior and exterior fluids. The function *u* is the concentration of one of the lipid phases, and *Λ ^{LL}*(

$${u}_{t}=\frac{1}{\mathit{Pe}}{\left(\frac{\delta {E}^{M}}{\delta u}\right)}_{\mathit{ss}},$$

(75)

where *Pe* is a dimensionless parameter, and the tangent angle *θ* to the interface evolves by

$${\theta}_{t}=-{V}_{s}+\kappa T,$$

(76)

where *V* and *T* are the normal and tangential velocities given by

$$V=\mathbf{u}(\mathbf{X}(\alpha ,t),t)\cdot \mathbf{n},$$

(77)

$$T=\mathbf{u}(\mathbf{X}(\alpha ,t),t)\cdot \mathbf{s}.$$

(78)

The arclength variation *s _{α}* evolves according to

$${\partial}_{t}{s}_{\alpha}=({T}_{s}+\kappa V){s}_{\alpha}.$$

(79)

Even though *T _{s}* +

The outline of the nonstiff numerical algorithm used to solve the system is as follows:

- For a given location of the membrane
*Σ*(*t*) at time_{n}*t*=*t*_{n}, use the iterative linear solver GMRES with the preconditioner described earlier to solve the Stokes equations to obtain**u**^{n}, the inextensible velocity field. - Update the membrane tangent angle
*θ*^{n+1}, the arclength variation ${S}_{\alpha}^{n+1}$ and the concentration*u*^{n+1}using the nonstiff time stepping algorithm outlined previously. - Reconstruct the new interface positions (
*x*(*α*,*t*_{n+1}),*y*(*α*,*t*_{n+1})) from ${S}_{\alpha}^{n+1}$,*θ*^{n+1}and a reference point (*x*_{0}(*t*),*y*_{0}(*t*)). - Goto step 1 and repeat.

In this section, we validate the numerical algorithm and investigate the influence of the line tension, bending coefficient, spontaneous curvature and external flow on two-dimensional membrane dynamics. We also vary the relative ratio of the surface lipid phase components. Unless otherwise specified, we use the following parameter and initial configuration setting: *a* = 100, spontaneous curvature (*u*) = 5(1 − *u*) + 0.1*u* and bending stiffness *b*(*u*) = (1 − *u*) + 0.5*u*. Also, we use a double well potential function *f*(*u*) = 0.25*u*^{2}(1 − *u*)^{2}. The initial condition is an ellipse given by:

$$\begin{array}{cc}\hfill x(\alpha ,0)& =\mathrm{cos}\left(\alpha \right),\hfill \\ \hfill y(\alpha ,0)& =0.3\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\left(\alpha \right),\hfill \end{array}$$

(80)

where (0 ≤ *α* < 2 *π*). With this initial condition, the reduced surface area *Δ* = *L/R* − 2*π* = 1.7244. The initial concentration *u* is assumed to be a mixture of both phases:

$$u(\alpha ,0)=\stackrel{-}{u}+\delta (\mathrm{cos}\left(\alpha \right)+\mathrm{cos}\left(3\alpha \right)+\mathrm{cos}\left(4\alpha \right)),$$

(81)

where the average concentration *ū* is varied between 0 and 1 to account for different ratios of the surface phases and the perturbation parameter *δ* = 0.001.

We demonstrate the effectiveness of our proposed semi-implicit method by comparing with a second-order accurate, explicit Adams–Bashforth (AB2) method used for all the equations. We consider the evolution from an unstable 50-50 mixture of the two lipid phases by taking *ū* = 0.5 with dimensionless parameters = 0.1 and *Pe* = 1. In Fig. 1, we present the Fourier transform |(*k, T*)| of the radius $r(\alpha ,T)=\sqrt{x{(\alpha ,T)}^{2}+y{(\alpha ,T)}^{2}}$ using different time steps and mesh sizes for the two approaches (semi-implicit method, AB2). Only the even Fourier modes are shown at a final time *T*.

Comparison of the Fourier spectra |(*k, T*) for an explicit Adams–Bashforth method (a and c) and our semi-implicit non-stiff algorithm (b and d) for *N* grid points, time step Δ*t* and final time *T* as indicated.

Using AB2 with *N* = 256 and Δ*t* = 1 × 10^{−10}, the simulation is only able to compute up to *T* = 2.6 × 10^{−8}. Around this time, there is significant and rapid growth of the high wavenumber amplitudes which then destabilize the simulation. See Fig. 1(a) (dashed). If the time step is reduced by a factor of two to Δ*t* = 0.5 × 10^{−10} (solid), simulation may run significantly longer indicating this is just below the critical stability threshold for the time step. In Fig. 1(b) the same simulation is performed using our semi-implicit, nonstiff algorithm. In contrast, a much larger time step may be used (Δ*t* = 1.0 × 10^{−4} (solid)) and the simulation may proceed indefinitely without losing stability (the time *T* = 4 × 10^{−3} is shown). As seen in the figure the high wavenumbers do not grow and that the solution is insensitive to the time step; the dashed corresponds to Δ*t* = 0.5 × 10^{−4}.

In Fig. 1(c), results are presented for AB2 using *N* = 512. In this case, the time step needs to be reduced by about a factor of 10 from that used with *N* = 256, which is consistent with the third-order time step constraint predicted by the small scale analysis. In the figure, the dashed curve corresponds to Δ*t* = 1 × 10^{−11}. Observe there is significant growth of high wavenumbers at *T* = 5.8 × 10^{−10} which destabilizes the simulation. Using a time step Δ*t* = 0.5 × 10^{−11}, the simulation may be run significantly longer. In contrast, performing the same simulation with our nonstiff algorithm [Fig. 1(d)] indicates that we can actually use the same time step for both *N* = 256 and *N* = 512 without loss of stability. At worst, we may have to reduce the time step by a factor of two for each factor of two decrease in the mesh size. Thus, it is clear that the use of AB2 is not practical and AB2 is significantly outperformed by our semi-implicit nonstiff method.

We test the numerical accuracy of our scheme by considering the evolution from an unstable 50–50 mixture of the two lipid phases by taking *ū* = 0.5 with varying *Pe* and . Consider first the resolution in time. Using *N* = 512 grid points, we perform simulations using time steps Δ*t* = 4 × 10^{−5}, 2 × 10^{−5}, 1 × 10^{−5}, 5 × 10^{−6} and 2.5 × 10^{−6}. The error in local arclength is max_{α} |*s _{α}*(

Convergence rate at time *t* = 0.02 for a 50–50 mixture of lipid phases with *ū* = 0.5, *N* = 512 grid points with *Pe* = 1 and varying (0.2, 0.1 and 0.05).

We note that the scheme is stable and still quite accurate using much larger time steps. However, as seen from Table 1, the convergence rate deteriorates as the time step increases although the actual errors are small. The use of a higher order accurate time stepping algorithm could be helpful in this regard (e.g. see [33]). This is currently under study. We thus chose to use a small time step to ensure the temporal errors are very small and that we are in the second-order convergence regime. We also investigated the rate of convergence with different choices of the Peclet number *Pe* = 1/ and *Pe* = (results not shown). Second-order convergence is also observed although the use of smaller values of require a smaller time step when *Pe* = for stability. This is not surprising since the initial phase decomposition occurs more rapidly when decreases with *Pe* = compared to that when *Pe* = 1.

To test the effect of spatial resolution, we fixed the time step Δ*t* = 2 × 10^{−5} and varied the spatial grid size using *N* = 256,512 and 1024. In each case, the area is preserved to more than 6 digits, indicating that the error is limited by the time step. The membrane positions for these resolutions are virtually indistinguishable.

Next, we continue the evolution described in Section 6.2 and study the surface phase separation of a vesicle containing a 50-50 mixture of lipid phases (*ū* = 0.5) for *Pe* = 1 and = 0.1. We use the time step Δ*t* = 2 × 10^{−5} and *N* = 512 grid points. The results are shown in Fig. 2(a)–(c). In (a), the total energy *E ^{M}* is shown together with insets of the corresponding vesicle morphologies and surface concentration

Multicomponent vesicle evolution with a 50–50 mixture of lipid phases: *ū* = 0.5. (a) The evolution of the total energy *E*^{M}. Insets: vesicle morphologies and surface phase concentration *u* at indicated times (color online: blue and red correspond **...**

The corresponding normal velocities are shown in Fig. 2(c). At early times the largest (in the magnitude) velocities occur near the vesicle tips and tend to round off the tip. Away from the tip, the negative normal velocities create the negatively curved regions. Since the velocity tends to zero, this confirms that the evolution tends to a steady state. In this simulation, the area is preserved to six digits and the maximum error in local arclength (as defined above) is 6 × 10^{−8}.

To see how the dynamics is affected by the Peclet number, we perform the same simulation as in Fig. 2 except using *Pe* = 1/ and *Pe* = . The results are shown in Fig. 3(a) for the concentration and (b) for the energies. As seen from the figure, the initial stages of phase separation occur at different rates, depending on the Peclet number. Once the phases separate, the evolution of the different cases is similar.

Multicomponent vesicle evolution with a 50–50 mixture of lipid phases: *ū* = 0.5 with fixed = 0.1 and varying *Pe* = (dash), *1* (dash dot) and 1/ (solid). (a) The surface phase concentration *u* at indicated times. **...**

Now we fix the dimensionless parameters to be the following: = 0.1 and *Pe* = 1. In Fig. 4 we consider a 70–30 mixture of lipid phases by taking *ū* = 0.7. As in the 50–50 case, the *u* ≈ 1 and *u* ≈ 0 phases initially separate at the top/bottom and tip regions of the vesicle, respectively. The energies *E ^{M}*,

Multicomponent vesicle evolution with a 70–30 mixture of lipid phases: *ū* = 0.7 with = 0.1 and *Pe* = 1. (a) The evolution of the total energy *E*^{M}. Insets: vesicle morphologies and surface phase concentration *u* at indicated times **...**

In Fig. 5, we increase the amount of the *u* = 1 component and consider a 80–20 mixture by taking *ū* = 0.8. At early times, the phases tend to separate more slowly compared to the 50–50 and 70–30 cases. This is because the driving force for phase separation is lower in this case. In fact, the initial concentration *ū* lies just outside the spinodal region of *f*(0.25 ≤ *u* ≤ 0.75) where negative diffusion (*ff*′ < 0) occurs. Here, phase separation is driven by the bending energy at early times. As can be seen from the Fig. 5(b), the early dynamics tend to reduce the bending energy *E ^{b}* at the expense of a slight increase in the line energy

Multicomponent vesicle evolution with a 80–20 mixture of lipid phases: *ū* = 0.8 with = 0.1 and *Pe* = 1. (a) The evolution of the total energy *E*^{M}. Insets: vesicle morphologies and surface phase concentration *u* at indicated times **...**

When the amount of the *u* = 1 phase is increased to yield a 85–15 mixture by taking *ū* = 0.85, the phases rearrange slightly but as seen in Fig. 6, there is no dramatic phase separation as the vesicle tends toward a discocyte shape. This is because the initial concentration lies far enough outside the spinodal region of *f*(0.25 ≤ *u* ≤ 0.75) that no decomposition occurs. In this simulation (Δ*t* = 2 × 10^{−5}, *N* = 512), the area is preserved to six digits of accuracy and the maximum error in local arclength is 3 × 10^{−8}.

Multicomponent vesicle evolution with a 85–15 mixture of lipid phases: *ū* = 0.85 with = 0.1 and *Pe* = 1. (a) The evolution of the total energy *E*^{M}. Insets: vesicle morphologies and surface phase concentration *u* at indicated times **...**

Next, in Fig. 7, we decrease the amount of the *u* = 1 phase and consider a 30–70 mixture by taking *ū* = 0.3. As in the previous cases, phase separation occurs at the vesicle top/bottom and tips, accompanied by a decrease in the energies *E ^{M}*,

Multicomponent vesicle evolution with a 30–70 mixture of lipid phases: *ū* = 0.3 with = 0.1 and *Pe* = 1. (a) The evolution of the total energy *E*^{M}. Insets: vesicle morphologies and surface phase concentration *u* at indicated times **...**

In Fig. 8 we repeat the simulation shown in Fig. 7 but with a smaller value of ; we take = 0.05. Since the initial concentration condition *ū* = 0.3 is far from equilibrium, using a smaller value of increases the initial magnitude of the line energy (via the coefficient *a*/) *E ^{T}* and correspondingly the total energy

We next consider the evolution of a multicomponent vesicle under an applied shear flow. We consider a 30–70 lipid phase mixture (*ū* = 0.3). The initial vesicle is defined by

$$x(\alpha ,0)=0.1\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\left(\alpha \right),$$

(82)

$$y(\alpha ,0)=\mathrm{sin}\left(\alpha \right),$$

(83)

where *α* [0,2*π*) and the simple shear flow **u**_{∞} = (25y,0) is used. All other parameters are the same as in the previous simulations.

With this initial condition, the reduced surface area is *Δ* = 6.5682 which is significantly larger than that used in the previous simulations. Experimentally, there is a critical value *Δ** such that for *Δ* ≤ *Δ** a vesicle will tank-tread in a shear flow (e.g. reach a steady shape) but for *Δ* > *Δ** the vesicle will tumble (e.g. evolve dynamically and the shape will remain unsteady). See [41]. As far as we are aware, the value *Δ** has not been determined for two-dimensional vesicles or multicomponent vesicles. Numerically, we found that for *Δ* = 6.5682 the vesicle tumbles (as seen below) while for *Δ* = 1.7244 (initial membrane given in Eq. (80) the vesicle tank-treads). A more precise determination of *Δ** and its dependence on heterogeneity and the associated parameters is currently under study.

In Fig. 9 we consider the evolution of the initially elliptical vesicle with *Δ* = 1.7244. For reference, the point (*x*(*α* = 0, *t*), *y*(*α* = 0, *t*)) is marked by ‘*’. At early times, the vesicle rotates, the tips round off and the phases separate on the flat regions of the vesicle, eventually accumulating on one of the vesicle sides, as seen in (a). Rather quickly, the vesicle reaches a nearly steady shape and the reference point travels around the interface indicating the persistence of tangential motion (tank-treading). The *u* = 1 phase, however, remains in a fixed spatial location and causes the membrane to bend slightly, and even develop a small negative curvature, due to its smaller bending stiffness and spontaneous curvature. This is precisely due to the interaction between the shape variation and the phase decomposition (e.g. presence of curvature terms in Eq. (13) and the presence of the *u*-dependent terms in the stress boundary condition in Eq. (72)). In Fig. 9(b) the concentration *u* is plotted versus *α* at different times throughout the simulation. Because the spatial location of the *u* = 1 domain is fixed, the concentration tends to a traveling wave (when plotted against *α* as the vesicle approaches a steady shape. In Fig. 10, the corresponding normal (a) and tangential (b) velocities are shown. The normal velocity tends to a small, apparently periodic perturbation about zero, indicating that the shape is nearly steady. The tangential velocity, on the other hand, approaches a small perturbation about negative constant. A careful examination of the shape (not shown) indicates that it is not quite steady and appears to vary periodically in time where the perturbation about the mean shape is very small, corresponding to the small perturbations in the normal and tangential velocities.

Multicomponent vesicle tank-treading under an applied shear with a 30–70 mixture of lipid phases: *ū* = 0.3 with *Δ* = 1.7244, = 0.1 and *Pe* = 1. (a) The evolution of the vesicle at the indicated times (color online: blue **...**

Next, we consider the evolution of the initially elliptical vesicle with *Δ* = 6.5682. In Fig. 11 the evolution of the resulting vesicle morphology (a) is shown together with the concentration *u* (b). As above, the reference point (*x*(*α* = 0, *t*), *y*(*α* = 0, *t*)) is shown as the point‘*’. At early times, the *u* = 1 phase decomposes on the flat regions of the vesicle as it begins to rotate. The *u* = 1 phase ultimately accumulates on one of the vesicle sides and as the vesicle continues to rotate, the decreased bending coefficient and spontaneous curvature of the *u* = 1 phase cause the vesicle to bend as it rotates (see times *t* = 0.4 and *0.7*) with the *u* = 1 region acquiring a negative curvature while the *u* = 0 phase retains a positive curvature. By time *t* = 0.4, the vesicle becomes nearly horizontal and acquires a stomatocyte-like shape analogous to that seen in Fig. 7. At this point the vesicle is again bent by the shear flow and the process repeats itself apparently indefinitely. As this process continues, the distance between the negatively curved *u* = 1 region and the positively curved *u* = 0 region on the opposite side of the vesicle decreases but does not apparently vanish in finite time as seen in Fig. 12. Indeed, the neck radius may even stabilize. This figure also shows the results of a resolution study in Δ*t* (Δ*t* is halved with *N*= 512) and in *N* (*N* is doubled to 1024). The plot shows that the all the resolutions are in agreement and thus the tumbling vesicle is well resolved using Δ*t* = 1 × 10^{−5} and *N* = 512. Further, the maximum error in local arclength is 1 × 10^{−4}.

A multicomponent vesicle tumbling under an applied shear with a 30–70 mixture of lipid phases: *ū* = 0.3 with *Δ* = 6.5682, = 0.1 and *Pe* = 1. (a) The evolution of the vesicle at the indicated times (color online: blue and **...**

The corresponding normal velocities are shown in Fig. 13(a). At early times the largest (in the magnitude) velocities occur near the vesicle tips and tends to round them off. Observe that the velocity roughly has a periodic or quasi-periodic structure that emerges in time (cf. compare times *t* ≈ 1.47, *t* ≈ 2.75 and *t* ≈ 6). Further evidence of this is shown in Fig. 13(b) where the maximal distance of the vesicle membrane from the origin is plotted as a function of time. As seen in Fig. 13(c) the temporal spacing between the consecutive valleys in the maximal distance decreases and may tend to a constant value as time increases. This is currently under study.

The morphologies observed in Fig. 11(a) critically depend on the heterogeneity of the vesicle. If one instead considers a homogeneous vesicle with *u* = 0 and all other parameters as before, then the tumbling morphologies are quite different. This is shown in Fig. 14. Observe that the homogeneous vesicle resists being bent by the flow and bends only slightly in the opposite direction as compared to the multicomponent case. This is consistent with the large positive spontaneous curvature (0) = 5 for the *u* = 0 phase. Interestingly, we found that by taking either different bending coefficients or different spontaneous curvatures for the *u* = 0 and *u* = 1 phases in a multicomponent vesicle also results vesicle morphologies analogous to those observed in Fig. 11.

We have developed and investigated a thermodynamically consistent model of two-dimensional multicomponent vesicles in an incompressible viscous fluid. The model is derived using an energy variation approach that accounts for variation among the surface phases and the associated the excess energy (line energy) between juxtaposed surface phases, bending energy, spontaneous curvature, local inextensibility and fluid flow via the Stokes equations. In particular, we derived the generalized Laplace–Young boundary condition for the normal stress jump across the vesicle membrane. The equations are high-order (fourth order) nonlinear and nonlocal due to incompressibility of the fluid and the local inextensibility of the vesicle membrane.

To solve the equations numerically, we developed a nonstiff, pseudo-spectral boundary integral method that relies on an analysis of the equations at small scales. The small scale decomposition is used to develop efficient, explicit time integration schemes for both the evolution of the membrane and the surface phases. In addition, the equation for the Lagrange multiplier introduced to enforce local inextensibility is reformulated and the small scale decomposition is used to motivate an efficient preconditioner. The algorithm is closely related to that developed very recently by Veerapaneni et al. [81] for homogeneous vesicles although we use a different and more efficient time stepping algorithm [30] and a reformulation of the Lagrange multiplier equation. The algorithm presented here can be improved in several ways. A higher order accurate time stepping algorithm could be employed (e.g. see [33]) to improve accuracy in time and a fast summation algorithm (e.g. see [45,81]) could be used to speed up the evaluation of the boundary integrals. These are currently under study.

We presented simulations of multicomponent vesicles in an initially quiescent fluid and in an applied shear flow. In the former, we investigated the effect of varying the average surface concentration of the lipid phases (e.g. ordered/disordered) where the initial concentration is a perturbation about the average and the initial vesicle is an ellipse. In addition, the *u* = 1 phase has smaller spontaneous curvature and bending stiffness coefficient. At early times, the phases decompose in a similar way with the *u* = 1 phase preferentially accumulating in flat regions of the vesicle. The phases then redistribute and alter the morphology of the vesicle and its dynamics according to the average concentration. When an applied shear is introduced, an initially elliptical vesicle tank-treads and attains a steady shape and surface phase distribution. A sufficiently elongated vesicle tumbles and the presence of different surface phases with different bending stiffnesses and spontaneous curvatures yields a complex evolution of the vesicle morphology as the vesicle bends in regions where the bending stiffness and spontaneous curvature are small.

The model presented here can be extended straightforwardly to systems with more surface components and more realistic thermodynamics, e.g. [64,73,9] and to three-dimensional vesicles. We are currently extending the numerical algorithm to axisymmetric three-dimensional vesicles. The presence of membrane proteins within the phases can also have an important influence on the vesicle shape and surface phase distribution, e.g. [39,54], and can be modeled by incorporating the concentration(s) of membrane proteins in the phase-field model described here via an additional protein energy and introducing coupling between the protein and the bending rigidity and spontaneous curvature, e.g. see [50,2,3]. Although the boundary integral model breaks down if vesicles fission or coalesce, and ad hoc interface surgery needs to be performed to transition through such events, the boundary integral approach is attractive because it provides highly accurate solutions that can be used to validate other more general approaches (e.g. level-set, phase-field) and the solutions provided by the boundary integral method are of intrinsic interest.

J.L. thanks Michael Shelley for valuable discussions. A.V. acknowledges support from DFG SFB 609 and VO-899/6-1. J.L. acknowledges support from the National Science Foundation Division of Mathematical Sciences (DMS) and from the National Institutes of Health through grant P50GM76516 for a Center of Excellence in Systems Biology at the University of California, Irvine. S.L. acknowledges support from NSF grant DMS-0914923. Finally, the authors thank the reviewers for their insightful comments that strengthened the paper.

^{1}Actually this system has a nullspace consisting of functions of the form **f** = *c***n** where *c* is a constant. A unique solution perpendicular to this nullspace may be found.

1. Alberts B, Bray D, Lewis J, Raff M, Roberts K, Watson JD. Molecular Biology of the Cell. Garland; New York: 1994.

2. Allain J-M, Ben Amar M. Biphasic vesicle: instability induced by adsorption of proteins. Phys. Rev. A. 2004;337:531–545.

3. Allain J-M, Ben Amar M. Budding and fission of a multiphase vesicle. Eur. Phys. J. E. 2006;20:409–420. [PubMed]

4. Andelman D, Kawakatsu T, Kawasaki K. Equilibrium shape of two-component unilamellar membranes and vesicles. Europhys. Lett. 1992;19(1):57–62.

5. Ayuyan AG, Cohen FS. Raft composition at physiological temperature and pH in the absence of detergents. Biophys. J. 2008;94:2654–2666. [PubMed]

6. Barrett JW, Garcke H, Nürnberg R. Parametric approximation of willmore flow and related geometric evolution equations. SIAM J. Sci. Comput. 2008;31(1):225–253.

7. Baumgart T, Das S, Webb WW, Jenkins JT. Membrane elasticity in giant vesicles with fluid phase coexistence. Biophys. J. 2005;89:1067–1080. [PubMed]

8. Baumgart T, Hess ST, Webb WW. Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension. Nature. 2003;435:821–824. [PubMed]

9. Beattie ME, Veatch SL, Stottrip BL, Keller SL. Sterol structure determines miscibility versus melting transitions in lipid vesicles. Biophys. J. 2005;89:1760–1768. [PubMed]

10. Biben T, Kassner K, Misbah C. Phase-field approach to three-dimensional vesicle dynamics. Phys. Rev. E. 2005;72:049121. [PubMed]

11. Biben T, Misbah C. Tumbling of vesicles under shear flow within an advected-field approach. Phys. Rev. E. 2003;67:031908. [PubMed]

12. Campelo F, Hernández-Machado A. Dynamic model and stationary shapes of fluid vesicles. Eur. Phys. J. E. 2006;20:37–45. [PubMed]

13. Campelo F, Hernández-Machado A. Model for curvature-driven pearling instability in membranes. Phys. Rev. Lett. 2007;99:088101. [PubMed]

14. D'Onofrio TG, Hatzor A, Counterman AE, Heetderks JJ, Sandel MJ, Weiss PS. Controlling and measuring the interdependence of local properties in biomembranes. Langmuir. 2003;19:1618–1623.

15. Du Q, Li ML, Liu C. Analysis of a phase field Navier–Stokes vesicle–fluid interaction model. Discr. Contin. Dyn. Syst. B. 2007;8(3):539–556.

16. Du Q, Liu C, Wang X. A phase field approach in the numerical study of the elastic bending energy for vesicle membranes. J. Comput. Phys. 2004;198(2):450–468.

17. Du Q, Liu C, Wang X. Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions. J. Comput. Phys. 2006;212:757.

18. Dutta S, Ray DS. Magnetic field induced pattern formation in reactive membranes. Phys. Rev. E. 2007;75:016205. [PubMed]

19. Funkhouser CM, Solis FJ, Thornton K. Coupled composition-deformation phase-field method for multicomponent lipid membranes. Phys. Rev. E. 2007;76:011912. [PubMed]

20. Gao LH, Lipowsky R, Shillcock J. Tension-induced vesicle fusion: pathways and pore dynamics. Soft Mat. 2008;4(6):1208–1214.

21. Gao LH, Shillcock J, Lipowsky R. Improved dissipative particle dynamics simulations of lipid bilayers. J. Chem. Phys. 2007;126(1):015101. [PubMed]

22. García-Sáez AJ, Chiantia S, Schwille P. Effect of line tension on the lateral organization of lipid membranes. J. Biol. Chem. 2007;282:33537–33544. [PubMed]

23. Góźdź WT, Gompper G. Shapes and shape transformations of two-component membranes of complex topology. Phys. Rev. E. 1999;59(4):4305–4316.

24. Góźdź WT, Gompper G. Shape transformations of two-component membranes under weak tension. Europhys. Lett. 2001;55(4):587–593.

25. Grafmüller A, Shillcock J, Lipowsky R. Pathway of membrane fusion with two tension-dependent energy barriers. Phys. Rev. Lett. 2007;98(21):218101. [PubMed]

26. Harden JL, MacKintosh FC. Shape transformations of domains in mixed-fluid films and bilayer membranes. Europhys. Lett. 1994;28(7):495–500.

27. Harden JL, MacKintosh FC, Olmsted PD. Budding and domain shape transformations in mixed lipid films and bilayer membranes. Phys. Rev. E. 2005;72:011903. [PubMed]

28. Helfrich W. Elastic properties of lipid bilayers-theory and possible experiments. Z. Nat. C. 1973;28:693–703. [PubMed]

29. Hess ST, Kumar M, Verma A, Farrington J, Kenworthy A, Zimmerberg J. Quantitative electron microscopy and fluorescence spectroscopy of the membrane distribution of influenza hemagglutinin. J. Cell Biol. 2005;169:965–976. [PMC free article] [PubMed]

30. Hou T, Lowengrub J, Shelley M. Removing the stiffness from interfacial flows with surface tension. J. Comput. Phys. 1994;114:312.

31. Hou T, Lowengrub J, Shelley M. Boundary integral methods for multicomponent fluids and multiphase materials. J. Comput. Phys. 2001;169:302–362.

32. Hou T, Shi Z. Removing the stiffness of elastic force from the immersed boundary method for the 2d Stokes equations. J. Comput. Phys. 2007;114:312.

33. Hou TY, Lowengrub JS, Shelley MJ. The long-time motion of vortex sheets with surface tension. Phys. Fluids. 1997;7:1933–1954.

34. Huang KC, Mukhopadhyay R, Wingreen NS. A curvature-mediated mechanism for localization of lipids to bacterial poles. PLOS Comput. Biol. 2006;2(11):1357–1364. [PubMed]

35. Jamet D, Misbah C. Toward a thermodynamically consistent picture of the phase-field model of vesicles: local membrane incompressibility. Phys. Rev. E. 2007;76:051907. [PubMed]

36. Jamet D, Misbah C. Toward a thermodynamically consistent picture of the phase-field model of vesicles: curvature energy. Phys. Rev. E. 2008;78:031902. [PubMed]

37. Jiang Y, Lookman T, Saxena A. Phase separation and shape deformation of two-phase membranes. Phys. Rev. E. 2000;61:R57–R61. [PubMed]

38. Jou HJ, Leo PH, Lowengrub J. Microstructural evolution in inhomogeneous elastic media. J. Comput. Phys. 1997;131:109–148.

39. Jülicher F, Lipowsky R. Shape transformations of vesicles with intramembrane domains. Phys. Rev. E. 1997;53:2670–2683. [PubMed]

40. Jülicher F, Lipwosky R. Domain-induced budding of vesicles. Phys. Rev. Lett. 1993;70:2964–2967. [PubMed]

41. Kantsler V, Steinberg V. Transition to tumbling and two regimes of tumbling motion of a vesicle in shear flow. Phys. Rev. Lett. 2006;96:036001. [PubMed]

42. Kaoui B, Ristow GH, Cantat I, Misbah C, Zimmermann W. Lateral migration of a two-dimensional vesicle in unbounded poiseuille flow. Phys. Rev. E. 2008;77:021903. [PubMed]

43. Kawakatsu T, Kawasaki K, Andelman D, Taniguchi T. Phase-transitions and shapes of 2-component membranes and vesicles 2. Strong segregation limit. J. Phys. II. 1993;3(7):971–997.

44. Kraus M, Wintz W, Seifert U, Lipowsky R. Fluid vesicles in shear flow. Phys. Rev. Lett. 1996;77:3685–3688. [PubMed]

45. Kropinski MCA. Numerical methods for multiple inviscid interfaces in creeping flows. J. Comput. Phys. 2002;180:1–24.

46. Kumar PBS, Gompper G, Lipowsky R. Budding dynamics of multicomponent membranes. Phys. Rev. Lett. 2001

47. Landau L. Statistical Physics. Butterworth-Heinemann; Oxford: 1984.

48. Laradji M, Kumar PBS. Dynamics of domain growth in self-assembled fluid vesicles. Phys. Rev. Lett. 2004;93(19):198105. [PubMed]

49. Laradji M, Kumar PBS. Domain growth, budding, and fission in phase-separating self-assembled fluid bilayers. J. Chem. Phys. 2005;123 [PubMed]

50. Leibler S. Curvature instability in membranes. J. Phys. 1986;47(3):507–516.

51. Lipowsky R. The conformation of membranes. Nature. 1991;349:475–481. [PubMed]

52. Lipowsky R. Budding of membranes induced by intramembrane domains. J. Phys. II France. 1992;2:1825–1840.

53. Lipowsky R, Sackman E, editors. Structure and Dynamics of Membranes – From Cells to Vesicles. 1A and B. Elsevier; Amsterdam: 1995.

54. Long MS, Cans A-S, Keating CD. Budding and asymmetric protein microcompartmentation in giant vesicles containing two aqueous phases. J. Am. Chem. Soc. 2008;130:756–762. [PubMed]

55. Lowengrub J, Xu J-J, Voigt A. Surface phase separation and flow in a simple model of multicomponent drops and vesicles. Fluid Dyn. Mater. Proc. 2007;3(1):1–19.

56. Lowengrub JS, Rätz A, Voigt A. Phase-field Modeling of the Dynamics of Multicomponent Vesicles: Spinodal Decomposition, Coarsening, Budding and Fission. Phys. Rev. E. 2009;79:031926. [PMC free article] [PubMed]

57. Markvoort AJ, Smeijers AF, Pieterse K, van Santen RA, Hilbers PAJ. Lipid-based mechanisms for vesicle fission. J. Phys. Chem. B. 2007;111:5719–5725. [PubMed]

58. McWhirter JL, Ayton GS, Voth GA. Coupling field theory with mesoscopic dynamical simulations of multicomponent lipid bilayers. Biophys. J. 2004;87:3242–3263. [PubMed]

59. Nelson DR, Piran T, Weinberg S, editors. Statistical Mechanics of Membranes and Surfaces. World Scientific; Singapore: 2004.

60. Noguchi H, Gompper G. Meshless membrane model based on the moving least-squares method. Phys. Rev. E. 2006;73:021903. [PubMed]

61. Pozrikidis C. Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press; Cambridge: 1992.

62. Pozrikidis C. Interfacial dynamics for stokes flow. J. Comput. Phys. 2001;169:250.

63. Du Q, Zhang J. Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations. SIAM J. Sci. Comput. 2008;30(3):1634–1657.

64. Ruiz-Argüello MB, Goñi FM, Alonso A. Vesicle membrane fusion induced by the concerted activities of sphingomyelinase and phospholipase c*. J. Biol. Chem. 1998;273(36):22977–22982. [PubMed]

65. Saad Y, Schultz MH. Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 1986;7:856–869.

66. Seifert U. Curvature-induced lateral phase segregation in two-component vesicles. Phys. Rev. Lett. 1993;70:1335–1338. [PubMed]

67. Seifert U. Configurations of fluid membranes and vesicles. Adv. Phys. 1997;46:13–137.

68. Semrau S, Idema T, Holtzer L, Schmidt T, Storm C. Accurate determination of elastic parameters for multicomponent membranes. Phys. Rev. Lett. 2008;100:088101. [PubMed]

69. Shillcock J, Lipowsky R. Tension-induced fusion of bilayer membranes and vesicles. Nat. Mater. 2005;4(3):225–228. [PubMed]

70. Smith KA, Upsal WE. Shear-driven release of a bud from a multicomponent vesicle. J. Chem. Phys. 2007;126:075102. [PubMed]

71. Solis FJ, Funkhouser CM, Thornton K. Conditions for overall planarity in membranes: applications to multicomponent membranes with lamellar morphology. Eur. Phys. Lett. 2008;82:38001.

72. Song P, Hu D, Zhang P. Numerical simulation of fluid membranes in two-dimensional space. Commun. Comput. Phys. 2008;3(4):794–821.

73. Stottrup BL, Stevens DS, Keller SL. Miscibility of ternary mixtures of phospholipids and cholesterol in monolayers, and application to bilayer systems. Biophys. J. 2005;88:269–276. [PubMed]

74. Sukumaran S, Seifert U. Influence of shear flow on vesicles near a wall: a numerical study. Phys. Rev. E. 2001;64:011916. [PubMed]

75. Taniguchi T. Shape deformation and phase separation dynamics of two-component vesicles. Phys. Rev. Lett. 1996;76:4444–4447. [PubMed]

76. Taniguchi T, Kawasaki K, Andelman D, Kawakatsu T. Phase transitions and shapes of two component membranes and vesicles: ii. Weak segregation limit. J. Phys. II. 1994;4(8):1333–1362.

77. Tian AW, Johnson C, Wang W, Baumgart T. Line tension at fluid membrane domain boundaries measured by micropipette aspiration. Phys. Rev. Lett. 2007;98:208102. [PubMed]

78. Tornberg A-K, Shelley MJ. Simulating the dynamics and interactions of flexible fibers in stokes flows. J. Comput. Phys. 2004;196:8–40.

79. Tsafrir I, Sagi D, Arzi T, Guedeau-Boudeville MA, Frette V, Kandel D, Stavans J. Pearling instabilities of membrane tubes with anchored polymers. Phys. Rev. Lett. 2001;86(6):1138–1141. [PubMed]

80. Veatch SL, Keller SL. Separation of liquid phases in giant vesicles of ternary mixtures of phospholipids and cholesterol. Biophys. J. 2003;85 [PubMed]

81. Veerapaneni SK, Gueyffier D, Zorin D, Biros G. A Boundary Integral Method for Simulating the Dynamics of Inextensible Vesicles Suspended in a Viscous Fluid in 2d. preprint.

82. Wallace E, Hooper NM, Olmsted PD. The kinetics of phase separation in asymmetric membranes. Biophys. J. 2005;88:4072–4083. [PubMed]

83. Wang XQ, Du Q. Modelling and simulations of multi-component lipid membranes and open membranes via diffuse interface approaches. J. Math. Biol. 2008;56:347–371. [PubMed]

84. Yuan J, Hira SM, Strouse GF, Hirst LS. Lipid bilayer discs and banded tubules: photoinduced lipid sorting in ternary mixtures. J. Am. Chem. Soc. 2008;130:2067–2072. [PubMed]

85. Zhou H, Pozrikidis C. Deformation of liquid capsules with incompressible interfaces in simple shear flow. J. Fluid Mech. 1995;283:175–200.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |