|Home | About | Journals | Submit | Contact Us | Français|
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.  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 ). 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 . 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 .
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  supplemented by a line energy associated with surface phase domain boundaries (e.g. Lipowsky , Seifert  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.  who developed a nonstiff algorithm by extracting dominant contributions to the equations at small spatial scales, following Hou et al. , Kropinski , Tornberg and Shelley  and Hou and Shi .
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.  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
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:
where l is a characteristic length of the membrane and τ = v2l3/ is the characteristic time, where v2 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 . Hereafter, we drop the prime notation and present only the nondimensional equations.
Define the stress tensor
where pi is the pressure, with i = 1, 2 denoting the interior and exterior fluids, vi are the nondimensional viscosities (e.g. v2 = 1 and v1 is the viscosity ratio), Di is the deformation tensor and ui is the velocity.
We will assume that the length and velocity scales are small so that Stokes flow governs the motion of the fluid
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
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:
In the remainder of the paper, we assume the viscosities of the fluids interior and exterior to the vesicle are matched: v1 = v2 = 1. An analysis of the effect of viscosity contrast is currently under study. Defining the stress componentwise to be F(s, t) = (f1,f2), the velocity u = (u1,u2) of the vesicle membrane is given by the boundary integral formulae (e.g. )
where r = |X(s) − X(s′)|, r1 = x(s) − x(s′), r2 = 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α,yα)/sα where . In the above, and in the remainder of the paper, subscripts denote partial derivatives. The outward normal vector is then n = (yα, − xα)/sα. Introducing the tangent angle θ = tan−1yα/xα, the tangent and normal vectors are s = (cosθ, sinθ) and n = (sinθ, − cosθ). Further, the total curvature is κ = ηα/sα = θs, where s is arclength. The Frenet formulas give ns = κs and ss = − κn. Using these relations, we may derive the following dynamical equations for the tangent angle θ, the arclength variation sα and the total curvature κ:
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 EM = Eb + ET. Here, Eb is the normal bending energy and ET 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 . 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
, where b is the nondimensional normal bending stiffness, κ is the total curvature and is the spontaneous curvature. The line energy ET is
, where f is a double-well potential, e.g , 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, ET converges to a finite constant times a. In 3D, this constant is proportional to the perimeter of the surface domain boundaries.
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:
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
should be invariant in time, i.e. = 0. Defining MΓ = ∫Γu ds then Γ = J|Γ, where Γ is any portion of Σ,J is the mass flux and J|Γ is the net mass flux into Γ. Taking the time derivative, we get
Thus, local mass conservation implies that
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
Note that an alternate form of Eq. (20) is s · us = 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 ELL, where
Note that ΛLL can also be interpreted as a tensile stress introduced to enforce inextensibility. Taking the time derivative, we get
Putting everything together, and assuming that ΛLL is chosen such that Ts = −κV, gives
Note that local inextensibility makes Eq. (18) become
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:
where Pe is a nondimensional Peclet number, and
A calculation shows that the stress F can be written as
This gives the jump in normal stress:
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:
where and D = (u + uT).
To close the system of equations governing the membrane evolution, the local Lagrange multiplier ΛLL is chosen such that
, 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 F. This system, which on the discrete level is nonsymmetric, can be solved using an iterative solver such as GMRES . Each iteration requires a solution of the Stokes equations. An efficient preconditioner was recently developed by Veerapaneni et al.  using an analysis of the Stokes integral operators at small spatial scales (small scale decomposition). While we take that approach here (see Section 4), we further reformulate Eq. (30) to obtain an even more efficient solver.
On the membrane, we decompose the Stokes velocity u as
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 v is easy to compute (no iteration is required) but suffers from the fact that it is not locally inextensible. Thus, the velocity ū is introduced to correct this. Numerically, we find that ū tends to be small which accelerates the convergence of the iterative solver. We now show how v is obtained.
To begin, we decompose the velocity v as
, where is the solution of the unconstrained Stokes equation:
and w = w1|Σ = w2|Σ satisfies
Further, ΛGL is chosen such that
which implies that
To ensure local inextensibility, we use the correction ū which satisfies:
where (α, t) = ΛLL(α, t) − ΛGL(t). To find (α, t), we rewrite Eq. (30) as
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 . 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 (v · s)s + κv · n = 0 by introducing an appropriate tangential velocity ζ. That is, replace Eq. (32) by
where ζ is chosen so that the velocity field v is actually locally inextensible:
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 ; 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  and later in  for immersed interface method and in  in the context of the dynamics of inextensible, homogeneous vesicles. Ref.  is most relevant to our work here although unlike  we follow  and use the SSD to develop an explicit, nonstiff time integration algorithm.
where the second term in Eq. (43) is smooth. For α ≈ α′, we have,
where O(α – α′) denotes a smooth function that vanishes as α′ → α . Thus, is a smooth function.
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
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:
Further, from Eq. (27)
Next, we use the SSD to design an efficient, nonstiff explicit time integration scheme and an efficient preconditioner for the local Lagrange multiplier.
where we have absorbed the ΛLL contributions into N1 which is given by
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
We use the trapezoidal rule to approximate
The arclength α(t) can be calculated using the Adams–Bashforth method:
To reconstruct the membrane interface (x(α, tn+1), y(α, tn+1)) from the updated θn+1(α) and , we first update a reference point (x(0, tn+1), y(0, tn+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 by 
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 , 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  and use the SSD and the lemma to obtain a spectrally accurate discretization using the discrete Fourier transform. Finally, following , 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.
Replacing θα by a characteristic value α, one may define the operator
which can be inverted in Fourier space using the symbol:
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
we may rewrite Eq. (24) as:
Taking the Fourier transform of Eq. (66) to get
we can use an analogous linear propagator method to perform the time discretization:
In this section, we summarize the governing equations and numerical procedure. For the fluid, the equations are:
where and i = 1, 2 denotes the interior and exterior fluids. The function u is the concentration of one of the lipid phases, and ΛLL(α, t) is the Lagrange multiplier (tension) that enforces inextensibility ((u · s)s + κu · n = 0 on the membrane Σ). The surface phase evolves according to
where Pe is a dimensionless parameter, and the tangent angle θ to the interface evolves by
where V and T are the normal and tangential velocities given by
The arclength variation sα evolves according to
Even though Ts + κV) ≈ 0 by inextensibility, we still solve this equation to provide an assessment of the actual arclength variation during the evolution.
The outline of the nonstiff numerical algorithm used to solve the system is as follows:
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.1u and bending stiffness b(u) = (1 − u) + 0.5u. Also, we use a double well potential function f(u) = 0.25u2(1 − u)2. The initial condition is an ellipse given by:
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:
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 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.
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α(α, t) − sα(α, 0)| ≈ 10−8 in all cases. We study convergence of the area enclosed by the interface, which is more sensitive to the discretization, by comparing the area at time t = 0.02 with its initial value. The area is calculated by where the integration is performed using the discrete Fourier transform. The results are shown in Table 1 where the convergence rate is estimated from the errors corresponding to the results using two consecutive time step sizes. This confirms that the scheme is second-order accurate in time.
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 ). 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 EM is shown together with insets of the corresponding vesicle morphologies and surface concentration u at the times indicated (color online: blue and red correspond to u = 0 and u = 1 phases, respectively). At early times phase separation rapidly occurs yielding two large regions of the u ≈ 1 at the vesicle top and bottom where the curvature is smallest and the u ≈ 0 phase is confined to the vesicle tips where the curvature is larger. This is consistent with the desired spontaneous curvature (u). Correspondingly, there is a large drop in the total energy EM and as seen in (b) both energy components ET and Eb also decrease rapidly. At later times (t > 0.04), the upper and lower regions of the vesicle become negatively curved and the vesicle tends to a discocyte (or dumbbell) shape. The u = 1 phase remains confined to the negatively curved neck region. During this phase of the evolution, EM, ET and Eb all decrease slowly and the system tends to a steady state.
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.
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 EM, ET and Eb correspondingly decrease rapidly. Note that the concentration u is still quite far from 0 at the vesicle tips at t = 0.012, however. Because of asymmetries in the initial concentration, the concentration at the right tip is slightly larger than that at the left. This drives the u = 0 phase to the left tip as a way to further reduce the line energy ET. While this reduces ET, and ultimately EM, the bending energy Eb increases because it is more energetically favorable to have the u = 1 phase located in regions where the curvature is small rather than having the configuration shown where the u = 1 phase is located at the large curvature right tip. The line energy ET decreases dramatically because the four phase interface points on the membrane are reduced to two by this process. At later times, negative curvature regions form and the membrane tends an asymmetric discocyte shape. The shape, surprisingly, has the u = 1 phase located at the larger curvature tip. This is in contrast to the desired spontaneous curvature (e.g. K̄(1) = 0.1 and K̄(0) = 5) and is what drives the rapid increase in Eb. This configuration is a result of the limited amount of the u = 1 phase. For example, having a lower curvature tip requires more of the u = 1 phase. The normal velocity, shown in Fig. 4(c), indicates that the evolution still has not quite reached a steady state even though the energy variation is very small at late times. 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 8 × 10−6.
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 Eb at the expense of a slight increase in the line energy ET. As the phase separation process continues, the u = 0 phase accumulates at the left tip of the vesicle. This results in a large drop in ET and EM. The energy Eb slightly increases. The resulting vesicle takes an asymmetric discocyte shape but now with the larger curved region at the left tip. This is consistent with the desired spontaneous curvature and occurs because there is now enough u = 1 phase to support a less curved right tip. The normal velocity shown in Fig. 5(c) shows that the evolution has nearly reached the steady state. 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 5 × 10−7.
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.
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 EM, ET and Eb. At later times, the vesicle develops negative curvature. Because of the initial asymmetry in u, the concentration in the lower region is less than that in the upper region of the vesicle. This drives the transport of the u = 1 phase to the upper region and results in a rapid decrease in ET which overcomes the increase in Eb that is primarily associated with having the u = 0 phase located in the negatively curved lower region. The vesicle responds later by evolving so that the lower region becomes convex thereby decreasing the bending energy. The vesicle accordingly tends to an asymmetric (stomatocyte) equilibrium shape. In this simulation (Δt = 2 × 10−5, N = 512), the area is preserved to 6 digits of accuracy and the maximum error in local arclength is 1 × 10−6.
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/) ET and correspondingly the total energy EM. The system responds by rapidly phase separating, with the u = 1 phase now placed asymmetrically on small curvature regions of the vesicle. Once the phase separation occurs the u = 1 phases develop negative curvature and the vesicle evolves to a slightly asymmetric discocyte shape. Unlike the previous case where the u = 1 phase is transported to the top of the vesicle, here the smaller value of prevents this transport since decreasing enhances the immiscibility of the u = 0 and u = 1 phases. Note that if we had considered locally equilibrated initial conditions (well-defined and separated u = 0 and u = 1 domains) then we would have observed convergence of the results in . We have confirmed this (results not shown). Also observe from the normal velocity in Fig. 8(c) that the shape at t = 15 is still evolving apparently towards a more symmetric shape. In this simulation (Δt = 2 × 10−5, N = 512), the area is preserved to five digits of accuracy and the maximum error in local arclength is 1 × 10−6.
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
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 . 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.
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.
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.  for homogeneous vesicles although we use a different and more efficient time stepping algorithm  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 ) 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.
1Actually this system has a nullspace consisting of functions of the form f = cn where c is a constant. A unique solution perpendicular to this nullspace may be found.