PHYSICAL REVIEW E 91, 012708 (2015)

Budding of domains in mixed bilayer membranes
Jean Wolff,1,3,* Shigeyuki Komura,2,† and David Andelman3,‡
1

Institut Charles Sadron, UPR22-CNRS, 23 rue du Loess, B.P. 84047, 67034 Strasbourg Cedex 2, France
2
Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397, Japan3
Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel(Received 14 October 2014; published 20 January 2015)
We propose a model that accounts for the budding behavior of domains in lipid bilayers, where each of thebilayer leaflets has a coupling between its local curvature and the local lipid composition. The compositional
asymmetry between the two monolayers leads to an overall spontaneous curvature. The membrane free energy
contains three contributions: the bending energy, the line tension, and a Landau free energy for a lateral phaseseparation. Within a mean-field treatment, we obtain various phase diagrams which contain fully budded, dimpled,and flat states. In particular, for some range of membrane parameters, the phase diagrams exhibit a tricriticalbehavior as well as a three-phase coexistence region. The global phase diagrams can be divided into three types
and are analyzed in terms of the curvature-composition coupling parameter and domain size.
DOI: 10.1103/PhysRevE.91.012708

PACS number(s): 87.16.D−, 64.75.St, 87.14.Cc

I. INTRODUCTION

Biological membranes are multicomponent assemblies
typically composed of lipids, cholesterol, glycosugars, andproteins, whose presence is indispensable to the normal
functioning of living cells [1]. Given the complexity of
biological membranes, studies of model membranes have been
conducted in vitro in order to gain insight into the structuraland physical behavior of biomembranes. Many studies, in particular, over the last two decades, have focused on simplified
artificial systems containing vesicles in solution, composed ofternary mixtures of lipids and cholesterol [2,3]. By decreasing
temperature, the ternary mixtures undergo a phase separation
between a liquid-ordered (Lo ) phase and a liquid-disordered(Ld ) one [4,5]. Depending on thermodynamical parameters,
the liquid domains show one of three distinct domain shapes:
flat, dimpled (partially budded), or fully budded [6].
A theoretical model for domain-induced budding of planar
membranes was proposed by Lipowsky [7,8] and later extended for the case of closed vesicles [9,10]. In the model, thecompetition between the membrane bending energy and thedomain line tension leads to a budding transition under
the constraint of a fixed domain area. Hu et al. [11] proposeda mechanism based on this interplay, which stabilizes patternsof several domains on closed vesicles without requiring any
osmotically induced membrane tension.
Returning to the case of planar membranes, Lipowsky’s
model [12] predicts that (i) an initially flat domain deforms
spontaneously into a completely spherical bud when theinitial domain size exceeds a critical size, and (ii) dimpleddomains are stable only when the spontaneous curvature of
the bilayer membrane is nonzero. The latter prediction waslater re-examined [13], because dimpled domains are observed
experimentally in vesicles with no apparent spontaneous curvature [14]. In order to resolve this discrepancy, Rim et al. [15]considered the effect of adding an overall lateral tension acting
on the membrane and used ideas about entropy-driven lateral*

jean.wolff@ics-cnrs.unistra.fr
komura@tmu.ac.jp‡
andelman@post.tau.ac.il†

1539-3755/2015/91(1)/012708(10)tension that were originally proposed by Helfrich and Servuss[16]. The resulting phase diagram [15] contains regions ofstability for the three domain morphologies mentioned above.
Moreover, the effect of the lateral tension on budding was
discussed by Lipowsky and coworkers [7,17] as well as byothers [18,19]. In particular, the budding process requiresthat the activation energy has to exceed the energy barrier
associated with the surface tension.
In this paper, we propose a model that describes domaininduced budding in bilayers composed of a binary mixture oflipids. We suggest that dimpled domains can be formed andremain stable due to a possible asymmetry between the twomonolayer compositions. We show that the dimpled structureappears when the line tension along the domain rim is not too
great. Global phase diagrams are calculated within mean-field
theory, and in some range of system parameters, we obtain atricritical behavior as well as a three-phase coexistence region.
We discuss different morphologies that characterize the phase
diagrams in terms of model parameters.
It was recognized long ago that such an asymmetry in themonolayer composition leads to a nonzero bilayer spontaneous
curvature [20,21]. The coupling between composition andmonolayer curvature was also considered [22] in order todescribe the transition between lamellar and vesicular phasesof bilayer membranes composed of two types of amphiphiles.
It is worthwhile mentioning related works by Harden et al.
[23,24] and Góźdź et al. [25], who studied budding and domain
shape transformations in bilayer membranes. In Refs. [23] and[24], the phase separation is assumed to occur in only oneof the monolayers, and the domain spontaneous curvature
due to the compositional asymmetry is kept constant. For
finite spontaneous curvatures, it was shown that dimpleddomains are obtained in equilibrium when both line andsurface tensions are low [24]. More recently, ring-shaped
domains were experimentally obtained in model membranes
using a bud-mimicking topography [26]. Such a ring-shaped
domain is located around a bud-neck region having a negative
curvature and is characterized by the composition asymmetry
between the two monolayers.
The outline of our paper is as follows. In the next section,
we present a model for bilayer domains. In Sec. III, various

012708-1©2015 American Physical Society

JEAN WOLFF, SHIGEYUKI KOMURA, AND DAVID ANDELMANPHYSICAL REVIEW E 91, 012708 (2015)

mean-field phase diagrams are obtained by changing the ratio
between the domain size and the invagination length, as wellas tuning the intermonolayer coupling parameter strength.
Finally, Sec. IV includes some discussion and interpretationof our results.

φ1
φ2

II. MODEL

φ1 ≈ φ 2

We model the membrane as a bilayer having two monolayers (leaflets), each composed of an A/B mixture of lipids
that can partition themselves asymmetrically between the twomonolayers. We consider the case where the lipids can undergo
a lateral phase separation creating domains rich in one of thetwo components. As discussed below, these domains can also
deform (bud) in the normal direction, and the deformations arecontrolled by the membrane curvature elasticity. We do notinclude any gradient terms in the free energy; this results inan unrealistic discontinuous jump of the membrane curvature
close to the bud edge. For the fully budded state this jump incurvature does not matter because the bud neck corresponds toa small length scale, of the order of the membrane thickness.
However, in the dimpled state, this jump occurs on a bigger
length scale and artificially affects the free energy. We further
introduce a coupling between local lipid composition and local
curvature [20–22], which can eventually drive the budding
process of the membrane.
We start by considering a single two-dimensional (2D)
circular domain of an initial and arbitrary radius L embedded
in an otherwise flat (2D) membrane, as shown in Fig. 1(a). The
area of the domain, S = π L2 , is assumed to stay constant evenwhen the domain buds into the third dimension. For simplicity,
we consider only budded domains whose shape is a sphericalcap of radius 1/C [Fig. 1(b)]. The total bending energy of thebudded domain is given by adding the curvature contributions
from the two monolayers [27,28],
Eb = 2π L2 κ[(C − C0 )2 + (C + C0 )2 ],(1)

where κ is the bending rigidity modulus and C0 the monolayer
spontaneous curvature. As shown in Fig. 2, the two monolayers
are fully coupled together, and their curvatures are given by+C and −C, respectively.

FIG. 1. (Color online) (a) The flat phase, (b) the dimpled phase,and (c) the fully budded phase. In (a) the circular flat domain (red)has radius L and area S = π L2 . In (b) a bud of the same area S forms
a spherical cap of radius 1/C, where C is the curvature, embedded
in an otherwise flat membrane. In (c) a fully budded domain of area
S has a spherical shape, just touching the flat membrane. The line
tension γ acts along the boundary [solid (blue) line] between thedomain and the flat membrane.

φ1

φ2

φ 1 = φ2

FIG. 2. (a) Flat bilayer domain when the relative A/B compositions in the two monolayers, φ1 and φ2 , are almost symmetric,
φ1 ≈ φ2 . (b) Curved bilayer domain when the compositions areasymmetric, φ1 = φ2 . The spontaneous curvature of each monolayer
is assumed to depend linearly on the composition, as given by Eq. (4).

The next contribution is the domain edge energy, which is
proportional to the edge length and its line tension, γ [7]:(2)
Eed = 2π Lγ 1 − (LC/2)2 .
In the extreme case, when the domain buds into a complete
spherical domain as in Fig. 1(c), C = ±2/L and Eed = 0.
For domains that are composed of two lipid types, therelative composition in each monolayer is defined as φi =φiA − φiB (i = 1,2), where φiA (φiB ) is the molar fraction of the
A lipid (B lipid) in the ith monolayer. We assume that each of
the monolayers is incompressible, hence, φiA + φiB = 1. For
the sake of simplicity, the molecular areas of A and B species
are taken to be the same, meaning that the molar fraction of
the lipids is indenting to their area fraction. As in any A/B
mixture, the possibility of a phase separation due to partial
incompatibility between the two species can be described bya phenomenological Landau expansion of the free energy inpowers of φi around the critical point, φi = 0. In our case, thisexpansion is done separately for each monolayer, and the freeenergy is the sum of the two contributions,2 

t 2 1 42 U
(3)φ + φi − μφi ,Eph = π L 2
 i=1 2 i
4where  ≡ κ/γ is the invagination length, U a parameter thatsets the energy scale, t ∼ (T − Tc )/Tc the reduced temperature(Tc being the critical temperature), and μ the chemical
potential, which fixes the A/B relative composition in eachlayer. In general, a different chemical potential can be assigned
to each of the two monolayers. However, since it is difficult
to control the average composition in each layer separately,
we introduce only one chemical potential μ that fixes thetotal relative composition φ1 + φ2 of the entire bilayer. Notethat we allow exchange of lipid molecules between the twomonolayers via a flip-flop process. Bilayers where each of the012708-2

BUDDING OF DOMAINS IN MIXED BILAYER MEMBRANES

monolayer compositions can be controlled independently will
be addressed in our future work.
As argued before [22], we do not include any φi gradientterm in Eq. (3) because we consider only a homogeneouscomposition within a single domain. The energy cost associated with a gradient term in composition is effectively
taken into account through the line tension γ in Eq. (2),which is regarded here as an external control parameter. This
assumption regarding γ can be justified for situations of strong
segregation (far from the critical point) between the domain
and the background, for which the domain boundary is sharp.
Hereafter, we use several dimensionless variables: arescaled curvature c ≡ LC, rescaled spontaneous curvaturec0 ≡ LC0 , and rescaled invagination length ξ ≡ /L. The
coupling between the spontaneous curvature c0 and the composition is taken into account by assuming a linear dependenceon φi [22] (see also Fig. 2),
c0 (φi ) = c̄0 − βφi ,PHYSICAL REVIEW E 91, 012708 (2015)
III. PHASE BEHAVIOR AND PHASE DIAGRAMS
A. Flat, dimpled, and fully budded states

The total free energy ε in Eq. (7) is first minimized withrespect to the curvature c, yieldingc4c − 4βφ− − = 0.
(8)4ξ 1 − c2 /4
The above equation indicates that the value of c, taken at theminimum of ε, uniquely determines the value of φ− , as long(4)

where all variables in Eq. (4) are dimensionless, c̄0 is
the spontaneous curvature of the monolayer at its criticalcomposition φi = 0, and β is a coupling constant. Since c̄0 is aconstant that merely shifts the origin of the chemical potential
μ, we can drop it without loss of generality.
The total free energy of the bilayer model is given by thesum of Eqs. (1), (2), and (3):Etot = Eb + Eed + Eph .

(5)

Denoting the average of and difference between the twomonolayer compositions, respectively, as
φ+ ≡φ2 + φ1
2

and

φ− ≡

φ2 − φ1,
2(6)

the dimensionless total free energy of one domain, ε =Etot /2π κ, is expressed as
ε = 2c2 − 4βcφ− + 2β 2 (φ+2 + φ−2 )
 
1 U
1
2
t(φ+2 + φ−2 )1 − c /4 + 2
+
ξ
ξ
2κ

1
+ (φ+4 + 6φ+2 φ−2 + φ−4 ) − 2μφ+ ,2(7)where we have dropped unimportant constant terms. Within amean-field theory, the equilibrium state of the system and thephase transitions are determined by minimization of the above
ε with respect to φ± and c.
We note that Eq. (7) depends on three dimensionlessparameters — β, ξ , and U/2κ — while the thermodynamic
variables are the temperature t and the three order parameters,
φ± and c. In the calculations presented hereafter, we set
U/(2κ) = 1 and vary the values of β and ξ . Since the total freeenergy is invariant under simultaneous exchange of β → −βand φ− → −φ− , it is sufficient to study only the β > 0 range.
Typical experimental values of flat domain size are in
the range of L  50–500 nm [29], the bending rigidity κ 
10−19 J ≈ 25kB T [28], and the line tension is in the range ofγ  0.2–6.2 × 10−12 J/m [14,30]. These parameter values
yield an invagination length, , of the order 0.01L to 10L(ξ  0.01–10), as used in the next section.

FIG. 3. (Color online) Plots of the free energy ε as a functionof the curvature c for (a) φ− = 1, (b) φ− = 0.5, and (c) φ− = 0.
The other parameter values are t = −0.5, φ+ = 0.4, β = 1, and ξ =0.25. The free-energy minimum is shown by the filled (red) circle
and corresponds to the fully budded state [c = 2; (a)], dimpled state[0 < c < 2; (b)], and flat state [c = 0; (c)].

012708-3

JEAN WOLFF, SHIGEYUKI KOMURA, AND DAVID ANDELMANPHYSICAL REVIEW E 91, 012708 (2015)

FIG. 4. (Color online) Plot of the free energy ε as a function of
φ+ for t = −0.5, β = 1, and ξ = 0.25. The solid gray (red) andblack lines of the free energy correspond to the flat and fully budded
states, respectively. The two dotted lines are the common tangent
constructions, which determine two sets of coexisting compositions
indicated by green dots.

as |c| < 2. The value of the curvature determines which of thedomain states is the equilibrium one: flat (F) with c = 0, fully
budded (B) with c = ±2, or dimpled (D) with 0 < |c| < 2. By
substituting back the above minimization condition for |c| < 2
into the total free energy, we obtain ε(φ+ ,φ− ) as a functionof φ+ and φ− . This free energy is further minimized withrespect to φ− , leading to the expression ε(φ+ ), which is only afunction of φ+ . We assume that the average composition φ+ is aconserved order parameter (while c and φ− are nonconserved)
and can be controlled by varying the conjugate chemical
potential μ acting as a Lagrange multiplier.
In order to illustrate this minimization process, we plot ε asa function of the curvature c in Fig. 3, for given values of t andφ± . We see that the free energy takes its minimum at different
curvature values [shown by filled (red) circles] for different
φ− values. Figures 3(a), 3(b), and 3(c) correspond to the fully
budded, dimpled, and flat states, respectively. In Fig. 4 the freeenergy ε, which was minimized with respect to both c and φ− ,is plotted as a function of φ+ for a fixed temperature. Different
colors of the free-energy plot correspond to different domainstates (F and B). The two dotted lines are the common tangentsthat determine the two sets of coexisting compositions. For the
chosen parameter values as in Fig. 4, the flat and fully budded
phases coexist (F + B and B + F).
B. Phase diagrams

The numerically computed phase diagrams are 3D ones
for fixed values of ξ and β. They can be plotted in either the(φ+ ,t,c) or the (μ,t,c) parameter space. Recall that throughout
this work we set, for simplicity, U/(2κ) = 1. In addition,note that the equilibrium φ− value is self-determined by
the equilibrium c value according to Eq. (8). As it is too
cumbersome to present 3D plots, we plot 2D phase diagramsin the (φ+ ,t) or (μ,t) planes, which represent a projection in
the c direction, or 2D cuts in the (c,t) plane for fixed values of
the conserved order parameter, φ+ (see Fig. 8).

FIG. 5. (Color online) Phase diagram (a) in the (φ+ , t) plane,
where φ+ is the average composition and t the reduced temperature;and (b) in the (μ, t) plane, where μ is the chemical potential. The
parameters are β = 1 and ξ = 0.25. Flat (F), dimpled (D), and fully
budded (B) phases. Coexistence regions are denoted “F + D”, etc.,
in (a). The black and darker gray (red) lines indicate first- andsecond-order phase transitions, respectively, while the lighter gray(blue) line indicates a first-order phase transition with a discontinuous
jump in both c and φ− . Filled circles represent the tricritical points
tcp(ttcp  0.011, φ+  ±0.094, μtcp  ±0.488), and open circles in (b)represent the triple points (ttri  −0.267, μtri  ±1.49), where threephases coexist, with φ+tri = ±0.484, ±0.066, and ±0.045.

In Fig. 5 we present the phase diagrams that are obtained
numerically for β = 1 and ξ = 0.25. In (a) the phase diagramis plotted in the (φ+ ,t) plane; in (b), in the (μ,t) plane. The
phase diagram in (a) is symmetric about φ+ = 0, and that in(b), about μ = 0. At high temperatures, only the flat phaseis stable. For lower temperatures, in the range −0.267 < t <0.042, the dimpled phase becomes stable. The phase diagrams
show a tricritical behavior, similar to the well-known tricriticalbehavior of the Blume-Emery-Griffiths spin-1 model [31].
The darker gray (red) line in Fig. 5 denotes a second-orderphase transition between the F and the D phases, occurring when c → 0. It terminates at two symmetric tricritical

012708-4

BUDDING OF DOMAINS IN MIXED BILAYER MEMBRANESPHYSICAL REVIEW E 91, 012708 (2015)

FIG. 6. (Color online) Enlarged middle zone (φ+ ≈ 0 and taround the triple-point value) from Fig. 5(a). The solid gray (blue)line is a first-order phase transition between dimpled and fully budded
phases, with a jump in the curvature value from c  0.75 to c = 2.
Solid black lines indicate the boundaries of the two-phase coexistence
regions: F + D and F + B. On the triple line, ttri  −0.267, threephases coexist with φ+tri = ±0.066 and ±0.045.

tcp

points (filled circles), ttcp  0.011, φ+  ±0.094 in (a) andμtcp  ±0.488 in (b). The tricritical points are also obtained
analytically using some approximations and their calculated
tcp
values, ttcp  0.014 and φ+  ±0.095, agree well with thenumerical ones. More details on the analytical derivations areprovided in the Appendix. For t < ttcp , the phase transition
between F and D becomes first-order (solid black line) withcoexistence lines in the (μ,t) plane and two coexistence
regions, labeled F + D and D + F, in the (φ+ ,t) plane. As
one crosses this phase transition line, there is a jump in φ+ , aswell as in c and φ− , and the jump in φ− is fully determined bya similar jump in c.
Two triple points are shown as open circles in Fig. 5(b) at
ttri  −0.267 and μtri  ±1.49 or, equivalently, as a horizontal
line in Fig. 5(a). At the triple point, the three phases (F, D, andB) coexist. In order to explain in more detail the phase behavior
close to the triple line, we show in Fig. 6 an enlarged sectionof Fig. 5(a) around the triple line. The tip of the middle [gray(blue)] line starts at about t  −0.261 and terminates at the
triple-point temperature. This is a first-order phase transitionline where both c and φ− show a discontinuous jump fromtheir dimpled values (c  0.75) to their fully budded values(c = 2). The other solid lines delimit the two-phase coexistence
regions: F + D above the triple line and F + B below it.
We plot the equilibrium values of c in Fig. 7(a) and thoseof φ− in Fig. 7(b), in order to view the phase transitions moreclearly. Both c and φ− are plotted in the (φ+ ,t) plane as acontour color plot. In Fig. 7(a) we see two parabola-like lines
delimiting different values of c. At the upper black line, thecurvature continuously tends towards 0, c → 0. The region
close to the curve tip (φ+ ≈ 0) coincides with the second-orderphase transition between phase F and phase D [darker gray(red) line in Fig. 5(a)], while the rest of the line lies inside
the two-phase coexistence region, and does not influence theequilibrium state of the system. The lower, lighter gray (blue)line represents a jump in c from c  0.75 (D phase) to c = 2(B phase). Its top region (close to φ+ = 0) coincides withthe first-order phase transition between phase D and phase B

FIG. 7. (Color online) (a) The curvature c as a contour plot inthe (φ+ , t) plane, with a color bar that corresponds to 0  c  2.
A jump from c = 2 to about c = 0.75 can be seen as t increases
its value and eventually crosses the lower black line, while along
the upper black line c vanishes continuously. (b) The compositional
asymmetry between the two monolayers, φ− , is plotted as a contour
plot in the (φ+ , t) plane. As t increases, a jump is seen from φ− = 0.7
to φ− = 0.55 along the lower black curve. For the upper black curveφ− vanishes continuously (just as c did). The chosen parameters are
β = 1 and ξ = 0.25.

[lighter gray (blue) lines in Figs. 5(a) and 6], and the rest of theline lies within the F + B coexistence region. In Fig. 7(b), a
similar contour plot is shown for φ− , as determined by Eq. (8).
The complementary plot is shown in Fig. 8 in the (c,t) plane
for several fixed values of φ+ , ranging from 0.1 to 0.4. We
recall that φ+ = (φ1 + φ2 )/2 = (φ1A + φ2A )/2 − (φ1B + φ2B )/2
is a conserved quantity determined by the total amount of the
A and B lipids in the domain. In the model we control it by thechemical potential μ. As t is lowered, the minimized curvature
c continuously increases from 0. This represents a continuous(second-order) phase transition from the flat state (F with c =0) to the dimpled one (D with c > 0). When the temperature
is lowered even further, the curvature discontinuously jumps
from c  0.75 (indicated by the arrow in Fig. 8) to c = 2. Thisis a first-order phase transition from the dimpled state (D) to

012708-5

JEAN WOLFF, SHIGEYUKI KOMURA, AND DAVID ANDELMANPHYSICAL REVIEW E 91, 012708 (2015)

FIG. 8. The reduced temperature t as function of the equilibrium
curvature c for fixed values of β = 1, ξ = 0.25 and several values of
φ+ : 0.1, 0.2, 0.27, 0.33, and 0.4. As t decreases, the curvature c firstcontinuously increases from 0 and then discontinuously jumps fromc  0.75 (marked by the arrow on the c axis) to c = 2. The former
corresponds to the second-order phase transition from the flat state to
the dimpled one; the latter, to the first-order phase transition from thedimpled state to the fully budded one.

the fully budded (B) one. Note that the maximum curvaturec  0.75 of the dimpled state does not depend on the averagecomposition φ+ .
When the ξ value is decreased, while β is kept fixed, theD phase disappears, and the only remaining stable phases are
F and B, with a phase transition between them. This is shownin Fig. 9, where the chosen parameter values are β = 1 andξ = 1/7  0.143. A second-order phase transition [gray (red)line] is seen between the F and the B phases in the proximity ofthe symmetric φ+ = 0 axis. This second-order line ends at two
tcp
tricritical points, located at ttcp  0.240 and φ+  ±0.063 in(a) or, equivalently, μtcp  ±2.64 in (b). Below the tricriticaltemperature, the coexistence region is between the F and theB phases (F + B) and is delimited by the solid black lines.
Note that as the D phase disappears there is no three-phasecoexistence at these parameter values. The disappearance of
the D phase can be understood in the following way. Smaller
values of the invagination length, ξ = κ/(Lγ ), correspond tolarger values of the line tension γ , and domains will fully bud
at lower temperatures without showing any D state.
At yet lower values of ξ , the line tension is large enoughso that only the B phase exists, while the F phase disappears.
In Fig. 10, we present such a phase diagram for β = 1 andξ = 0.125. The only coexistence regions are between different
fully budded phases, denoted B1 + B2 and B2 + B3. Each ofthese coexistence regions terminates at critical points (filledsquares), tc  −0.028, φ+c  ±0.251, and μc  ±8.
C. Effects of ξ and β on the phase behavior

By exploring the entire parameter range of ξ and β, wefind the crossover among the three types of phase diagramsas represented in Fig. 5 (type I), Fig. 9 (type II), and Fig. 10(type III). This is shown in Fig. 11, where we present the

FIG. 9. (Color online) Phase diagram (a) in the (φ+ , t) planeand (b) in the (μ, t) plane for β = 1 and ξ = 1/7  0.143. The
meaning of the lines and symbols is the same as in Fig. 5. A
critical [gray (red)] line separates the F and B phases and terminatestcp
at two tricritical points (filled circles) with (ttcp  0.240, φ+ 
±0.063, μtcp  ±2.64). For t < ttcp , coexistence regions, B + F and
F + B, separate the F and B phases.

stability regions for each of these three phase behaviors inthe (ξ , β) plane. Type I is characterized by the existence ofa dimpled phase and has a triple point where the D, B, andF phases coexist. In type II, the tricritical points exist but the
triple points and the D phase disappear. Type III is dominatedby various B phases, with coexistence regions between themthat terminate at a critical point. The crossover line between
type II and type III behaviors is almost a straight line, while
the crossover line from type II to type I is almost linear forβ < 1.25 and then saturates at about β  1.8. This saturation
occurs when the coupling is strong (large β) and/or the domain
size is small (large ξ ). At these values, budding is promoted
because of the large spontaneous curvature.
When ξ decreases, for a fixed value of the coupling
parameter β, the B phase swells and the D phase disappears,
signaling the crossover between type I and II. Upon further
decrease of ξ , only the B phase stays, i.e., crossover between

012708-6

BUDDING OF DOMAINS IN MIXED BILAYER MEMBRANESPHYSICAL REVIEW E 91, 012708 (2015)

FIG. 11. (Color online) Behavior diagram as a function of ξ =/L and the monolayer coupling parameter β. Examples of types I,
II, and III phase-diagram morphologies are shown in Figs. 5, 9, and10, respectively. Filled (red) squares delimit the calculated borderline
between type I and type II behavior, while filled (green) circles delimit
the calculated crossover between type II and type III. Dotted lines
serve only as a guide for the eye, and the error bar of the data points
is about ±0.1 in β.

FIG. 10. Phase diagram (a) in the (φ+ , t) plane and (b) in the(μ, t) plane for β = 1 and ξ = 0.125. The meaning of the lines andsymbols is the same as in Fig. 5. Only the B domain is stable, afirst-order phase transition separates B1 and B2, and another one
separates B2 and B3. Each of the coexistence regions terminates at a
critical point (tc  −0.028, φ+c  ±0.251, μc  ±8). Filled squares
correspond to critical points.

type II and III. On the other hand, the larger β is, the larger isthe spontaneous curvature that favors the fully-budded state.
For this reason, at higher values of β, the system buds at lower
temperatures for the same value of ξ .
IV. DISCUSSION

We have proposed a model that accounts for domain
budding of lipid bilayers, where each of the bilayer leaflets
has a coupling between its local curvature and the local A/B
lipid composition. The composition asymmetry between the
two leaflets is equivalent to the introduction of a membranespontaneous curvature. This spontaneous curvature is not taken
to be fixed (as assumed in previous works) but is calculated anddepends on the asymmetry in leaflet composition. Hence, due
to this extra mechanism of generating a spontaneous curvature,

dimpled domains can be stabilized even for bilayers with anominal zero spontaneous curvature. Our free-energy model
contains three contributions: bending energy, accounting fordomain deformation in the normal direction; line tension along
the rim of the budded or flat domain; and a Landau free-energy
expansion, which accounts for a lateral phase separation of thebinary lipid mixture. We assume that the domain area remains
constant during the budding process.
Our model predicts three states for domain, as observed
experimentally: fully budded (B), dimpled (D), and flat (F)states. In particular, in some ranges of parameters, the Dstate is found to be the most stable one. The obtained results
indicate that for a certain range of temperatures, monolayercomposition, domain size, and coupling between curvatureand composition, a triple point can appear. At the triple point,
the B, D, and F phases coexist, each with its own composition.
Such a triple point has been reported by Harden et al. [24].
Moreover, we also found a tricritical point, which corresponds
to the intersection of a critical (second-order) line, which joinsa first-order phase transition region between F and D. Finally,
three types of phase diagram morphologies are found andanalyzed in terms of the coupling parameter β and domain
size ξ (see Fig. 11).
The formation of domains in membranes and understanding
them remain an open and ever-challenging problem, even after
intense research in the last two decades. Many hypotheses havebeen proposed to explain the appearance of domains and their
possible structure and function [2,3,32].
One of the important assumptions in our model is thatthe domain size L and area S remain fixed during thebudding process. Domains of fixed size can be obtained inthermodynamical equilibrium for a binary mixed membranethat undergoes a lateral microphase separation and forms a
2D modulated phase with an equilibrated spatial periodicity

012708-7

JEAN WOLFF, SHIGEYUKI KOMURA, AND DAVID ANDELMAN[33]. This microphase separation can be driven by a coupling
between the local lipid composition and the membrane curvature, leading to a curvature instability [34–37], as discussed, inparticular, in Refs. [38] and [39]. When the A/B average lipid
composition is off-critical, circular domains rich in one of thelipids can form spontaneously and be arranged in a hexagonalarray, embedded in a background rich in the second lipid.
These circular domains are characterized by their equilibrium
fixed size and can undergo a budding process as explored in
the present work.
Formation of finite-size domains in equilibrium can also be
explained by the presence of hybrid lipids having one saturated
tail and a second unsaturated one [40]. Such hybrid lipids
decrease the domain line tension [41–43] and offer anotherpotential mechanism to induce micro phase separation. In aprevious work [44,45], we considered a model that includes
a coupling between a compositional scalar field and a 2D
vectorial order parameter. This coupling yields an effective 2D
free energy that exhibits micro phase separation and results
in a modulated phase. A somewhat different viewpoint of
membrane domains was recently discussed by Shlomovitz
et al. [46], who investigated a general phenomenologicalmodel capable of producing macro phase separation, microphase separation, and microemulsion-like phases. In theseworks, the characteristic length of compositional modulations
is responsible for the origin of finite-size domains that are
equilibrium structures. These types of domains can undergo
the budding transition as discussed in this paper.
On the other hand, when macro phase separation takes
place in mixed membranes, the domain size grows to
macroscopic sizes as a function of time, and the assumption of
fixed domain size becomes more questionable. However, if the
shape transformation of domains occurs on time scales much
faster than the time required for domain coarsening, one canstill use our equilibrium argument for domain morphologies,
whereby we regard the domain size L and, hence, theinvagination length ξ to be time dependent. In fact, theslowdown of dimpled-domain coarsening was experimentallyobserved [5] and theoretically discussed [13]. According tothese works, the suppression of the phase separation may
be caused by membrane-mediated elastic interactions and/orhydrodynamic interactions acting between domains. Ourmodel, with its assumption of fixed domain size, can also be
applied in such situations.
A dynamical growth of the budding domain size wasproposed [7] to occur in two steps, when the spontaneous
curvature is not too large. In the early stage the domains
are small and the diffusion-aggregation phenomenon induces
a growing dimpled domain until its size becomes unstable,
whereas in the later stage, the domains are large enough tofully bud into a sphere that detaches completely from theplanar membrane at the neck point. In other words, the budded
domain curves mainly during the second step. Our results
are in qualitative agreement with these predictions. We show
in Fig. 11 that, for small values of β and for small domain
sizes (large ξ ∼ 1/L), the typical phase diagram is of typeI, for which the dimpled state appears as a stable phase. As
the domain size becomes larger (smaller ξ ), the typical phase
diagram will change to either type II or type III so that themembrane can bud easily. In contrast, for large values of β,

PHYSICAL REVIEW E 91, 012708 (2015)

the dimpled state cannot be stable. In this case, the domain willretain a highly curved state already in the early stage of thephase separation, and the budding can occur in only one step.
Finally, our model may be applied to the description of
the formation and growth of vesicles in mixed amphiphilic
systems [47]. For example, it was observed in experimentsthat mixtures of anionic and cationic surfactants in solution
form disk-like bilayers for some range of relative surfactant
compositions. As these disk-shaped bilayers grow in size,
they transform into spherical caps and eventually become
spherically closed vesicles. Such a sequence of morphologicalchanges was indeed observed by cryo–transmission electron
microscopy [47]. In such a setup, it is likely that the
spontaneous curvature of bilayer membranes is induced dueto the compositional asymmetry between the two monolayers.
Hence, one would expect that disks, caps, and vesicles canbe analyzed similarly to the flat, dimpled, and fully budded
phases in our model.
ACKNOWLEDGMENTS

We thank T. Kato, B. Palmieri, and S. A. Safran for
useful discussions and numerous suggestions. J.W. acknowledges support from the Service de Coopération Scientifiqueet Universitaire de l’Ambassade de France en Israël, theFrench ORT association, and ORT School of Strasbourg. S.K.
acknowledges support from the Grant-in-Aid for Scientific Research on Innovative Areas “Fluctuation & Structure” (Grant
No. 25103010) from the Ministry of Education, Culture,
Sports, Science, and Technology of Japan, the Grant-in-Aid
for Scientific Research (C) (Grant No. 24540439) from theJapan Society for the Promotion of Science (JSPS), and the
JSPS Core-to-Core Program “International Research Network
for Non-equilibrium Dynamics of Soft Matter”. D.A. acknowledges support from the Israel Science Foundation under Grant
No. 438/12 and the United States–Israel Binational Science
Foundation under Grant No. 2012/060.
APPENDIX: THE TRICRITICAL POINT

It is possible to compute analytically the location of thetricritical point (ttcp ), corresponding to the intersection of thefirst- and second-order transition lines in the phase diagram in
Fig. 5(a). The left side of the binodal line corresponds to c = 0(F phase), while its right side corresponds to c > 0 (D phase).
Using the fact that the F phase with c = 0 has two symmetricmonolayers, and that φ− = 0 from Eq. (8), we calculate thefree energy for c = φ− = 0 by substituting c = 0 in Eq. (7):ε(φ+ ,c = 0) = 12 ξ −2 φ+4 + (2β 2 + tξ −2 )φ+2 + ξ −1 .

(A1)

The free-energy expression can then be expanded up to fourthorder in c (valid close to the tricritical point where c  1),
yielding
ε(c,φ+ ,φ− ) = 12 ξ −2 (φ+4 + φ−4 + 6φ+2 φ−2 )+ (2β 2 + tξ −2 )(φ+2 + φ−2 ) − 4cβφ−

1 4
(A2)
+ 2c2 + ξ −1 1 − 18 c2 − 128c .
From Eq. (8) we can expand c up to third order in φ− : c 
aφ− + bφ−3 . Substituting this c expression back into Eq. (A2)012708-8

BUDDING OF DOMAINS IN MIXED BILAYER MEMBRANESPHYSICAL REVIEW E 91, 012708 (2015)

and retaining terms up to fourth order in φ− , we can expand
ε(φ+ ,φ− ), obtaining a fourth-order polynomial in both φ+ and
φ− ,


ab −1
a 4 −1 1 −2ε(φ+ ,φ− ) = φ−4 4ab − 4bβ −
ξ −
ξ + ξ
4
128
2

a2
+ φ−2 2a 2 + 2β 2 − 4aβ − ξ −18
1+ tξ −2 + 3φ+2 ξ −2 + ξ −2 φ+42+ 2β 2 + tξ −2 φ+2 + ξ −1 ,(A3)where the coefficients a and b are defined as a = 16βξ/(16ξ − 1) and b = a 3 /(128ξ − 8).
The free energy, Eq. (A3), is then minimized with respectto φ− , yielding(A4)
φ−2 = δ − tξ −2 − 3φ+2 ξ −2 /η,with new coefficients δ and η defined as δ = −2a 2 −2β 2 + 4aβ + a 2 ξ −1 /8 and η = 8ab − 8bβ − abξ −1 /2 −a 4 ξ −1 /64 + ξ −2 . Substituting the expression of φ− in Eq. (A4)
into Eq. (A3), we get1 2δ − tξ −2 − 3ξ −2 φ+2 + 2β 2 + tξ −2 φ+2ε (φ+ )  −2η
1+ ξ −2 φ+4 + ξ −1 ,(A5)2which is valid in the limit c  1.
At the tricritical point, the free energies for c = 0 and c  1
are equal for the same values of φ+ and t. Comparing Eqs. (A5)and (A1), we obtain the condition for the tricritical point:δ − tξ −2 − 3ξ −2 φ+2 = 0.

(A6)

In addition, the spinodal line is obtained from the requirementthat ε (φ+ ) = 0 in Eq. (A5):93δ − tξ −2 + 2β 2 ξ 2 + t = 0. (A7)
3φ+2 1 − ξ −2 +ηη

[1] B. Alberts, A. Johnson, P. Walter, J. Lewis, and M. Raff,
Molecular Biology of the Cell (Garland Science, New York,
2008).
[2] S. L. Veatch and S. L. Keller, Biochim. Biophys. Acta 1746, 172(2005).
[3] S. Komura and D. Andelman, Adv. Coll. Interface Sci. 208, 34(2014).
[4] S. L. Veatch and S. L. Keller, Phys. Rev. Lett. 94, 148101 (2005).
[5] M. Yanagisawa, M. Imai, T. Masui, S. Komura, and T. Ohta,
Biophys. J. 92, 115 (2007).
[6] T. S. Ursell, Bilayer Elasticity in Protein and Lipid Organization(VDM Verlag, Berlin, 2009).
[7] R. Lipowsky, J. Phys. II (France) 2, 1825 (1992).
[8] R. Lipowsky, Biophys. J. 64, 1133 (1993).
[9] F. Jülicher and R. Lipowsky, Phys. Rev. Lett. 70, 2964 (1993).
[10] F. Jülicher and R. Lipowsky, Phys. Rev. E 53, 2670 (1996).

FIG. 12. Tricritical temperature ttcp as function of ξ , calculated
from Eq. (A8), for a fixed value of the coupling parameter, β = 1.
The analytical expansion is valid only for ξ values that are not too
small. In particular, below ξ = 0.2 the overall change in the phase
diagram is from type I to type II, as shown in Fig. 11.

By combining Eq. (A6) with Eq. (A7), we obtainttcp = ξ 4 6δξ −2 − δη − 2β 2 η /6,(A8)

which depends on the values of ξ and β. In Fig. 12 we plot thevariation of ttcp with respect to ξ for β = 1 and U/(2κ) = 1.
tcp
Substituting β = 1 and ξ = 0.25 in Eq. (A8), we obtain φ+ 
±0.095 and ttcp  0.014. These tricritical values are in good
agreement with the numerical ones, as can be read off from thetcpphase diagram in Fig. 5(a), φ+  ±0.094 and ttcp  0.011.
The present argument is valid only when c  1 and canbe applied to a type I phase diagram, for which the domain
size should be small (large ξ ). For a type II phase diagram,on the other hand, the tricritical points are located close to thefully budded phase with c = 2 (see Fig. 9). In this case, theexpansion in terms of c cannot be justified, and ttcp (ξ ) is nolonger valid for ξ < 0.2 when β = 1.

[11] J. Hu, T. Weikl, and R. Lipowsky, Soft Matter 7, 6092 (2011).
[12] R. Lipowsky, Biol. Chem. 395, 253 (2014).
[13] T. S. Ursell, W. S. Klug, and R. Phillips, Proc. Natl. Acad. Sci.
USA 106, 13301 (2009).
[14] T. Baumgart, S. T. Hess, and W. W. Webb, Nature 425, 821(2003).
[15] J. E. Rim, T. S. Ursell, R. Phillips, and W. S. Klug, Phys. Rev.
Lett. 106, 057801 (2011).
[16] W. Helfrich and R. M. Servuss, Nuovo Cimento 3, 137 (1984).
[17] R. Lipowsky, M. Brickmann, R. Dimova, C. Haluska,
J. Kierfeld, and J. Schillcock, J. Phys.: Condens. Matter 17,
S2885 (2005).
[18] P. Sens and M. S. Turner, Biophys. J. 86, 2049 (2004).
[19] P. Sens and M. S. Turner, Phys. Rev. E 73, 031918 (2006).
[20] S. A. Safran, P. Pincus, and D. Andelman, Science 248, 354(1990).

012708-9

JEAN WOLFF, SHIGEYUKI KOMURA, AND DAVID ANDELMAN[21] S. A. Safran, P. A. Pincus, D. Andelman, and F. C. MacKintosh,Phys. Rev. A 43, 1071 (1991).
[22] F. C. MacKintosh and S. A. Safran, Phys. Rev. E 47, 1180
(1993).
[23] J. L. Harden and F. C. MacKintosh, Europhys. Lett. 28, 495(1994).
[24] J. L. Harden, F. C. MacKintosh, and P. D. Olmsted, Phys. Rev.
E 72, 011903 (2005).
[25] W. T. Góźdź and G. Gompper, Europhys. Lett. 55, 587(2001).
[26] Y.-S. Ryu, I.-H. Lee, J.-H. Suh, S. C. Park, S. Oh, L. R. Jordan,
N. J. Wittenberg, S.-H. Oh, N. L. Jeon, B. Lee, A. N. Parikh,and S.-D. Lee, Nature Comm. 5, 4507 (2014).
[27] W. Helfrich, Z. Naturforsch. 28c, 693 (1973).
[28] S. A. Safran, Statistical Thermodynamics of Surfaces,
Interfaces, and Membranes (Addision Wesley, Reading, MA,
1994).
[29] K. Simons and E. Ikonen, Nature 387, 569 (1997).
[30] A. Tian, C. Johnson, W. Wang, and T. Baumgart, Phys. Rev.
Lett. 98, 208102 (2007).
[31] M. Blume, V. J. Emery, and R. B. Griffiths, Phys. Rev. A 4, 1071
(1971).
[32] D. Lingwood and K. Simons, Science 327, 46 (2010).
[33] M. Seul and D. Andelman, Science 267, 476 (1995).

PHYSICAL REVIEW E 91, 012708 (2015)[34] S. Leibler and D. Andelman, J. Phys. (Paris) 48, 2013 (1987).
[35] D. Andelman, T. Kawakatsu, and K. Kawasaki, Europhys. Lett.
19, 57 (1992).
[36] T. Taniguchi, K. Kawasaki, D. Andelman, and T. Kawakatsu,
J. Phys. II (France) 4, 1333 (1994).
[37] T. Kawakatsu, D. Andelman, K. Kawasaki, and T. Taniguchi,
J. Phys. II (France) 3, 971 (1993).
[38] H. Kodama and S. Komura, J. Phys. II (France) 3, 1305
(1993).
[39] P. B. Sunil Kumar, G. Gompper, and R. Lipowsky, Phys. Rev.
E 60, 4610 (1999).
[40] R. Brewster, P. A. Pincus, and S. A. Safran, Biophys. J. 97, 1087
(2009).
[41] T. Yamamoto, R. Brewster, and S. A. Safran, Europhys. Lett.
91, 28002 (2010).
[42] B. Palmieri and S. A. Safran, Langmuir 29, 5246 (2013).
[43] B. Palmieri and S. A. Safran, Langmuir 30, 11734 (2014).
[44] Y. Hirose, S. Komura, and D. Andelman, ChemPhysChem 10,2839 (2009).
[45] Y. Hirose, S. Komura, and D. Andelman, Phys. Rev. E 86,
021916 (2012).
[46] R. Shlomovitz, L. Maibaum, and M. Schick, Biophys. J. 106,
1979 (2014).
[47] A. Shioi and T. A. Hatton, Langmuir 18, 7341 (2002).

012708-10