PHYSICAL REVIEW E 94, 062414 (2016)

Relaxation dynamics of a compressible bilayer vesicle containing highly viscous fluid
T. V. Sachin Krishnan,1,2 Ryuichi Okamoto,1 and Shigeyuki Komura1,*
1

Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397, Japan2
Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India(Received 12 October 2016; published 30 December 2016)
We study the relaxation dynamics of a compressible bilayer vesicle with an asymmetry in the viscosity of theinner and outer fluid medium. First we explore the stability of the vesicle free energy which includes a coupling
between the membrane curvature and the local density difference between the two monolayers. Two types of
instabilities are identified: a small wavelength instability and a larger wavelength instability. Considering the bulk
fluid viscosity and the inter-monolayer friction as the dissipation sources, we next employ Onsager’s variational
principle to derive the coupled equations both for the membrane and the bulk fluid. The three relaxation modes
are coupled to each other due to the bilayer and the spherical structure of the vesicle. Most importantly, a higher
fluid viscosity inside the vesicle shifts the crossover mode between the bending and the slipping to a larger value.
As the vesicle parameters approach the unstable regions, the relaxation dynamics is dramatically slowed down,
and the corresponding mode structure changes significantly. In some limiting cases, our general result reduces tothe previously obtained relaxation rates.
DOI: 10.1103/PhysRevE.94.062414
I. INTRODUCTION

Macromolecules such as proteins, nucleic acids, andpolysaccharides are present in high concentration in cells andphysically occupy up to 40% of the total cell volume [1,2].
This feature, referred to as “macromolecular crowding,” has asignificant effect on important processes taking place inside thecell. This is due to the fact that any reaction which depends on
the available volume and the diffusivity of reactants is affectedby crowding [3]. The increase in the effective viscosity of thecrowded environment is considered to be responsible for theobserved anomalous diffusion of proteins inside cells [4,5]. It
is known that crowding influences the conformation, stability,
and the function of nucleic acids [6] and proteins [7]. However,most of the in vitro biochemical and biophysical studies
involving macromolecules are carried out in highly diluteenvironments which do not reflect the natural conditions inside
the cells [8].
Artificial cells with crowded internal environment are being
developed in order to study the macromolecules in their
physiological conditions [9,10]. Giant unilamellar vesicles(GUVs), which are cell-sized spherical shells consisting ofa lipid bilayer, are commonly used as containers for artificial
cells in which many synthetic chemical reactions take place.
In a recent paper by Fujiwara and Yanagisawa [11], they used
GUVs filled with macromolecules in order to investigate theeffect of higher viscosity inside the vesicles. When these GUVs
were subjected to osmotic pressure, they exhibited different
shape deformations depending on the inside viscosity. This
result suggests that the asymmetry in the viscosity affects theshape dynamics of the bounding lipid membrane [11].
One of the early quantitative works on dynamics ofmembranes is by Brochard and Lennon, who studied theflickering phenomena in red blood cells (RBCs) [12]. They
identified the observed flickering of RBC as equilibrium
thermal fluctuations of flexible lipid bilayers [13]. Analysis*

komura@tmu.ac.jp

2470-0045/2016/94(6)/062414(14)of flickering phenomena is a useful method to estimate elastic
properties of fluid membranes. Furthermore, measurementsof membrane flickering can be used to identify diseases
affecting biological cells. On the other hand, recent experiments revealed that the flickering is also controlled by active
ATP-dependent processes taking place inside the cell [14–16].
Even though the relaxation dynamics of fluid vesicles has
been studied theoretically for several decades [17–21], theeffect of fluid viscosity difference between the inside andoutside of the vesicle has not received much attention. Seki andKomura explicitly took into account the viscosity difference,and obtained an expression for the bending relaxation modeof quasispherical vesicles [22–24]. However, since the bilayer
nature of the vesicle was not considered in their calculations,the relaxation of the lipid densities on the monolayers mediatedby the inter-monolayer friction was not taken into account.
The relaxation mode governed by this inter-monolayer friction(called the “slipping mode”) affects the small wavelength
dynamics of planar membranes [25–28].
In this work, we investigate the relaxation dynamics ofa vesicle by explicitly taking into account the effect of
fluid viscosity contrast between the inside and outside of
the vesicle. We consider a quasispherical, single-component,
unilamellar, and compressible bilayer vesicle. The vesicle
shape deformation and the difference in the local densities ofthe monolayers are governed by the bending rigidity and themembrane compressibility, respectively. The dissipations arise
from the viscosity of the bulk fluids and the inter-monolayer
friction. We employ Onsager’s variational principle to derive
the relaxation equations for the system that can be described by
low Reynolds number hydrodynamics [29,30]. We show thatthe bending and the slipping modes are coupled to each otherdue to the bilayer architecture. The viscosity contrast results inslowing down of the bending mode, and it significantly affects
the slowest relaxation mode. We also discuss the stability ofthe vesicle free energy, and identify two distinct instabilities.
One of them is associated with the slowing down of the largest
wavelength relaxation mode, while the other is related to theslowing down of the shorter modes.

062414-1©2016 American Physical Society

SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)

In previous works, exact expressions and approximations
for the coupled relaxation modes of a compressible vesicle
have been discussed [19,20]. Our goal is different fromthese works in that we account for the viscosity asymmetry
between the inside and outside fluid media. In addition,
we discuss the conditions for stability of the vesicle freeenergy. Our calculation also follows a different method basedon Onsager’s variational principle in order to derive theappropriate dynamical equations. We also comment that theeffect of different viscosities on the dynamics of vesicles
under shear flow was studied previously [31–33]. It was foundthat such a viscosity contrast causes a transition between
tank-treading and tumbling motions [32,33]. While vesicles
under flow in these works attain one of the nonequilibrium
steady states in the long time limit, vesicles in our study
eventually relax to a mechanical equilibrium described by aperfect sphere.
The outline of this paper is as follows. Section II introduces the parametrization to specify the configuration of
the vesicle and its free energy functional. In Sec. III, we
examine the statics of the system near equilibrium, and derive
the conditions for thermodynamic stability. In Sec. IV, weconstruct the Rayleighian functional and derive the equations
for the relaxation dynamics by extremalizing the Rayleighian.
Section V presents the main results of the work, including theeffect of viscosity contrast on the relaxation dynamics of thebilayer vesicle. Finally, some discussion on the relation of ourwork to the previous theoretical calculations are provided inSec. VI.
II. FREE ENERGY OF A COMPRESSIBLE BILAYERVESICLEA. Variables

A vesicle is composed of two opposing layers of lipid
molecules whose tails meet at the bilayer midsurface. For anearly spherical vesicle, this surface can be mathematicallydescribed by an angle-dependent radius field r(θ,ϕ). We
assume that the deformation of the vesicle is small enoughso that the function r(θ,ϕ) is a continuous single-valued one.
As shown in Fig. 1(a), we define a dimensionless deviation inradius u by using the radius of a reference sphere r0 as [17,34]u(θ,ϕ) =r(θ,ϕ)
− 1.
r0

(1)

A consequence of the membrane bilayer structure is thatbending deformation always accompanies stretching of one
monolayer and compression of the other, as shown in Fig. 1(b).
Hence the membrane monolayers are weakly compressible as
was found in the experiments [35,36]. To take into account
this compressibility, the local lipid density on each monolayer
is allowed to vary. Let the number of lipids in the outer andinner monolayers be N + and N − , respectively. Throughout thispaper, we assign the superscript “+” to variables associatedwith the outer monolayer and the fluid outside the vesicle.
Similarly, we assign the superscript “−” to those associatedwith the inner monolayer and the fluid enclosed by the vesicle.
Since the spherical vesicle is stable only if the outer monolayer
has more lipids than the inner one, we impose the restriction
N + > N − . We assume that the flip-flop motion of the lipids

FIG. 1. (a) Schematic picture of a vesicle showing the reference
sphere (dashed line) of radius r0 with fluids of viscosity η− insideand η+ outside. The solid red line represents a deformed vesicle
configuration whose shape is parametrized with u(θ,ϕ). (b) Crosssection of a curved region in the bilayer showing the lipid moleculesin each monolayer. Local densities in the monolayers are indicated as
φ + and φ − which are defined on the bilayer midsurface represented
by the dashed line. The inter-monolayer friction w also acts at thebilayer midsurface.

between the two monolayers does not occur within the timescale of experimental observation. This means that both N +and N − are conserved quantities. We denote a reference lipid
density by ρ0 = (N + + N − )/(2A0 ), where A0 = 4π r02 is thetotal area of the vesicle in the reference state. Let ρ + (θ,ϕ)
and ρ − (θ,ϕ) be the variables representing local lipid density
in each of the monolayers. We then define dimensionless local
density deviation as
φ ± (θ,ϕ) =ρ ± (θ,ϕ)
− 1,
ρ0

(2)

in each monolayer. Notice that we have defined the densities
on the bilayer midsurface rather than the neutral surface ofa monolayer. On the neutral surface, stretching and bending
deformations are decoupled by its definition.
Further calculations can be simplified by using the density
difference and the density sum, defined as

062414-2

φ  (θ,ϕ) =φ+ − φ−,2

φ (θ,ϕ) =φ+ + φ−,2(3)

RELAXATION DYNAMICS OF A COMPRESSIBLE BILAYER . . .

respectively. At equilibrium, the vesicles take a sphericalshape with uniform lipid density in each monolayer. In thishomogeneous state, the two variables defined above becomeφ0 =N+ − N−,
N+ + N−φ0 = 0,(4)

respectively. Notice that φ0 > 0 due to the constraint mentioned above. The configuration of the vesicle at time t is
fully specified using the three dimensionless variables, namely,u(θ,ϕ,t), φ  (θ,ϕ,t), and φ (θ,ϕ,t).
B. Spherical geometry

Using the formalism of differential geometry, we briefly
introduce the definitions of the relevant quantities [37,38]. Letêr , êθ , and êϕ be the orthonormal unit vectors in the sphericalcoordinates. The surface of the vesicle is described using avector function r(θ,ϕ), written asr = r0 [1 + u(θ,ϕ)]êr ,(5)where u is defined in Eq. (1).
The two-dimensional (2D) metric tensor of the surface gij(i,j = θ,ϕ) is defined as gij = ∂i r · ∂j r, where the notation ∂iimplies a partial derivative with respect to the coordinate i.
The inverse of the metric tensor g ij is defined by g ij gj k = δki ,where δki is the Kronecker delta. Hereafter we use Einstein’ssummation convention and sum over all repeated indices. The
unit normal vector pointing outward on the surface is given by√n̂ = (∂θ r × ∂ϕ r)/ g, where g = Det gij . The curvature tensor
that quantifies the extrinsic curvature of the surface is further
defined as Lij = ∂i ∂j r · n̂. Two coordinate-independent scalar
curvatures are obtained from Lij : the mean curvature H =g ij Lij /2, and the Gaussian curvature K = (Det Lij )/g.
√
For a small shape deformation u, we express g and H up
to the second order as(∇⊥ u)2√,(6)
g ≈ r02 sin θ 1 + 2u + u2 +
2


1
∇2 u
H ≈−
1 − u − ⊥ + u2 + u∇⊥2 u ,(7)r0
2
respectively, where we have used
 2
 2
∂u
∂u
1
(∇⊥ u)2 =+,2∂θsin θ ∂ϕ1 ∂ 2u∂ 2 u cos θ ∂u+∇⊥2 u = 2 +.
∂θsin θ ∂θsin2 θ ∂ϕ 2

(8)
(9)

Notice that ∇⊥ is a 2D nabla operator defined on the unit
sphere.

PHYSICAL REVIEW E 94, 062414 (2016)

where σ is the surface tension, κ the bending rigidity, λ thecurvature-density coupling parameter, k the compression modulus, and P = P + − P − the pressure difference between
the outside and the inside of the vesicle in equilibrium. The
first term represents the energy associated with the change inthe total area of the vesicle. The second term describes theenergy cost due to the bending deformation [39]. The third
term represents a coupling between the mean curvature H andthe local density difference φ  . The energy associated withcompression of the bilayer is described by the fourth term. The
first integral is performed over the whole area A of the vesicle,√
and hence the areal element is given by dA = dθ dϕ g. The
volume integral includes the entire volume V of the deformed
vesicle. We assume that the membrane is impermeable, and thefluids inside and outside are incompressible. This means thatthe volume of the vesicle is fixed, which is taken into accountby the last term in Eq. (10). We note that exactly the same freeenergy as given by Eq. (10) was previously employed by Miaoet al. [19].
The physical meaning of the coupling parameter λ can beunderstood such as for a lipid bilayer of finite thickness as
discussed by Seifert and Langer for planar membranes [26].
Since the lipid densities are defined on the bilayer midsurface[see Eq. (2)], the local curvature and density are necessarily
coupled to each other. In Ref. [26], they showed that thestrength of this coupling is proportional to the distancee between the bilayer midsurface and the neutral surface,
and the coupling parameter can be interpreted as λ = ke.
Moreover, the bending rigidity κ in Eq. (10) corresponds to
the renormalized bending rigidity κ̃ = κ + 2ke2 in Ref. [26].
The parameter λ can be controlled either by changing
the length of the lipid tails, or by attaching a layer of othermolecules to the membrane. However, this mechanism alone
cannot significantly affect the λ value. On the other hand,recent work has shown that an asymmetry in the nature of
fluids surrounding the bilayer can also influence the curvaturedensity coupling, even in the case of a single-component
membrane [40]. Moreover, spontaneous curvature can begenerally taken into account using a term similar to thecoupling term in Eq. (10) [41,42]. In such a case, a strong
coupling between the curvature and the membrane internal
degree of freedom causes a curvature instability. Since λ canbe varied by multiple mechanisms, we consider it to be a freeparameter in our study.
A bilinear coupling between the curvature H and φ is
not allowed because the assignment of ± label is arbitrary. In
accordance with the Gauss-Bonnet theorem, we neglect theGaussian curvature term K in the free energy because it isrelevant only when the vesicle topology is allowed to change.

C. Vesicle free energy

The thermodynamic free energy of a compressible vesicle
includes a curvature elastic energy and a stretching elastic
energy, and can be written as [19]


κ
F =
dA σ + (2H )2 + 2λH φ 
2
 
 2

2
+ k (φ ) + (φ )+ dV P ,(10)

III. STABILITY ANALYSIS
A. Free energy variations

Here we investigate the stability of the vesicle by considering small deviations of its configuration from the reference
state. A spherical vesicle of radius r0 with a lipid density
difference φ0 and a density sum φ0 = 0 is the reference
state. As defined in Eq. (1), the quantity u is a first-orderdeviation from the reference radius. For the density variables,

062414-3

SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)

the deviations are given byδφ = φ− φ0 ,

δφ = φ ,(11)

respectively. In terms of these deviations, the zeroth orderterms of the free energy in Eq. (10) describe the ground stateenergy of the vesicle,
4
F (0) = π r03 P + 4π r02 σ + 8π κ3+ 4π r02 k φ0

2− 8π r0 λφ0 .

This second-order free energy will be used to analyze thevesicle stability.
B. Expansions in spherical harmonics

Since the variables specifying the configuration are allfunctions of the angular coordinates θ and ϕ, we expand themin terms of surface spherical harmonics Ynm (θ,ϕ) as(12)u(θ,ϕ,t) =Next, the first-order terms in the free energy are


2λ 
(1)
 2
F = 2k φ0 −φ + 2σ + P r0
dA0 u
r0 0


2λ
+ 2kφ0 −(13)
dA0 δφ  ,r0where the areal element of the reference sphere is defined
here by dA0 = r02 sin θ dθ dϕ and should be distinguishedfrom dA. The condition for the vesicle to be in mechanicalequilibrium is that the first-order variation should vanish, i.e.,
F (1) = 0. For this to hold for an arbitrary deviation u, we
impose the condition
2λ 
2
2k φ0 −φ + 2σ + P r0 = 0.
(14)
r0 0
This relation is a special case of the more general stability
condition of a vesicle when λ is interpreted as a spontaneous
curvature [38]. In the absence of the first two terms, this is
the Young-Laplace equation for surface tension dominatedinterfaces [43]. The second term in Eq. (13) vanishes becauseof the assumption that the total number of lipids is conservedin each monolayer separately, namely, dA0 φ  = A0 φ0 .
After eliminating P in the free energy with the use of
Eq. (14), we obtain the second-order terms in free energy as


κ
1
(2)
F =
dA0 2 (2u∇⊥2 u + (∇⊥2 u)2 )
2
r0


2λ 
2
φ0 (2u2 + u∇⊥2 u)− σ + k φ0 −
r0


4λ
2λ  2
δφ  u ++ 8kφ0 −δφ ∇⊥ ur0
r0

+ 2k (δφ  )2 + (δφ )2 .
(15)

δφ  (θ,ϕ,t) =

λ =

λr0
,
κ

k =

kr02
.
κ


φnm
(t)Ynm (θ,ϕ),(17)

φnm (t)Ynm (θ,ϕ).

(18)



δφ (θ,ϕ,t) =n,mHere, n,m indicates the sum over all spherical harmonic
∞
n

modes
n=0
m=−n , whilen,m in Eqs. (17) and (18)∞ndenotes n=2 m=−n . The n = 0 term is omitted in the latter
two equations because the total number of lipids is conserved.
The n = 1 terms are also omitted as they are not coupled toshape changes for small perturbations [19].
The spherical harmonics are eigenfunctions of the Laplacian operator on the unit sphere ∇⊥2 , and obey the relation∇⊥2 Ynm (θ,ϕ) = −n(n + 1)Ynm (θ,ϕ).

(19)

By using the spherical harmonics representation, the secondorder free energy is written in the bilinear form as⎛F (2) =

unm

⎞∗

κ
⎜ ⎟

φnm A⎝φnm
unm φnm
⎠ ,
2 n,m
φnm

(20)

where ∗ denotes the complex conjugate. We have used
the orthogonality property of spherical harmonics to obtain
Eq. (20). Since the three configurational variables are real,
one has u∗nm = un,−m . The 3 × 3 free energy matrix A is
made dimensionless by dividing all the components by κ. The
symmetric matrix A is given by
4k  φ0 − 4λ − λ (n − 1)(n + 2)

⎞

0

2k 

⎟
0⎟
⎠.

0



2k

(21)

equilibrium fluctuation spectrum for the spherical harmoniccoefficients as

In the above, the dimensionless parameters are defined by
σ r02
,
κ

(16)

n,m

⎛

2
(n − 1)(n + 2) n(n + 1) + σ  − 2λ φ0 + k  φ0
⎜A=⎜
4k  φ0 − 4λ − λ (n − 1)(n + 2)
⎝
0

σ =

unm (t)Ynm (θ,ϕ),
n,m


(22)

The term proportional to |unm |2 in Eq. (20) represents theenergy of the bending mode. Using the equipartition theorem
and taking the limits of λ → 0 and φ0 → 0, we obtain the|unm |2 =1kB T,κ (n − 1)(n + 2)[n(n + 1) + σ  ](23)where kB is the Boltzmann constant and T the temperature [44]. This expression is commonly used to estimate the062414-4

RELAXATION DYNAMICS OF A COMPRESSIBLE BILAYER . . .

PHYSICAL REVIEW E 94, 062414 (2016)

bending rigidity κ and the surface tension σ of RBCs orvesicles in flicker analysis [34,45].
C. Stability of the free energy

The stability of the compressible bilayer vesicle is ensured
by the positivity of the eigenvalues of the free energy matrix.
Since A is a real symmetric matrix, all its eigenvalues are real.
The matrix A is further decomposed into one diagonal element
A33 = 2k  and 2 × 2 sub-matrix a given by
a=


A11
A21


A12
.
A22

(24)

The two eigenvalues ε1 and ε2 (ε1 < ε2 ) of the matrix a canbe calculated in terms of its trace and determinant. Together
with ε3 = 2k  , they constitute the three eigenvalues of thematrix A. The inverse of these three eigenvalues correspondto the “susceptibilities” of a compressible bilayer vesicle, andquantify the proportionality between externally applied forces
and the changes in the three variables u, φ  , and φ . In otherwords, the smaller the susceptibility values are, the larger theforces are that are required to induce the variable changes.
We numerically obtain the two eigenvalues ε1 and ε2 using
the following values for the parameters. The bending rigidityof membranes has been estimated in several experiments tobe about κ ≈ 10−19 J [34,35,45,46]. For diacyl phosphatidylcholine (PC) lipids in their fluid phase, the compression
modulus is estimated to be k ≈ 10−1 J/m2 [35,36]. The rangeof values of the surface tension σ changes over four ordersof magnitude between 10−11 –10−7 J/m2 . By considering atypical vesicle radius r0 = 10 μm, we obtain the parametervalues k  = 108 and σ  = 10−2 [see Eq. (22)].
Figure 2 shows the stability diagram of the free energy
in the parameter space (λ , φ0 ). The white region indicatesthe parameter ranges where all the eigenvalues are positive
and the vesicle is stable. Both the yellow and the orange
regions represent unstable regions where at least one of theeigenvalues is negative. In the yellow regions, the eigenvalue
ε1 for n = 2 mode is negative (and hence unstable), while allother eigenvalues are positive. The condition for this instability
to occur is approximately obtained by requiring that Det aFIG. 2. Stability diagram in the (λ , φ0 ) plane when σ  = 10−2
and k  = 108 (see the text). The vesicle is stable in the white region,
while it is unstable in the yellow and orange regions. In the yellow
region, only the n = 2 mode is unstable, whereas larger-n modes are
unstable in the orange region. We investigate the stability and the
dynamics represented by the three points inside the stable region: (i)(λ ,φ0√) = (104 ,3 × 10−4 ) (black), (ii) (104 ,5.6 × 10−4 ) (blue), and(iii) ( 2 × 104 ,4 × 10−4 ) (green). The case (i) is well inside thestable region, while the cases (ii) and (iii) are very close to the unstableregion.

becomes negative for n = 2:

  3λ 
φ −
<
 0
k 



λ 2
6 + σ
+
.
k 2
k

(25)

In the orange regions, the larger n modes become unstable.
The condition for this instability can be also obtained bynoticing that the sign of ε1 is determined by that of Det a. Since
Det a ≈ (2k  − λ2 )n4 √for large n, the eigenvalue becomes
negative when |λ | > 2k  . Notice that the instability in theorange region is induced only by the parameter λ .
In Fig. 3, we show the two inverse susceptibilities ε1 and ε2
as a function of n for three sets of parameters (i), (ii), and (iii)

FIG. 3. The two eigenvalues ε1 and ε2 of the energy matrix A as a function of the spherical harmonic√ mode n [see Eq. (21)]. The chosenparameters are (a) (λ ,φ0 ) = (104 ,3 × 10−4 ) [case (i) in Fig. 2], (b) (104 ,5.6 × 10−4 ) [case (ii)], and (c) ( 2 × 104 ,4 × 10−4 ) [case (iii)]. Other
parameter values are k  = 108 and σ  = 10−2 as in Fig. 2. The numbers in the plots indicate the slopes of different curves.
062414-5

SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)

marked in the stable region of Fig. 2. The point (i) is locatedwell inside the stable region, while the points (ii) and (iii) aremarginal cases with respect to the two instabilities mentioned
above. We do not plot the third eigenvalue ε3 as it is merely
a constant. Figure 3(a) corresponds to the case (i) where theparameters are λ = 104 and φ0 = 3 × 10−4 . This value of λcorresponds to a membrane with e = 1 nm between the bilayer
midsurface and the neutral surface. According to Eq. (21),the energy associated with the pure bending deformation is
proportional to n4 , which is clearly seen in Fig. 3(a). There isalso a mode crossing behavior between the n4 bending mode
and a constant compression mode around n ≈ 100.
Figure 3(b) shows the eigenvalues for case (ii) where theparameters are λ = 104 and φ0 = 5.6 × 10−4 , which areclose to the yellow region in Fig. 2 where the n = 2 instability
occurs. The eigenvalue ε1 for n = 2 is much smaller than thatin Fig. 3(a). In the range 2 < n < 10, ε1 is slightly larger thanin (a), while the eigenvalues remain almost unchanged forlarger n. Figure 3(c) shows the√eigenvalues for the case (iii)
where the parameters are λ = 2 × 104 and φ0 = 4 × 10−4 ,which are close to the orange region in Fig. 2. Compared to(a), the n dependence of ε1 is significantly altered because weobserve n2 dependence for smaller n, and ε1 even decreases
for larger n values.
The asymptotic value of the smaller eigenvalue is ε1 ≈2k  − λ2 for large n. In all the three cases (a), (b) and (c), thelarger eigenvalue ε2 remains almost the same. It behaves as
102 and ε2 ≈ n4 for n102 .
ε2 ≈ 2k  for n
IV. MEMBRANE HYDRODYNAMICS
A. Onsager’s variational principle

Onsager’s variational principle is a fundamental framework
to describe the dynamical behavior of soft matter [29]. It has
been applied in a variety of problems including hydrodynamics
of liquid crystals, diffusion, and sedimentation of colloidal
particles [30]. In the context of Stokesian hydrodynamics,
Onsager’s principle states that the viscous forces are balancedby the potential forces. In this approach, a “Rayleighian”
function is defined as the sum of dissipation function  andthe time derivative of the free energy as
dF
.
(26)
dt
Note that  is half the free energy dissipated into heat per unit
time. The dynamical equations are obtained by extremalizing
R with respect to all the dynamical variables. Any constraint
is taken into account by introducing an additional termwith a corresponding Lagrange multiplier. A more detailed
description on Onsager’s variational principle is given in
Refs. [28–30].
For a compressible bilayer vesicle, we need to obtain the
dynamical equations for the three variables defined in Sec. II,
i.e., u, φ  , and φ . As we have assumed before, the fluid
medium on either side of the membrane is incompressible.
However, the inside fluid viscosity η− is taken to be different
from the outside viscosity η+ . The components of the fluid
velocity in these regions are represented by Vα± , and thehydrodynamic pressure fields are denoted by P ± . Hereafter,
greek indices are used for components (r,θ,ϕ) of vectors in
R=+three-dimensional (3D) space. The two lipid monolayers are
regarded as ideal 2D fluids, and their velocity fields are denotedby vi± .
We take into account two dissipation mechanisms: thedissipation due to the bulk fluid viscosity and the dissipation
due to the inter-monolayer friction. The dissipation functionals
for the outer and inner fluid medium are [28],

±
±
±
 =
dV (η± Gσ α Gρβ DαβDσρ− P ± Gσ α ∂α Vσ± ), (27)±±where Dαβ= (∂α Vβ± + ∂β Vα± )/2 is the rate of deformation
αβtensor in the bulk fluid and G is the inverse of 3D (rather
than 2D) metric tensor Gαβ . The + (−) subscript of the integral
indicates that it is performed over the volume of the outside(inside) bulk fluid. The dissipation functional due to the intermonolayer friction is given by
w
w = dA0 (vi+ − vi− )2 ,(28)
2where w is the friction constant.
Next we discuss the constraints and the hydrodynamic
boundary conditions. Since the monolayer lipid numbers N ±are conserved, the continuity equation should hold in eachmonolayer separately:∂ρ ±ρ ± ∂g+ g ij ∂j (ρ ± vi± ) +
= 0.
∂t
2g ∂t

(29)

The last term exists due to the curved geometry of the vesicle,
and it describes the change in the lipid density when thelocal area element is varied [23,47]. After the linearizationof Eq. (29), we obtain the following two equations in terms ofsmall variables:∂  1 ijδφ + (g ∂j vi+ − g ij ∂j vi− ) = 0,
∂t
2
∂
1 ij
∂u
δφ + (g ∂j vi+ + g ij ∂j vi− ) + 2
= 0.
∂t
2
∂t

(30)
(31)

The no-slip boundary condition requires that the velocityof the fluid coincide with that of the vesicle at the interface.
Since u1, we approximate that this interface lies at r =r0 . While equating the velocities, the correction due to thefinite membrane thickness can be neglected because e/r0u.
Hence the no-slip condition becomes
vi± = Vi± ,

r0

∂u
±
= V(r)
.
∂t

(32)

±
is the physical component of the fluid velocity
Here V(r)along the radial direction, and should be distinguished from
Vr± . Physical components are the components of a tensor
quantity with a correct dimension [37]. In orthogonal coordinate systems, the physical components are related tothe covariant components of a vector through the relationhα V(α) = Vα . The scale factors hα for the spherical coordinatesare (hr ,hθ ,hϕ ) = (1,r,r sin θ ). Additionally, unlike the tensor
indices, the quantities in the bracket are not summed over.
Using Eqs. (10), (27), (28), (30), (31), and (32), the total
Rayleighian functional for the vesicle and the surrounding

062414-6

RELAXATION DYNAMICS OF A COMPRESSIBLE BILAYER . . .

PHYSICAL REVIEW E 94, 062414 (2016)

by extremalizing R with respect to ∂t u:

fluid is given bydF
R = + +  − +  w +
dt




∂u
 ij




− V(r)+
dA0 μi g vj − Vj + ξ r0
∂t
=±
 


∂  1 ij
+ dA0 ζ 
δφ + [g ∂j vi+ − g ij ∂j vi− ]
∂t
2


∂
1
∂u
δφ + [g ij ∂j vi+ + g ij ∂j vi− ] + 2
.
+ζ
∂t
2
∂t
(33)

r0 (ξ + + ξ − ) + 2ζ − k φ0 (2u + ∇⊥2 u) + 4kφ0 δφ 
2

2λ 
λ
2λ  24λ 
δφ + ∇⊥2 δφ  +φ 0 ∇⊥ u +
φ u
r0
r0
r0
r0 0
κ
+ 2 (2∇⊥2 u + ∇⊥4 u) − σ (2u + ∇⊥2 u) = 0.
r0−

Finally, we extremalize R with respect to vi± and obtain thelateral force balance on each of the monolayers as
1
±
∓

μ±
i + w(vi − vi ) ∓ ∂i (ζ ± ζ ) = 0.
2

Here we have introduced the Lagrange multipliers μiandξ to take into account the boundary conditions in Eq. (32).
The other Lagrange multipliers ζ  and ζ are used for thecontinuity equations in Eqs. (30) and (31), respectively. The
full expression of the Rayleighian is given in the Appendix A.
B. Basic equations

Using the Rayleighian in Eq. (33), we derive the force
balance conditions and hence the dynamical equations byextremalizing it with respect to Vα± , P ± , ∂t (δφ  ), ∂t (δφ ),∂t u, and vi± , where ∂t indicates the time derivative.
Extremalizing Eq. (33) with respect to Vα± yields
±

−η G

σα

∂σ ∂α Vβ± + ∂β P ± = 0,Gσ α ∂α Vσ± = 0.

Next we use Lamb’s solution for Stokes equation in Eq. (34)in spherical coordinates [23,48,49]. It is represented in termsof the pressure field P ± (r,t) and the newly introduced function
ψ ± (r,t) that is the solution of the homogeneous Stokes
equation Gσ α ∂σ ∂α Vβ± = 0. We expand these functions in termsof solid spherical harmonics as


∓η

1
±
±
+ ∂(θ) V(r)−∂(r) V(θ)
r0

±
V(θ)

ψ (r,t) =P (r,t) =±∓2η± ∂(r) V(r)± P ± − ξ ± = 0.

ζ + 2kδφ = 0,



r
r0
r
r0

n
Ynm (θ,ϕ),(42)

Ynm (θ,ϕ).

(43)

n





(n − 2)
P + r 2∇+ n(2n − 1) nm
2η
n,m
 
(n + 1)
r0 n+1++ +Pnm r
Ynm (θ,ϕ), (44)
η n(2n − 1)r (n + 3)
−
− 2
Pnm
ψnm
V− (r) =
∇+ −
r ∇
2η
(n
+
1)(2n
+
3)
n,m
 n
n
r
−
− −
Pnm
r
Ynm (θ,ϕ).
η (n + 1)(2n + 3)
r0(45)
V+ (r) =

(37)



4λ
λ
ζ  + 2kδφ  + 4kφ0 −u + (2u + ∇⊥2 u) = 0,
r0
r0
(38)

±
Pnm
(t)



Then Lamb’s solution can be written as(36)

In these equations, we have also included the surface terms
which arise from the volume integral in R after performing
the integration by parts. Equation (36) is obtained by settingi = θ after the extremalization. Although a similar equation
for the coordinate ϕ can be also obtained, it is unnecessary forthe present calculation. Equations (36) and (37) are later used±
to eliminate the Lagrange multipliers μ±i and ξ .
Extremalizing R with respect to ∂t (δφ ) and ∂t (δφ ), we
obtain

±
ψnm
(t)

n,m

(35)

− μ±
(θ) = 0,

n,m


±



r0



±

Extremalizing R with respect to Vi± leads to the following
equations:±(41)

C. Solutions of hydrodynamic equations(34)

which is the Stokes equation for the bulk fluid inside andoutside the vesicle. By extremalizing R with respect to P ± ,we recover the incompressibility condition of the bulk fluid

(40)

+
ψnm
∇−

The velocity field of the lipid monolayer v± is expanded interms of vector spherical harmonics. Here, only the component
along ∇Ynm is relevant for our calculation, and it is denoted as±. Using the boundary conditions in Eq. (32), we equate thevnm
components along ∇Ynm to obtain
1(n − 2)
+
+
+ 2
ψnm
(46)
Pnm
vnm
=
− +
r0 ,
r0
2η n(2n − 1)1(n + 3)
−
−
− 2
vnm
=
+ −
r0 . (47)
ψnm
Pnm
r0
2η (n + 1)(2n + 3)

(39)

respectively, which are also used to eliminate ζ  and ζ . The
force balance condition along the radial direction is obtained

These equations are used in the force balance conditions
Eqs. (40) and (41) to derive the final dynamical equations
as derived in the Appendix B.

062414-7

SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)

D. Dynamical equations

V. RELAXATION MODESAs shown in the Appendix B, the dynamical equation for
the variable u is given by

A. Three relaxation modesn −n∂unm = 2 ψnm
P− ,
+ −
∂t
2η (2n + 3) nm
r0

(48)

[see Eq. (B1)]. The dynamical equations for the density
variables φ  and φ are derived from Eqs. (30) and (31)as
(n − 2)r02∂ n(n + 1)+ψnm
δφnm =P+− +2∂t2η n(2n − 1) nm
2r0

(n + 3)r02
−
−
,
Pnm
− ψnm
− −
2η (n + 1)(2n + 3)∂
n(n + 1) +(n + 1)(n − 2) +ψnm −δφnm =
P2∂t
4η+ (2n − 1) nm
2r0
 −

ψ
n(n + 1)− 2n nm
+
2
r02


n
n(n + 3)−P− ,+
4η− (2n + 3) η− (2n + 3) nm

(49)

(50)

respectively.
±±
and Pnm, as obtained in the
By using the solutions for ψnm
Appendix B, the final dynamical equations are expressed in
the matrix form as
⎛
⎛
⎞
⎞
unm
unm
∂⎜ ⎟
1 ⎜ ⎟(51)⎝φnm ⎠ = − B⎝φnm
⎠.
∂t
τ
φnm
φnm

The coupled equations for the relaxation dynamics of a
weakly compressible bilayer vesicle using the three variablesu, φ  , and φ are obtained in Eq. (51). The three eigenvaluesof the dynamical matrix B, denoted as γ1 < γ2 < γ3 , givethe relaxation rates as a function of n. The presence of theoff-diagonal elements in B indicates that these variables arecoupled to each other. This is a fundamental and unique featureof a spherical vesicle, because the mode associated with φis always decoupled from the other two modes for planar
bilayers [26].
As for the dynamical parameters, the inter-monolayer
friction is estimated to be about w ≈ 109 N s/m3 according
to the vesicle fluctuation analysis [50]. The viscosity of theexternal fluid is set to be η+ = 10−3 Pa s for water. Hence,in addition to the dimensionless static parameters given inSec. III, we use here w  = 107 . With these parameter values,
the characteristic time scale in Eq. (52) becomes τ = 10 s. We
first consider the equal viscosity case E = 1 and numericallycalculate the three eigenvalues in Fig. 4 when λ = 104 and
φ0 = 3 × 10−4 [case (i) in Fig. 2]. Notice that all the relaxation
rates are made dimensionless here by using τ . In order toobserve the mode crossing behavior, we have computed up ton = 104 which corresponds to few nm wavelength for a vesicleof size r0 = 10 μm.
From Fig. 4, we find that the two eigenvalues γ1 and γ2
are asymptotically proportional to n2 and n3 , respectively, forlarge n. In the later subsection, we show that the n2 dependence
predominantly corresponds to the slipping relaxation mode,
while the n3 dependence is associated with the bending

Here the 3 × 3 matrix B is made dimensionless by using thecharacteristic time scaleτ=η+ r03.
κ(52)

We also define the viscosity contrast E as the ratio between
the viscosities inside and outside of the vesicle
E=

η−
.
η+

(53)

We further introduce a dimensionless friction coefficient
given by
w =

wr0
.
η+

(54)

The dynamical matrix B in Eq. (51) is expressed as a productof two symmetric matrices
B = C · A,

(55)

where A was defined before in Eq. (21) as a free energy matrix.
On the other hand, the matrix C depends only on E and w which characterize the dissipation of the whole system. Since
the full analytical expression for the matrix C is lengthy, we
present each of its components in Appendix C, which provides
the main result of our analytic calculation.

FIG. 4. The three relaxation rates γ1 , γ2 , and γ3 as a function of thespherical harmonic mode n obtained from the dynamical matrix B [seeEq. (55)]. The chosen parameter values are (λ ,φ0 ) = (104 ,3 × 10−4 )[case (i) in Fig. 2], k  = 108 , σ  = 10−2 , and w = 107 . The fluid
viscosities are taken to be symmetric between inside and outside of
the vesicle, i.e., E = η− /η+ = 1. All the relaxation rates are scaled by
the characteristic time τ = η+ r03 /κ. The numbers indicate the slopesof different curves.

062414-8

RELAXATION DYNAMICS OF A COMPRESSIBLE BILAYER . . .

relaxation mode. At around n ≈ 20, these two modes start
to overlap each other, which is a mode crossing behavior. The
presence of the slipping mode and the crossover between thetwo relaxation modes were first considered by Seifert andLanger for planar bilayers [26]. Indeed neutron spin echo
and flicker spectroscopy studies of GUVs have confirmed thatlarge-n relaxation is dominated by the slipping mode [20,51].
The three relaxation rates for n = 2 mode are estimated tobe 0.296 s−1 , 3.097 s−1 , and 1.257 × 107 s−1 . The smallest
relaxation rate obtained here is close to that measured using
flicker analysis [45]. We also notice that the largest eigenvalue
γ3 is proportional to n throughout the range of n studied, andcorresponds to the relaxation of the total density φ [19,26].
This mode is much (about eight orders of magnitude) fasterthan the other two slower ones.

PHYSICAL REVIEW E 94, 062414 (2016)

to eliminate it from the relaxation equations. Hence weassume that the variable φ equilibrates immediately, and set∂φ /∂t ≈ 0. Within this approximation, we can eliminate φfrom the coupled dynamical equations Eq. (51) to obtain




∂ unm
1 unm
=− b  .
∂t φnmτ φnm(56)Here the 2 × 2 matrix b can be again written as a product of
two symmetric matrices as
b = c · a,(57)

B. Elimination of the fastest mode

Since the relaxation of the total density is much fastercompared to the other two slower modes, we are allowed n(n + 1) 4w + (2n + 1)(1 + E)c=d−(n + 2) + (n − 1)E−(n + 2) + (n − 1)E(n + 2)(2n2 − n + 2) + (n − 1)(2n2 + 5n + 5)E

In the above, the common denominator d is given by
d = 4w  [(n − 1)(2n2 + 5n + 5)E + (n + 2)(2n2 − n + 2)]+ 2(n − 1)(n + 1)(2n2 + 4n + 3)E 2+ [8n2 (n + 1)2 − 5]E + 2n(n + 2)(2n2 + 1).

where a was defined in Eq. (24), and the matrix c takes the
form

(59)

In Fig. 5, we plot the two eigenvalues of the matrix b (black
symbols) for the same parameter values used in Fig. 4. The tworelaxation rates obtained here completely coincide with γ1 andγ2 in Fig. 4. Such a coincidence justifies our approximation toeliminate the fastest relaxation mode associated with φ .
C. Effects of different viscosities

.

(58)

As shown in Fig. 5, the viscosity contrast shifts the modecrossing point between the slipping and the bending modes.
To better understand this effect, we approximate one of theeigenvalues of the matrix b = c · a as follows. Keeping only
the most dominant term, one can approximate Eqs. (58), (59),and (24) asn(n + 1) 4w 
0
,
(60)
c≈
0
2n3 (1 + E)dd ≈ 8w n3 (1 + E), 4−λ n2n.
a≈−λ n2
2k 

(61)
(62)

Then the larger eigenvalue of b is approximated as

In this subsection, we discuss the effect of viscosity contrast
E = η− /η+ on the relaxation dynamics of a compressible
bilayer vesicle. Figure 5(a) shows the relaxation rates when
the inner viscosity is 10 times larger than that of the outside,i.e., E = 10 (red symbols), while all the other parameters
are same as the case of E = 1 (black symbols). The modecrossing behavior between the bending and the slipping modes
at around n ≈ 200 is more apparent in the higher viscosity
contrast case. When E = 10, the bending mode becomesslower than that of E = 1 case, while the slipping moderemains almost unaffected. The effect of viscosity contrast
is more pronounced in the case of E = 100 (red symbols)
as shown in Fig. 5(b). Here the slowest relaxation modeis dominated by the bending one up to n ≈ 2000. Such arelaxation behavior is significantly different from the equal
viscosity case, where the slipping mode dominates most of therelevant n range. In Fig. 5(b), we find that the relaxation rate γ1for small n values differs by two orders of magnitude between
the cases E = 1 and 100.

w  n3 + k  n2 (1 + E)
2w  (1 + E)
⎧  2
kn
⎪
⎪
⎨ ,
n
n∗ ,
2w
≈
n3
⎪
⎪
⎩
, n
n∗ ,
2(1 + E)

γ2 ≈

where

(63)

k
(1 + E)(64)w
is the crossover mode and is linearly proportional to theviscosity contrast E.
The above limiting expressions clearly show that the n2
dependence corresponds to the slipping mode, while the n3
dependence reflects the bending mode. Moreover, only thebending mode is affected by the viscosity contrast E as it
appears only for nn∗ in Eq. (63). The above limiting
expressions of the relaxation rates coincide with those for062414-9n∗ =SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)

FIG. 5. The two slower relaxation rates γ1 and γ2 as a function of
the spherical harmonic mode n obtained from the dynamical matrix b[see Eq. (57)]. The chosen parameter values are (λ ,φ0 ) = (104 ,3 ×
10−4 ) [case (i) in Fig. 2], k  = 108 , σ  = 10−2 , and w = 107 . (a)The comparison between E = η− /η+ = 1 (black) and E = 10 (red).
(b) The comparison between E = 1 (black) and E = 100 (red). All
the relaxation rates are scaled by the characteristic time τ = η+ r03 /κ.
Arrows indicate the crossover mode n∗ given by Eq. (64).

a planar bilayer [26] when E = 1, and n/r0 is regarded as awave number in the limit of r0 → ∞.
D. Dynamics close to the unstable regions

As we have discussed in Sec. III, the two distinct instabilities are identified in the (λ ,φ0 ) parameter plane. Here
we focus on the relaxation dynamics of a vesicle when theparameters are close to these unstable regions as marked by
the points (ii) and (iii) in Fig. 2. In Fig. 6(a), we plot the tworelaxation rates from the matrix b as a function of n when
the parameter values are that of point (ii), i.e., λ = 104 and
φ0 = 5.6 × 10−4 . Compared to the case (i), the relaxation rate
γ1 with n = 2 is significantly reduced. The relaxation rate γ1

FIG. 6. The two slower relaxation rates γ1 and γ2 as a function of
the spherical harmonic mode n obtained from the dynamical matrix b[see Eq. (57)]. The chosen parameter values are k  = 108 , σ  = 10−2 ,w = 107 , and E = η− /η+ = 1. (a) The comparison between thestable case (i) (black) and the marginal case (ii) (blue) in Fig. 2. (b)The comparison between the stable case (i) (black) and the marginalcase (iii) (green) in Fig. 2. All the relaxation rates are scaled by thecharacteristic time τ = η+ r03 /κ.

in the range 2 < n < 10 is slightly larger than that of case (i),and they are almost the same for larger n. The relaxation rates
γ2 does not vary between the cases (i) and (ii), except for asmall increase in the case of (ii) when n = 2.
In Fig. 6(b), we plot the two relaxation rates from the matrixbwhen
the parameter values are that of the case (iii), i.e., λ =√2 × 104 and φ0 = 4 × 10−4 . As compared to the case ofpoint (i), the large-n relaxation of the slipping mode is slowed
down by more than four orders of magnitude. Moreover, theslipping mode is almost independent of n in the intermediate-n
range (150 < n < 500) for the case (iii). The relaxation rate
γ2 is not significantly different between the cases (i) and (iii).

062414-10

RELAXATION DYNAMICS OF A COMPRESSIBLE BILAYER . . .
E. Limiting expressions

Here we make a connection of our work to that by Seki
and Komura [23], where they investigated the relaxation
dynamics of a compressible vesicle without considering abilayer structure. In the present work, this corresponds tothe case in which the friction between the monolayers islarge enough such that the two monolayers move together
as a single layer, i.e., w  → ∞. Moreover, we also take thelimit of λ → 0, φ0 → 0 in order to recover the previous
result.

PHYSICAL REVIEW E 94, 062414 (2016)

In this limit, all the coefficients of φ  in the dynamical equationin Eq. (51) vanish, resulting in a reduced 2 × 2 dynamical
matrix:




∂ unm
1 unm
= − b̃
.
(65)∂t φnmτ φnm
Again, the matrix b̃ can be expressed as a product of twosymmetric matricesb̃ = c̃ · ã.

(66)

Here the matrix c̃ is
(2n + 1)(1 + E)n(n + 1)c̃ =d̃−(n + 2) + (n − 1)E−(n + 2) + (n − 1)E(n + 2)(2n2 − n + 2) + (n − 1)(2n2 + 5n + 5)E

with a common denominatord̃ = [2(n − 1)(n + 1)E + 2n + 1]
2

× [(2n2 + 4n + 3)E + 2n(n + 2)].
On the other hand, the matrix ã is given by(n − 1)(n + 2)[n(n + 1) + σ  ]
ã =
0

0

(68)



2k 

.

(69)

Following the work of Ref. [23], we obtain two limiting
expressions for the relaxation rates. When the membranecompressibility is large compared to the bending rigidity, i.e.,n(n + 1) + σ k  , we getγn =[n(n + 1) + σ  ](n − 1)n(n + 1)(n + 2).
(n − 1)(2n2 + 5n + 5)E + (n + 2)(2n2 − n + 2)
(70)

If we further set E = 1, this expression reduces to that obtained
first by Milner and Safran [17]. When the bending rigidity
dominates over the compressibility k n(n + 1) + σ  , weobtainγn =2k  n(n + 1).
(2n + 1)(1 + E)

(71)

The limiting expressions Eqs. (70) and (71) exactly coincidewith Eqs. (D.6) and (D.7) in Ref. [23], respectively.
VI. SUMMARY AND DISCUSSIONSIn this paper, we have investigated the statics and dynamics
of a weakly compressible bilayer vesicle. First we have calculated a free energy matrix in terms of linear perturbations inlocal curvature and local densities. Calculating the eigenvaluesof the free energy matrix, we performed the stability analysis
by varying the curvature-density coupling parameter λ , andthe lipid density difference between the two monolayers φ0 .
As a result, we identified two different instabilities: one
affecting the small-n modes and the other influencing thelarge-n modes.
Onsager’s variational principle offers an universal framework for obtaining the dynamical equations for soft matter,


,

(67)

and it has been used to derive the dynamical equations of acompressible bilayer vesicle. The eigenvalues of the dynamical
matrix B in Eq. (55) give the relaxation rates for a vesicle
whose inside and outside fluids are characterized by different
viscosities. The three relaxation modes are coupled to eachother as a consequence of the bilayer and the spherical structureof the vesicle. Assuming that one of the relaxation modes ismuch faster than the other two, we derived the dynamical
equation for the slower modes in Eq. (56). We focused on theeffect of viscosity contrast E = η− /η+ on the relaxation rates,
and found that it linearly shifts the crossover n mode between
the bending and the slipping relaxations [see Eq. (64)]. As E isincreased, the relaxation rate of the bending mode decreases,
while that of the slipping mode remains almost unaffected.
For parameter values close to the unstable region, some of
the relaxation modes are dramatically reduced. For example,
as they approach the region of large-n instability, we find an
unusual behavior of the relaxation rate as shown in Fig. 6(b).
Although the viscosity contrast in vesicles and RBCs has
been discussed in some experiments [11,16], we have derived
here the exact relaxation rates of a compressible bilayer
vesicle for arbitrary E values. In experiments, a commonly
encountered case is where the inside viscosity is slightly
larger than the outside E > 1, for which the bending mode isslowed down and crossover mode n∗ becomes larger. For evenlarger viscosity contrast (E100), the relaxation is entirely
dominated by the bending mode. For E < 1, on the other hand,n∗ does not depend on E and the relaxation rate of the bending
mode is slightly increased. Our result clearly shows that theviscosity contrast significantly affects the dynamical behaviorof a vesicle.
In the present work, the vesicle stability was analyzed interms of the curvature-density coupling parameter λ and thelipid density difference between the monolayers φ0 . Whereas
the large-n instability is induced only by λ , the small-n
instability can be triggered by changing either λ or φ0 .
Since the large-n instability corresponds to small wavelength
deformations, it leads to the stabilization of small buds or
tubelike deformations. For parameter ranges close to theyellow region in Fig. 2, we predict that the small-n modes
slow down significantly. Although it would be experimentally

062414-11

SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)

challenging, the control of the number of lipids in eachmonolayers enables the change of the parameter φ0 .
Here we shall briefly mention the relation of our result to theprevious theoretical works. Extending the work by Schneider
et al. [18], Milner and Safran derived the bending relaxation
rate in vesicles and microemulsion droplets [17]. In their work,
however, the membrane was assumed to be an incompressible
2D sheet immersed in fluids having the same viscosity on eitherside of the membrane. This case is equivalent to a bilayerthat moves together as a single entity with no difference inthe number of molecules in the upper and lower leaflets, i.e.,φ0 → 0. The bending relaxation rate obtained by Milner andSafran corresponds to a limiting case of Eq. (70) when E = 1.
The slipping relaxation rates for planar bilayers was originallydiscussed by Seifert and Langer who included the dissipation
due to inter-monolayer friction [26]. As mentioned before,
all the results obtained for the planar membrane case can be
reproduced from our results by setting n = qr0 (q being a wave
number) and taking the limit of r0 → ∞.
It is worthwhile to mention that, even though we havediscussed only the case when the vesicle size is r0 = 10 μm,
our results can be readily used to predict the relaxation
behavior of vesicles of other sizes, provided that the parametersk  , σ  , λ , and w  are scaled appropriately. For smaller vesicles
such as r0 < 10 μm, the physically relevant n range is reduced,
and even small changes in viscosity contrast significantly
affects the relaxation behavior. Vesicles in biological systems
usually belong to this size range.
The effects of surface tension on the relaxation dynamics
of membranes were investigated in some of the previous

works [27,28]. In general, the surface tension makes the smalln relaxation of the bending mode faster. Arroyo and DeSimone
studied the dynamics of fluid membranes with general curved
shapes, not restricted to a quasispherical vesicle, and also amembrane with free or internal boundaries [21]. They took into
account the 2D viscosity of the membrane monolayers, which
has been neglected in our work. Based on the previous results,
we expect that including the membrane 2D viscosity will slow
down only the large-n slipping relaxation modes [21,26].
As a future work, it is interesting to consider the situation
where the internal material is a viscoelastic fluid as investigatedin the experiment [52]. Although, the viscoelasticity of themembrane itself has been taken into account in some of theprevious works [53,54], the viscoelasticity of the fluid inside
has not yet been studied. One can expect that the existence ofdifferent timescales due to the viscoelasticity of the inner fluid
will make the dynamics of the vesicle much richer [55–57].

ACKNOWLEDGMENTS

T.V.S.K. thanks Tokyo Metropolitan University for the
support provided through the co-tutorial program. S.K. acknowledges support from the Grant-in-Aid for Scientific
Research on Innovative Areas “Fluctuation and Structure”(Grant No.25103010) from the Ministry of Education, Culture,
Sports, Science, and Technology of Japan, and from the Grantin-Aid for Scientific Research (C) (Grant No.15K05250) fromthe Japan Society for the Promotion of Science (JSPS).

APPENDIX A: RAYLEIGHIAN FUNCTIONAL

The complete expression for the Rayleighian functional used to derive the basic equations in Sec. IV B is given by

R=
=± 

=±




+

dA0 ζ 



+







dV η Gσ α Gρβ DαβDσρ− P  Gσ α ∂α Vσ +

dA0





w +
∂u
− 2
ij  
− V(r) + (vi − vi )dA0 g μi vj − Vj + ξ r0
∂t
2



 
∂  1 ij∂1∂u
δφ + (g ∂j vi+ − g ij ∂j vi− ) + dA0 ζ
δφ + (g ij ∂j vi+ + g ij ∂j vi− ) + 2
∂t
2
∂t
2
∂t
2

− k φ0 (2u + ∇⊥2 u) + 4kφ0 δφ  +λ− σ (2u + ∇⊥2 u) + ∇⊥2 δφ 
r0κ2λ  24λ 
2λ (2∇⊥2 u + ∇⊥4 u) +φ 0 ∇⊥ u +φ0 u −
δφ
2
r0
r0
r0
r0




∂u
λ
∂ 
∂


2
+ 2kδφ + 4kφ0 u − (2u − ∇⊥ u)
δφ + (2kδφ ) δφ .
∂t
r0
∂t
∂t

APPENDIX B: FORCE BALANCE EQUATIONSIn this appendix, we derive a set of four linear equations toeliminate the variables introduced in the solution of Stokes
equations. These are the force balance equations and theboundary conditions at the membrane surface. We expand all
the variables in terms of the spherical harmonics, and equate
the components of Ynm . We use Eqs. (36)–(39) for eliminating
the Lagrange multipliers.

(A1)

Using the boundary condition Eq. (32), we canwrite

062414-12−(n + 1) +n+1+
r0 Pnmψnm + +
r02η (2n − 1)
∂
n −
n
−
r0 Pnm= ψnm+ −= r0 unm .
r02η (2n + 3)∂t
(B1)

RELAXATION DYNAMICS OF A COMPRESSIBLE BILAYER . . .

The force balance along the radial direction in Eq. (40)expressed in the spherical harmonics is
2η+ (n + 1)(n + 2) +2η− n(n − 1) −
−
ψnm +
ψnm
r0
r0




n(n + 1)n(n + 1)+−
+ 1 r0 Pnm− 1 r0 Pnm++
2n − 12n + 3
2λλ+ 4kφ0 −− n(n + 1) δφnm
− 4kδφnm
r0
r0


κ
2λ 
2
+ k φ0 + 2 n(n + 1) −φ0 + σr0
r0
× [n(n + 1) − 2]unm = 0.

and(n + 3)wr0
P − = 0,
η− 2(n + 1)(2n + 3) nm

Here, we list the elements of the symmetric matrix C
defined in Eq. (55).
2n(n + 1)(2n + 1)[w  (1 + E) + (2n + 1)E], (C1)D
3n(n + 1)
C12 = C21 = −(2n + 1)E,(C2)
DC13 = C31

(B2)

(B3)



k
η−
wr0

−
δφnm − δφnm
+ 2 2n − 2 + − ψnm
r0
η
r0


2k  2λλ−
φ0 − 2 − 2 (n(n + 1) − 2) unm
r0
r0
2r0




1 −
wr0
(n + 3)
1− − −P+ −2(n + 1)(2n + 3)
η
2 nm
w +
wr0 (n − 2)
P + = 0,
− ψnm
+ +
r0
η 2n(2n − 1) nm

APPENDIX C: ELEMENTS OF THE MATRIX C

C11 =The lateral force balance along the θ -direction yields


k
η+
wr0

+
δφnm + δφnm + 2 2n + 4 + + ψnm
r0
η
r0


2k  2λλ+
φ0 − 2 − 2 (n(n + 1) − 2) unm
r0
r0
2r0




(n − 2)
w −
1 +
wr0
+ −
1+ + −
Pnm − ψnm
2n(2n − 1)
η
2
r0
−

PHYSICAL REVIEW E 94, 062414 (2016)

2n(n + 1) [w [(n + 1) + (n − 1)E] − (2n + 1)2 E],D(C3)n(n + 1)(n2 − 1)(2n2 + 4n + 3)E 2
C22 =
D


3
+ 2n(2n3 + 4n2 + n − 1)
+E
2

2
(C4)
+ n(n + 2)(2n + 1) ,n(n + 1)(n2 − 1)(2n2 + 4n + 3)E 2
C23 = C32 =D92+ (2n + 1)E − n(n + 2)(2n + 1) , (C5)2n(n + 1) 2[(n − 1)(2n2 + 4n + 3)E 2
C33 =D+ 2w  (n − 1)(2n2 + 5n + 5)E + [5 − 8n2 (n + 1)2 ]E=−+ (n + 2)[2w (n(2n − 1) + 2) + n(2n2 + 1)]].
(C6)
The common denominator in all the components is(B4)

for the outer and inner monolayers, respectively.
Equations (B1)–(B4) constitute a set of four simultaneous±±
linear equations, which can be used to eliminate ψnm
and Pnm.
Then we obtain the dynamical equations for the variables u,
φ  , and φ .

[1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, andP. Walter, Molecular Biology of the Cell (Garland Science,
New York, 2008).
[2] R. Phillips, J. Kondev, J. Theriot, and H. Garcia, Physical
Biology of the Cell (Garland Science, New York, 2012).
[3] R. J. Ellis and A. P. Minton, Nature (London) 425, 27(2003).
[4] S. B. Zimmerman and A. P. Minton, Annu. Rev. Biophys.
Biomol. Struct. 22, 27 (1993).
[5] D. S. Banks and C. Fradin, Biophys. J. 89, 2960 (2005).
[6] D. Miyoshi and N. Sugimoto, Biochimie 90, 1040 (2008).
[7] I. M. Kuznetsova, K. K. Turoverov, and V. N. Uversky, Int. J.
Mol. Sci. 15, 23090 (2014).

D = 2{2E 2 (n − 1)(n + 1)(2n2 + 4n + 3)[w + (2n + 1)]+ E[w (8n4 + 16n3 + 4n2 − 4n + 3)
+ 2n(n + 2)(2n + 1)(2n2 + 1)]+ 2w  n(n + 2)(2n2 + 1)}.

(C7)[8] R. J. Ellis, Trends Biochem. Sci. 26, 597 (2001).
[9] K. Kurihara, M. Tamura, K. Shohda, T. Toyota, K. Suzuki, andT. Sugawara, Nat. Chem. 3, 775 (2011).
[10] Y. Kuruma and T. Ueda, Nat. Protocols 10, 1328 (2015).
[11] K. Fujiwara and M. Yanagisawa, ACS Synth. Biol. 3, 870(2014).
[12] F. Brochard and J. F. Lennon, J. Phys. France 36, 1035
(1975).
[13] D. Boss, A. Hoffmann, B. Rappaz, C. Depeursinge, P. J.
Magistretti, D. V. de Ville, and P. Marquet, PLoS ONE 7, e40667(2012).
[14] T. Betz, M. Lenz, J.-F. Joanny, and C. Sykes, Proc. Natl Acad.
Sci. U.S.A. 106, 15320 (2009).

062414-13

SACHIN KRISHNAN, OKAMOTO, AND KOMURAPHYSICAL REVIEW E 94, 062414 (2016)[15] Y. Park, C. A. Best, T. Auth, N. S. Gov, S. A. Safran, G. Popescu,
S. Suresh, and M. S. Feld, Proc. Natl Acad. Sci. U.S.A. 107, 1289
(2010).
[16] H. Turlier, D. A. Fedosov, B. Audoly, T. Auth, N. S. Gov,
C. Sykes, J.-F. Joanny, G. Gompper, and T. Betz, Nat. Phys. 12,513 (2016).
[17] S. T. Milner and S. A. Safran, Phys. Rev. A 36, 4371 (1987).
[18] M. B. Schneider, J. T. Jenkins, and W. W. Webb, J. Phys. (Paris)
45, 1457 (1984).
[19] L. Miao, M. A. Lomholt, and J. Kleis, Eur. Phys. J. E 9, 143(2002).
[20] M. Mell, L. H. Moleiro, Y. Hertle, I. López-Montero, F. J. Cao,
P. Fouquet, T. Hellweg, and F. Monroy, Chem. Phys. Lipids 185,
61 (2015).
[21] M. Arroyo and A. DeSimone, Phys. Rev. E 79, 031915 (2009).
[22] S. Komura and K. Seki, Physica A 192, 27 (1993).
[23] K. Seki and S. Komura, Physica A 219, 253 (1995).
[24] S. Komura, in Vesicles, edited by M. Rosoff (Marcel Dekker,
New York, 1996), p. 198.
[25] E. Evans, A. Yeung, R. Waugh, and J. Song, in The Structure
and Conformation of Amphiphilic Membranes, Proceedings ofthe International Workshop on Amphiphilic Membranes, editedby R. Lipowsky, D. Richter, and K. Kremer (Springer, Berlin,1992), p. 148.
[26] U. Seifert and S. A. Langer, Europhys. Lett. 23, 71 (1993).
[27] R. Okamoto, Y. Kanemori, S. Komura, and J.-B. Fournier, Eur.
Phys. J. E 39, 52 (2016).
[28] J.-B. Fournier, Int. J. Nonlinear Mech. 75, 67 (2015).
[29] M. Doi, J. Phys.: Condens. Matter 23, 284118 (2011).
[30] M. Doi, Soft Matter Physics (Oxford University Press, Oxford,2013), Chap. 7.
[31] S. R. Keller and R. Skalak, J. Fluid Mech. 120, 27 (1982).
[32] V. Kantsler and V. Steinberg, Phys. Rev. Lett. 95, 258101(2005).
[33] V. Kantsler and V. Steinberg, Phys. Rev. Lett. 96, 036001(2006).
[34] J. F. Faucon, M. D. Mitov, P. Méléard, I. Bivas, and P. Bothorel,
J. Phys. (Paris) 50, 2389 (1989).
[35] W. Rawicz, K. C. Olbrich, T. McIntosh, D. Needham, and E.
Evans, Biophys. J. 79, 328 (2000).

[36] E. Sackmann, in Structure and Dynamics of Membranes, Vol. 1A
From Cells to Vesicles, edited by R. Lipowsky and E. Sackmann(Elsevier, Amsterdam, 1995), p. 213.
[37] R. Aris, Vectors, Tensors, and the Basic Equations of Fluid
Mechanics (Dover, New York, 1962).
[38] O.-Y. Zhong-can and W. Helfrich, Phys. Rev. A 39, 5280
(1989).
[39] W. Helfrich, Z. Naturforsch. C 28, 693 (1973).
[40] B. Różycki and R. Lipowsky, J. Chem. Phys. 142, 054101(2015).
[41] S. Leibler, J. Phys. 47, 507 (1986).
[42] S. Leibler and D. Andelman, J. Phys. 48, 2013 (1987).
[43] L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Pergamon,
Oxford, 1987).
[44] S. A. Safran, J. Chem. Phys. 78, 2073 (1983).
[45] P. Méléard, M. D. Mitov, J. F. Faucon, and P. Bothorel, Europhys.
Lett. 11, 355 (1990).
[46] G. Niggemann, M. Kummrow, and W. Helfrich, J. Phys. II
France 5, 413 (1995).
[47] F. Youhei, Physica A 203, 214 (1994).
[48] H. Lamb, Hydrodynamics (Cambridge University Press, London, 1975).
[49] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Kluwer Academic, The Hague, 1973).
[50] T. Pott and P. Méléard, Europhys. Lett. 59, 87 (2002).
[51] R. Rodrı́guez-Garcı́a, L. R. Arriaga, M. Mell, L. H. Moleiro, I.
López-Montero, and F. Monroy, Phys. Rev. Lett. 102, 128101(2009).
[52] A. Viallat, J. Dalous, and M. Abkarian, Biophys. J. 86, 2179
(2004).
[53] A. J. Levine and F. C. MacKintosh, Phys. Rev. E. 66, 061606(2002).
[54] S. B. Rochal, V. L. Lorman, and G. Mennessier, Phys. Rev. E.
71, 021905 (2005).
[55] S. Komura, S. Ramachandran, and K. Seki, Europhys. Lett. 97,
68007 (2012).
[56] S. Komura, S. Ramachandran, and K. Seki, Materials 5, 1923
(2012).
[57] S. Komura, K. Yasuda, and R. Okamoto, J. Phys.: Condens.
Matter 27, 432001 (2015).

062414-14