PHYSICAL REVIEW E 96, 012416 (2017)

Dynamics of a bilayer membrane coupled to a two-dimensional cytoskeleton: Scale transfers ofmembrane deformations
Ryuichi Okamoto,1 Shigeyuki Komura,1 and Jean-Baptiste Fournier2
1

2

Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397, Japan
Université Paris Diderot, Sorbonne Paris Cité, Laboratoire Matière et Systèmes Complexes (MSC), UMR 7057 CNRS, F-75205 Paris, France(Received 23 February 2017; published 31 July 2017)
We theoretically investigate the dynamics of a floating lipid bilayer membrane coupled with a two-dimensional
cytoskeleton network, taking into account explicitly the intermonolayer friction, the discrete lattice structure of thecytoskeleton, and its prestress. The lattice structure breaks lateral continuous translational symmetry and couples
Fourier modes with different wave vectors. It is shown that within a short time interval a long-wavelength deformation excites a collection of modes with wavelengths shorter than the lattice spacing. These modes relax slowly witha common renormalized rate originating from the long-wavelength mode. As a result, and because of the prestress,
the slowest relaxation is governed by the intermonolayer friction. Conversely, and most interestingly, forces
applied at the scale of the cytoskeleton for a sufficiently long time can cooperatively excite large-scale modes.
DOI: 10.1103/PhysRevE.96.012416

I. INTRODUCTION

In biological materials, much attention has been paid tothe dynamics from the viewpoint of nonequilibrium physics,
because of the high complexity of composition, hydrodynamic
interactions, active components, etc. [1–3]. In particular, shaperelaxation and fluctuations of lipid bilayer membranes haveintensively been studied both in, or near, equilibrium [4–6]
and far from equilibrium [7–12]. The dynamics of a bilayer
membrane is determined by many factors, such as the viscosityof the surrounding fluid [4], the membrane bending rigidity,
the intermonolayer friction caused by relative lateral motionsof two monolayers [6,13], and possibly active inclusions [7].
Red blood cell (RBC) membranes have further complexity
because of the cytoskeleton which is attached to the lipid
bilayer. It has been argued that the cytoskeleton plays crucial
roles in both the statics and the dynamics, e.g., drastic effective
tension increase in equilibrium [14–17], tension decrease in
the presence of adenosine triphosphate (ATP) [8,12], andenhanced nonequilibrium fluctuations on the scale of thecytoskeleton mesh size [9].
In RBCs, the cytoskeleton consists of spectrin filaments
forming a prestressed [12] two-dimensional (2D) triangularlattice with a protein at each vertex embedded in the bilayer
membrane. The lattice spacing a ≈ 100 nm is quite large, andtherefore what matters in the membrane collective dynamics
is not only the modes whose wavelengths are much larger
than a, but also those having wavelengths smaller than a. To
understand simultaneously the dynamics on such a wide rangeof spatial scales, we need to take into account explicitly thediscrete nature of the lattice structure. The latter breaks lateralcontinuous translational symmetry, giving rise to a coupling
between modes on different length scales and thus to a rich
dynamical behavior.
This paper is organized as follows. In Sec. II, our model
free energy is constructed on the basis of the previous theories
for a bilayer membrane (without a cytoskeleton) [6] and for aprestressed 2D cytoskeleton coupled to a membrane [16]. In
Sec. III, following Ref. [6], the hydrodynamic equations are
introduced, where we take into account the hydrodynamic
flows of the surrounding fluid and of the monolayers and2470-0045/2017/96(1)/012416(10)intermonolayer friction between the monolayers. Then we
obtain the coupled equations for the membrane variables byintegrating out the flow velocity fields. In Sec. IV, we discuss
how the cytoskeleton alters the dynamics of the membrane
when initially a large-scale deformation is imposed and whenforce(s) are applied to small-scale mode(s) for a long time( 10 ms). Section V is devoted to a discussion and summary.
We also present our detailed calculations in the appendixes.
II. FREE ENERGY

We consider out-of-plane deformations of a RBC membrane patch described by its height h(x,y) above a reference
plane z = 0. Our model free energy is given by F = F0 + Fc
as follows. We take into account the areal compression which isnecessarily coupled with h due to the finite thickness d ≈ 1 nmof the monolayers [6]. Then the bilayer membrane free energy
F0 is given by


κ
σ
2
F0 = d x (∇ 2 h)2 + (∇h)2
2
2

k 
2
2
+
[ρ + d(∇ h)] ,
(1)
2 =±
where κ is the bare bending rigidity, σ the bare tension, k theareal compression modulus, and ∇ = (∂x ,∂y ). In the above, ρ +(ρ − ) denotes the dimensionless projected excess lipid density
in the upper (lower) monolayer [6].
Note that in our paper the surface tension is not considered
as a constant. The quantity σ in Eq. (1) is only the backgroundtension for a flat membrane with homogeneous, reference
lipid density (ρ ± = 0). The actual tension fluctuates about
the zeroth-order tension σ according to a model that isclosely related to the area-difference-elasticity model [18].
Indeed, the term proportional to k in Eq. (1) that involves
the variables ρ + and ρ − is the excess energy associatedwith a local compression or dilation of the lipids in eachmonolayer. The actual tension (without the cytoskeleton) is
σ + k(ρ + + d∇ 2 h) and σ + k(ρ − − d∇ 2 h) in the upper andlower monolayers, respectively.

012416-1©2017 American Physical Society

OKAMOTO, KOMURA, AND FOURNIERPHYSICAL REVIEW E 96, 012416 (2017)

can have different fluid velocities, v + and v − , respectively[Fig. 1(a)]. The full dynamic equations consist of (i) lateralforce balance for each monolayer, (ii) force balance normal tothe bilayer, and (iii) the continuity equation for lipids in eachmonolayer.
We use the Stokes equation for the solvent velocity field V
and the pressure field p, with the shear viscosity η,(a) Red blood cell membraneˆ = 0,
η∇ˆ 2 V − ∇p(b) Cytoskeleton lattice∇ˆ · V = 0,

where ∇ˆ = (∂x ,∂y ,∂z ) is the 3D nabla operator. No-slip
boundary condition is employed at the membrane surface,vi± = Vi (i = x,y) and Vz = ∂h/∂t at z → 0± . We also
impose V → 0 and p → p0 as z → ±∞. The 2D viscousstress tensors in the monolayers are given by(c) Reciprocal lattice

FBZ

τij± = μ(∂i vj± + ∂j vi± ) + (ζ − μ)δij ∇ · v ± ,FIG. 1. (a) Illustration of the model red blood cell (RBC)
membrane. (b) Cytoskeletal network. (c) Reciprocal lattice. The
yellow region represents the first Brillouin zone (FBZ).

The other contribution, Fc , arises from the membranecytoskeleton coupling. We assume, to simplify, that thecytoskeleton network is a uniform triangular lattice without
defects. An anchoring protein at each lattice site is embeddedin the membrane and interacts with its nearest neighborthrough an effective spring of relaxed length a0 and stiffness
ks [Fig. 1(a)]. In the ground state (ρ = h = 0), the network
forms a regular triangular lattice, and the lattice points {R } areexpressed in terms of the primitive lattice vectors eα (α = 1,2)
as R = Rα eα , with Rα ∈ Z positive or negative integers [seeFig. 1(b)]. The lattice spacing is a = |eα |. If the out-of-plane
deformation of the membrane is sufficiently small, Fc is givenby [16]ν 
Fc =[h(R ) − h(R + n)]2 ,
(2)
4  n
where ν = ks (1 − a0 /a) is the effective stiffness of theharmonic potentials
associated with the out-of-plane defor
mations, and n denotes the sum over the nearest-neighborsites. Note that ν is nonzero only if the lattice is prestressed(a = a0 ). In Ref. [12], as a result of fitting their experimental
data, it has been shown that the cytoskeleton in healthy
RBCs is naturally stretched (by about 4%), while the baremembrane tension can be negative. Let us introduce thein-plane Fourier transform as Fq [· · · ] ≡ d 2 x (· · · )e−iq·x andthe reciprocal vectors eα satisfying eα · eβ = δβα , the Kroneckerδ [see Fig. 1(c)]. Because Fc breaks the lateral continuoustranslational symmetry, modes with different wave vectors arecoupled to one another. As shown in Appendix A, the modes
coupled to a given q belong to the subset
Qq = {q + 2π mα eα |mα ∈ Z}.

(4)

(3)

III. DYNAMIC EQUATIONS

Following Seifert and Langer [6], we regard each monolayer as a compressible 2D fluid having the shear viscosity
μ and the bulk viscosity ζ . The upper and lower monolayers(5)where the superscript “+” (“−”) denotes the upper (lower)monolayer. Then the lateral force balance equation in eachmonolayer readsδF+ ∂j τij± ± Tiz± ∓ b(vi+ − vi− ) = 0,
(6)
−∂i
δρ ±
where Tij+ (Tij− ) is the stress tensor Tij = −pδij + η(∂i Vj +∂j Vi ) in the solvent fluid evaluated at z → 0+ (z → 0− ). The
last term is due to the intermonolayer friction, with the friction
coefficient b [6]. In the normal direction, the force exerted by
the surrounding fluid is balanced with the restoring force ofthe membrane,
δF
.
(7)
δh
At linear order in v ± and ρ ± , which are both considered to besmall, the continuity equation in each monolayer is given by
Tzz+ − Tzz− =

∂ρ ±
(8)
−∇ · v ± .
∂t
The velocities, V and v ± , can be eliminated from the
dynamic equations by integrating the Stokes equations alongz for each mode q (see Appendix B ). This yields coupledlinear equations for ĥ ≡ h/d and ρ ≡ (ρ + − ρ − )/2:
δF
∂ ĥ(q,t)+ uh (q,t),
= −Fq
∂t
δ ĥ
δF
2c(q) ∂ρ(q,t)= −Fq+ uρ (q,t),
q2
∂t
δρ

4ηd 2 q(9)(10)

where c(q) = 2b + 2ηq + (μ + ζ )q 2 , with q = |q|. We
have added uh and uρ , representing external forces appliedmechanically (e.g., by active molecules) which act on thevariables ĥ and ρ, respectively. Since Fq [δF /δh] includesĥ(q ) for ∀ q ∈ Qq , these equations actually consist of setsof coupled equations for the variables {ĥ(q ),ρ(q )} in eachset Qq (see Appendix A ).
Without the cytoskeleton, the modes for different wave
vectors are not coupled in Eqs. (9) and (10). Then, ĥ(q)and ρ(q) exhibit two relaxation rates, γ+(0) (q) > γ−(0) (q),associated with some linear combinations of ĥ(q) and ρ(q).
Seifert and Langer discussed these relaxation modes forvanishing tension [6]. They found a crossover wave number

012416-2

DYNAMICS OF A BILAYER MEMBRANE COUPLED TO A . . .

PHYSICAL REVIEW E 96, 012416 (2017)

TABLE I. List of the parameter values used in numerical calculations.
σ
(N/m)

κ
(J)

k
(N/m)

d
(m)

ν
(N/m)

a
(m)

η
(J s/m3 )b(J s/m4 )μ+ζ
(J s/m2 )10−11

2 × 10−20

7 × 10−2

10−9

10−6

10−7

10−3

2 × 108

2 × 10−9

qc = 2ηk/(bκ̃) ≈ 4.4 × 106 m−1 , at which the relaxation behavior of the membrane changes qualitatively. Here we set κ =2 × 10−20 J as in [15] (the value of κ measured in experiments
lies in quite a wide range, 1 to 30 × 10−20 J [4,8,12,17],
but the following results remain almost unchanged even withthese different values). For large scales satisfying q
qc , therates correspond to ρ relaxing quickly followed by h relaxing
slowly, with ρ being slaved [19]. For small scales, q  qc ,conversely, they correspond to h relaxing quickly followed byρ relaxing slowly, with h being slaved. Hence, the dynamics on
the small scales is dominated by the inter-monolayer friction,
whereas that on the large scales is dominated by the solvent
viscosity. In the presence of tension, their results hold for σσc ≡ (2ηk)2 /(κ̃b2 ), except at very large scales (see Refs. [19–21] and Appendix C). However, for σ  σc the dynamics isdominated at all scales by the intermonolayer friction, with
γ+(0) (σ q + κq 3 )/(4η) > γ−(0) ≈ kq 2 /(2b) [19].

√in the direction q/q = ( 3/2,1/2), and we integrate numerically the dynamical equations up to the cutoff 20 × π/a, of
the order of the inverse membrane thickness. Experiments
indicate σeff ≈ ν ≈ 10−7 –10−5 N/m, with σ very small oreven negative [5,10,12,17]. Accordingly, besides the values
already given, we set σ = 10−11 N/m and ν = 10−6 N/m,yielding κeff = 1.9 × 10−20 J and σeff = 1.73 × 10−6 N/m.
We study the coupled evolution of ĥ and ρ for q and forsmall-scale modes q ∈ Qq . In Fig. 2, we present as an example

1
10 -2
10 -4
10 -6

IV. RESULTS

10 -8

A. Relaxation of a large-scale deformation

10 -10The cytoskeleton shifts the mode relaxation rates by anamount that depends on the prestress ∼ν, and at the same times
it couples all the modes belonging to a common set Qq . Let us
first discuss how the rates of the large-scale modes, with qqc
2π/a, are shifted by the cytoskeleton. For such modes,
the dependence on the direction of q is negligible. In thefollowing, analytical expressions will be given systematicallyat first order in a perturbative expansion in a power series of ν(see Appendix C ). The parameter values used in the following
numerical calculations are summarized in Table I.
From the dynamic equations (9) and (10), we find that therates of the large-scale modes, γ+ > γ− , are shifted according
to
√ 2

3d b
σeff q + κeff q 3
k
, γ−
+
γ+
ν q 2 , (11)
4η
2b
4η2
√
√
where σeff = σ + 3ν and κeff = κ − 3νa 2 /16 are thetension and the bending rigidity renormalized by thecytoskeleton [12,16]. Note that the fast and slow rates
have been exchanged with respect to their bare value
γ+(0) kq 2 /(2b) > γ−(0) (σ q + κq 3 )/(4η) because weanticipate σeff  σc > σ for large-enough ν. The other rates
γ± (q ) of the subset q ∈ Qq associated with q are also
shifted from their bare values γ±(0) (q ); they correspond towavelengths comparable to or smaller than the cytoskeletonmesh size and are much faster than those in Eq. (11). Ourdetailed calculations show that the shifts of these rates are
small (see Appendix C), but this does not mean that thesmall-scale modes are not affected by the cytoskeleton
The relaxation of a large-scale mode q excites all the smallscale modes in Qq (Fig. 2). To investigate this effect, we set theinitial condition ĥ(x) = eiq·x and ρ(x) = 0, with q = 106 m−1(a)

10 -12 -8
10

10 -7

10 -6

10 -5

1
10

10 -4

(b)
-2

10 -4
10 -6
10 -8
10 -10
10 -12
0

0.01

0.02

0.03

(c)

FIG. 2. Relaxation of a (normalized) large-scale deformationĥ(x) = eiq·x . (a) Short-time behavior of the large-scale modes ĥ(q)and ρ(q), and that of the small-scale modes ĥ(q ) and ρ(q ) withq = q + 2π e2 . (b) Long-time behavior of the same variables.
(c) Schematic pictures of the process. Although the amplitude ofthe excited small-scale mode is small, a large effective tension
σeff changes qualitatively the large-scale dynamics such that it isdominated by the intermonolayer friction.

012416-3

OKAMOTO, KOMURA, AND FOURNIERPHYSICAL REVIEW E 96, 012416 (2017)

only q = q + 2π e2 2π e2 , as we find the other small-scalemodes in Qq also exhibit a similar behavior. In the shorttime interval 0 < t1/γ± (q ) ≈ 1/γ±(0) (q ), the small-scalemodes ĥ(q ) and ρ(q ), that are initially zero, are excited, whileĥ(q) almost remains unchanged [Figs. 2(a) and 2(c)]. All the
excited small-scale modes rapidly approach their respective
quasiequilibrium states ĥqe and ρqe , which minimize the freeenergy for a fixed value of ĥ(q), given byρqe (q ; ĥ(q))
d 2q 22ν ĥ(q)Kq−√,3κa 2 q 4

10 8
(a)
10

6

10 4
10 2

(12)

1


where Kq = n (1 − eiq·n ). Figure 2(b) illustrates the longt  γ+−1 , ĥ(q),time evolution of the system. For 1/γ±(0) (q )ĥ(q ), and ρ(q ) decay with the common rate γ+ . Then, aroundt ≈ γ+−1 , ĥ(q) follows the dynamical quasiequilibrium valueĥρqe (q; ρ(q,t)) that minimizes the free energy at fixed ρ(q,t).
Finally, for t  γ+−1 , all the modes decay with the common
rate γ− , with ĥ(q ), ρ(q ), and ĥ(q) following their respectiveρdynamical quasiequilibrium values ĥ(q ) ĥqe (q ; ĥqe (q)),ρρρ(q ) ρqe (q ; ĥqe (q)), and ĥ(q) ĥqe (q; ρ(q,t)) [Figs. 2(b)and 2(c)]. Suppose in Fig. 2 the initial amplitude of h(q)is comparable with the mode wavelength, 2π/q ≈ 6.3 μm.
Then the amplitude of the excited small-scale mode h(q )is about 2.4 × 10−5 μm, which is much smaller than themode wavelength 2π/q ≈ 8.6 × 10−2 μm and may not be
observable in experiments. This is because an energy cost to
make a deformation with amplitude q −1 at the small-scale q
is larger than to make a deformation with amplitude q −1 at thelarge-scale q. Nevertheless, we notice that the cytoskeleton
alters qualitatively the large-scale dynamics; because of
the cytoskeleton that yields σeff ∼ σc ≈ 3 × 10−6 N/m, theslowest relaxation process is dominated by the large-scale
compression mode ρ(q) limited by the intermonolayer friction.

10 -2ĥqe (q ; ĥ(q)) =B. Large-scale deformation induced by small-scale deformation
via the cytoskeleton

We have seen that large-scale deformations excite themodes whose scales are comparable to or smaller than thecytoskeleton mesh. Now a question arises. Can small-scale deformations excite large-scale ones? If yes, is the amplitudeof the excited modes large enough to be observable? The
answer is no, if there are no applied forces. This is because thesmall-scale modes are much faster than the large-scale ones,
so that the small-scale modes rapidly relax before large-scalemodes are excited. However, if we keep applying forces only
to the small-scale modes for a time longer than the relaxation
time of the large-scale modes, the latter will be excited via
the cytoskeleton. Furthermore, if forces are applied to manysmall-scale modes, the amplitude of the excited large-scalemode can be noticeably large. In fresh RBCs, active molecules
could be the source of these forces, as their characteristic timeis of the order of 1 s [8,12] which is much larger than the
typical relaxation time of the modes for q ∼ 106 m−1 . It was
further proposed that the active force is particularly enhancedon the scales of the cytoskeleton mesh [9].
To study this, we choose again q = 106 m−1 , oriented as
before, and we apply a constant force uh (q 1 ) = ū only to

10 -7

10 -5

10 -3

10 -1

10

(b)
Force

FIG. 3. Time evolution of the membrane shape under a constantforce ū applied to the small-scale mode ĥ(q 1 ), with q 1 = q + 2π e1 .
All√ values of ĥ are normalized by ū. The set Qq is chosen as q/q =( 3/2,1/2), with q = 106 m−1 . (a) ĥ as a function of t for q, q 1 , andq 2 = q + 2π e2 . (b) Schematic pictures of the process.

ĥ(q 1 ), with q 1 = q + 2π e1 ∈ Qq (Fig. 3). The corresponding
wavelengths of q and q 1 are then 2π/q ≈ 6.3 μm and
2π/|q 1 | ≈ 0.87a, respectively. We investigate the response ofĥ at the large scale q but also at another small scale, q 2 =q + 2π e2 ∈ Qq . With the initial condition ĥ = ρ = 0, all themodes will be proportional to ū. For t  1/γ±(0) (q1 ), the smallscale deformation ĥ(q 1 ) is excited and reaches the stationary
value ĥst (q 1 ) ū/(κd 2 q14 ) minimizing F − (2π )−2 ūĥ(q 1 ).
Then, for 1/γ±(0) (q1 )  t  γ±−1 , the large-scale mode ĥ(q)gets excited by ĥ(q 1 ) via the cytoskeleton deformation, andthen for t  γ±−1 , ĥ(q) reaches the stationary valueĥst (q)−σeff + κeff q 2 ū,σeff + κeff q 2 κd 2 q14(13)

where σeff = σeff − σ and κeff = κeff − κ. With our choiceof parameters, (σeff + κeff q 2 )/(σeff + κeff q 2 ) 1, and thusĥst (q) −ĥst (q 1 ), consistent with Fig. 3(a). We find that theother small-scale modes, such as ĥ(q 2 ), are also excited, but
not significantly [Fig. 3(a)].
When a force distribution is applied to multiple small-scalemodes, the magnitude of the excited large-scale mode ĥ(q)can become much larger than in the case examined above. To
show this, let us consider at each lattice site of the cytoskeleton
active forces inducing some local curvature. Such forces
can formallybe derived by adding a “fictitious” potentialU = −  w d 2 xδ(x − R )∇ 2 h to the freeenergy. This
yields u(q ) = q 2 wq , ∀ q ∈ Qq , with wq =  w e−iq·R .
Note that some force is applied also to the large-scale mode
q. By linearity, the effect of the small-scale modes q on
the large-scale mode q, denoted by δ ĥ(q), will be enhanced

012416-4

DYNAMICS OF A BILAYER MEMBRANE COUPLED TO A . . .


by a factor r = m1 ,m2 q12 /|2π mα eα |2 with respect to thecase shown in Fig. 3, where only the mode q 1 was excited.
The sum is taken up to the high wave vector cutoff while
excluding m1 = m2 = 0. With the parameters given above, wefind r ≈ 20. Assuming that the microscopic forces can produce
a deformation of amplitude comparable to the mesh size, i.e.,
h1 50 nm [8,9], we expect the large-scale response δh to beabout rh1 , which is a sizable deformation of the order of 1 μm.
Note that this scale transfer of membrane deformation requires
applying the small-scale forces for at least about 10 ms (Fig. 3).
V. DISCUSSION AND SUMMARY

In this paper, since we assume only vertical motion of thequasiplanar RBC patch, and equilibrium dynamics, we neglect
the dissipation due to tangential motion of the cytoskeleton as
well as the cytoskeleton activity that was studied in Ref. [12].
Nevertheless, in Appendix D, we have considered the friction
between the tangential monolayer flow and the anchored proteins, yielding an extra contribution to the lateral force balance
equation. However, in our detailed calculation, it is shown to be
negligible. As for the viscous drag of the spectrin filaments dueto the surrounding fluid, it was also shown to be negligible [12].
For simplicity, we have assumed a quasiplanar membrane,i.e., small deformations about a flat reference shape. However,real RBCs are intrinsically curved objects that fluctuate about
a curved reference shape [12,22]. For such cases, not only
ρ = (ρ + − ρ − )/2 but also ρ̄ = (ρ + + ρ − )/2 is coupled tothe membrane deformation h [21]. Furthermore, for a curved
membrane, the tangential deformation of the cytoskeleton is
also coupled to h [12]. For the relaxation of ρ̄, we can showthat the intermonolayer friction is not a dissipation source,
while the friction between the anchored proteins and themonolayers is one of the major dissipation sources for large
λ/a 2 . However, in thescales satisfying 2ηq + (μ + ζ )q 2
relaxation of ρ̄ for a bilayer without the cytoskeleton, theinertia effect of the surrounding fluid cannot be neglected [6],so that a more careful study is necessary in the future.
In summary, we have studied the dynamics of RBC membranes modeled as bilayers coupled to a prestressed discrete
elastic network, and subject to viscous dissipation in thesolvent, in each monolayer and between the monolayers. Given
the mesh size of the cytoskeleton (≈100 nm), it is importantfrom the biological point of view to address the dynamics atscales both larger and smaller than the cytoskeleton. Because
the latter breaks lateral translational symmetry, each mode
is coupled to all the modes that are congruent modulo awave vector of the cytoskeleton’s reciprocal lattice. We havecharacterized how the small modes renormalize the relaxation
rates of the large modes. We have found that, because of thelarge renormalized tension σeff , the shape relaxation dynamics
on the large scales is dominated by the intermonolayer frictionthat has regularly been neglected in the previous theories
on RBC dynamics [12,14,22]. It has been also shown that
applying forces on the small-scale modes for a sufficiently
long time can excite large-scale deformations.
It is informative to measure the correlations between
different Fourier modes, h(q,t)h(q ,t ), in order to know the
precise dynamical processes where modes in different scales
are coupled due to the cytoskeleton and also to understand the
behavior of active forces in fresh RBCs.

PHYSICAL REVIEW E 96, 012416 (2017)
ACKNOWLEDGMENTS

R.O. and S.K. acknowledge support from a Grant-inAid for Scientific Research on Innovative Areas Fluctuationand Structure (Grant No. 25103010) from the Ministry of
Education, Culture, Sports, Science, and Technology of Japanand from a Grant-in-Aid for Scientific Research (C) (Grant
No. 15K05250) from the Japan Society for the Promotion of
Science (JSPS).
APPENDIX A: COUPLING BETWEEN DIFFERENTFOURIER MODES

Notice that the Fourier modes between different wave
vectors are coupled with each other because the cytoskeleton network breaks the continuous translational symmetry.
More precisely, while Fq [δF /δρ] = 2k[ρ(q) − dq 2 h(q)] in
Eq. (B11) does not couple different Fourier modes, Fq [δF /δh]in Eq. (B9) couples the Fourier modes in the set Qq ={q + 2π mα eα |mα ∈ Z}. To see this, we rewrite Eq. (2) as1ν 
δ(x − R )[h(x) − h(x + n)]2 . (A1)
Fc =
d 2x
2
2  n
Then we can calculate the Fourier-transformed functional
derivative as
δF= (σ q 2 + κ̃q 4 )h(q) − 2kdq 2 ρ(q)
Fq
δh

+ν
e−iq·R [h(R ) − h(R + n)]. (A2)n

Now we use the identity(2π )2 
e−iq·R = √δ(q ),
g q ∈Q(A3)

q

where g = |eα · eβ | is the determinant of the metric tensor in
the primitive vector frame {eα } and is given by g = 3a 4 /4.
Using this identity, we can rewrite Eq. (A2) as
Fq

δF
δh

= (σ q 2 + κ̃q 4 )h(q) − 2kdq 2 ρ(q)ν 
h(q )Kq ,+√g q ∈Q(A4)q
where Kq = n (1 − eiq·n ). We can clearly see from Eq. (A4)that in Eq. (9) the Fourier modes in the common set Qq are
coupled to each other. Note that without the cytoskeleton ν = 0(and without the force terms uh and uρ ), Eqs. (9) and (10)reduce to the equations studied by Seifert and Langer [6].
APPENDIX B: ELIMINATION OF THE VELOCITY FIELDSIn accordance with the in-plane Fourier transform with awave vector q, we introduce q ⊥ = ez × q, where ez is theunit vector pointing towards the z direction. Then unit vectors
q̂ and q̂ ⊥ are defined as q̂ = q/q and q̂ ⊥ = q ⊥ /q = ez × q̂,
respectively. The longitudinal and transverse components of
the Fourier transform of V and v ± are defined by V (q,z) =q̂ · V (q,z), V⊥ (q,z) = q̂ ⊥ · V (q,z), v± (q) = q̂ · v ± (q), and±v⊥(q) = q̂ ⊥ · v ± (q). The Fourier transform of Eq. (4) is012416-5

OKAMOTO, KOMURA, AND FOURNIERPHYSICAL REVIEW E 96, 012416 (2017)

written asη ∂z2 − q 2 Vz − ∂z p = 0,(B1)

η ∂z2 − q 2 V − iqp = 0,(B2)

η ∂z2 − q 2 V⊥ = 0,

(B3)

iqV + ∂z Vz = 0.

(B4)

These equations are solved to obtain
± ∓qz

p=p e

∂
|(t)q = −(SL + νh )|(t)q ,(C2)∂twhere the linear operators SL and h are defined by kq 2
kd 2 q 4|ρ; q ρ; q | −|ρ; q h; q |SL (q) =
c(q )
c(q )
q ∈Q
q

,

(B5)

± ∓qz
e ,
V⊥ = v⊥(B6)+= [v± − (i ḣ ± v± )qz]e∓qz ,(B8)

where the upper and the lower signs indicate the solutions forz > 0 and for z < 0, respectively.
Substituting Eqs. (B5) and (B7) into the Fourier transformof Eq. (7), we obtain
4ηq

δF
∂h(q)
= −Fq
,
∂t
δh

(B9)

which yields Eq. (9) (without the force term uh ). Next, we
use Eq. (B8) to eliminate V and Vz from the longitudinalcomponent of Eq. (6), and obtain
0 = [2ηq + (μ + ζ )q 2 ]v± (q) + iqFq

δF
δρ ±

± b[v+ (q) − v− (q)].

1   Kqh (q) = √|h; q h; q |.
g q ∈Q q ∈Q 4ηqq(SL + νh )|e = γ |e,

δF
∂ρ(q)
2
δF
c(q)
+ iqFq
− − ,
iq
∂t
δρ +
δρ

(B11)

where c(q) = 2b + 2ηq + (μ + ζ )q 2 . Using the identityδ(· · · )/δρ = δ(· · · )/δρ + − δ(· · · )/δρ − , we obtain Eq. (10)(without the force term uρ ).
APPENDIX C: OPERATOR REPRESENTATION ANDPERTURBATION EXPANSION

We seek the relaxation rates and their associated eigenmodes of Eqs. (9) and (10) as a power series of ν, sincefor ν = 0 they can be obtained analytically. It is convenient
to introduce for each Qq a Hilbert space Sq spanned by anorthonormal set {|h; q , |ρ; q  | q ∈ Qq }, and its dual space
Sq∗ spanned by an orthonormal set {h; q |, ρ; q | | q ∈ Qq }.
We define a state vector in Sq ,|(t)q ≡[ĥ(q ,t)|h; q  + ρ(q ,t)|ρ; q ],
(C1)

(C5)

where γ and |e are the eigenvalue and the eigenvector ofSL + νh , respectively. We expand γ and |e in powers ofν as γ = γ (0) + νγ (1) + · · · and |e = |e(0)  + ν|e(1)  + · · · .
The zeroth- and the first-order equations read
SL |e(0)  = γ (0) |e(0) ,(h − γ(1))|e  = (γ
(0)

(0)

(C6)

− SL )|e .
(1)

(C7)

For the following perturbation calculation, it is convenient tointroduce the zeroth- and first-order left eigenvectors e(0)† |and e(1)† |, respectively. These satisfy

+

0=−

(C4)

q

Notice that, in the absence of the cytoskeleton (ν = 0),
˙ q + SL |q , which was disEq. (C2) reduces to 0 = |cussed by Seifert and Langer for vanishing tension σ = 0 [6].
Let us consider the eigenvalue problem(B10)

The Fourier transform of Eq. (8) relates ∂ρ/∂t = (∂/∂t)(ρ −ρ − )/2 with the longitudinal velocity v as ∂ρ/∂t = −iq(v+ −v− )/2. Then Eq. (B10) yields

σ q + κ̃q 3
kq|h; q h; q | −|h; q ρ; q | ,
4η
2η
(C3)

±

p
z e∓qz = [ḣ ± (ḣ ∓ v± )qz]e∓qz , (B7)
2η

p±
1 p±
∓ q ḣ +z e∓qzV = −iq 2η
2η
Vz = ḣ +where ĥ = h/d. Then Eqs. (B9) and (B11) [or Eqs. (9) and(10) without the force terms uh and uρ ] are written ase(0)† |SL = e(0)† |γ (0) ,(C8)e(0)† |(h − γ (1) ) = e(1)† |(γ (0) − SL ).

(C9)

Here γ (0) and γ (1) are common to Eqs. (C6) and (C7),respectively. Note that e(0)† | (e(1)† |) is not the Hermitian
conjugate of |e(0)  (|e(1) ), because neither SL nor h is a
Hermitian operator.
1. Zeroth order

Since SL does not couple the Fourier modes of different
wave vectors, we readily obtain the zeroth-order eigenvalues,q[(σ + κ̃q 2 )c(q)γ±(0) =
8ηc(q)
+ 4ηkq ± g(q)2 + 32ηk 2 d 2 q 3 c(q)] (C10)with g(c) = (σ + κ̃q 2 )c(q) − 4ηkq. Their associated righteigenvectors are given by(0)|e±; q = |h; q + e± (q)|ρ; q,(C11)
where e± (q) = [g(q) ∓ g(q)2 + 32ηk 2 d 2 q 3 c(q)]/[4kc(q)].
The corresponding left eigenvectors areq ∈Qq

012416-6(0)†e± ; q| =1†
1 + e± (q)e± (q)†[h; q| + e± (q)ρ; q|],

(C12)

DYNAMICS OF A BILAYER MEMBRANE COUPLED TO A . . .

PHYSICAL REVIEW E 96, 012416 (2017)

TABLE II. Approximate expressions of zeroth-order eigenvalues and eigenvectors for sufficiently small bare tension, σ(i) bσ/(2kη)(ii) q  qc

qqc

γ+(0) (q)γ−(0) (q)e+ (q)e− (q)

kq 2 /(2b)κ̃q 3 /(4η)(σ q + κq 3 )/(4η)
kκq 2 /(2bκ̃)−ηq/b−2kηd 2 q/(bκ̃)(qd)2
κ̃q 2 /(2k)†with e± (q) = c(q)e± (q)/(2ηd 2 q 3 ). In the above, the left
(0)†(0)
eigenvector is normalized such that e± ; q|e±; q = 1. Since
a contraction of left and right eigenvectors associated withdifferent eigenvalues vanishes, we obtain
 (0)†
e ; q |e(0) ; q = δ δq q .
(C13)
†

σc .

Similarly, we can calculate the first-order left eigenvector(1)†e± ; q| from Eq. (C9). One can also show that Eq. (C13)is generalized to the first order in ν ase(0)† ; q | + νe(1)† ; q | |e(0) ; q  + ν|e(1) ; q = δ δq q + O(ν 2 ).

†This yields e+ (q)e− (q) = e− (q)e+ (q) = −1, which can also†
be confirmed directly from the definition of e± and e± .
In Table II, we present approximate expressions of γ±(0) ande± for different length scales classified by the characteristic
wave numbers bσ/(2kη) and qc = 2ηk/(bκ̃), where a small
bare tension σ
σc is assumed [19]. Typical parametervalues quoted in the main text yield σc ≡ (2ηk)2 /(κ̃b2 ) ≈3 × 10−6 J/m2 . The behavior in the two regimes (i) and (ii)in Table II is essentially the same as those in Ref. [6] forσ = 0. In the presence of tension, there appears another regimeqbσ/(2kη), where the dynamics is again dominated by theintermonolayer friction [19]. For σ = 10−11 N/m chosen in
the main text, the characteristic wave number corresponds
to the length 4π kη/(bσ ) ≈ 0.44 m, which is too large to bemeasured in experiments. However, with σ = 4 × 10−7 N/m,
for instance, we have 4π kη/(bσ ) ≈ 10 μm, which is relevantfor giant unilamellar vesicles. For large tension σ  σc , on theother hand, the dynamics is dominated by the intermonolayer
friction in all length scales [19].
2. First order

We assume for simplicity that the zeroth-order eigenvalue
γ±(0) (q) is not degenerated in the Hilbert space Sq . This
assumption is always valid for scales much larger than thelattice spacing of the cytoskeleton, qa −1 . Equation (C7)(1)
yields the first-order correction νγ± (q) to the zeroth-ordereigenvalue γ±(0) (q) as
Kq
(0)†
(0)
.
; q = √γ±(1) (q) = e± ; q|h |e±†
4 gηq[1 + e± (q)e± (q)](C14)
To obtain the shifted eigenvectors, we expand the first-ordereigenvectors in terms of the zeroth-order eigenvectors,
 
(1)
|e±
; q =s± (q ; q)|e(0) ; q .
(C15)
q ∈Qq =+,−(C17)

Case (i): q  2ηk/(bκ̃)

Let us suppose a wave vector q L satisfies qL = |q L |2ηk/(bκ̃), i.e., case (i) in Table II. We then study how thecytoskeleton alters the rates γ±(0) (qL ). We notice that the zerothorder eigenvalue γ±(0) (qL ) is not degenerated in Sq L , because,
for such a small wave vector, there is no other wave vector q ∈Qq L satisfying |q | = qL . Using Eq. (C14) and approximate
expressions in Table II, we find
νbd 2 Kq L
,
√
4 gη2

γ+(1) (q L )

νKq
√ L ,
4 gηqLγ−(1) (q L )(C18)when qL2ηk/(bκ̃). In general, K(q L ) depends on the
direction of the wave vector q̂ L . However, it is approximatedas
Kq L

3(qL a)43(qL a)2
−
,
2
32

(C19)

for qLa −1 , which is isotropic. We can always use thisapproximation in Eq. (C18), which is valid for qL
2ηk/(bκ̃).
This is because 2ηk/(bκ̃) = 4.37 × 106 m−1 for typicalparameter values chosen in the main text, and it is smaller
than the reciprocal of the lattice spacing a ∼ 10−7 m ofthe cytoskeleton. Substitution of Eq. (C19) into Eq. (C18)
yields Eq. (11). Note for ν = 10−6 N/m chosen in the main
text, γ−(0) + νγ−(1) is larger than γ+(0) + νγ+(1) , and therefore, bydefinition, γ+ γ−(0) + νγ−(1) and γ− = γ+(0) + νγ+(1) .
(1)
Next we discuss the shifted eigenvector |e±; q L , qL
2ηk/(bκ̃). In the present case of qL
2ηk/(bκ̃), we canassume q  2π/a  2ηk/(bκ̃) for ∀ q ∈ Qq L \ {q L } (here“\” indicates set difference). Thus, for q = q L , we canset γ±(0) (qL ) − γ(0) (q ) −γ(0) (q ) in the denominator of
Eq. (C16). Using the expressions in (ii) of Table II, we obtainfor q = q L ,Here we set s±± (q; q) = 0 (as for the perturbation theory inquantum mechanics). Operation of e(0)† ; q | to the both sidesκ −s (q ; q L )2kd 2 of Eq. (C7) yieldss+ (q ; q L )Kq−√ L 4 .
g κ̃q

(C20)

For q = q L , we can assume γ+(0) (qL )  γ−(0) (qL ) in the
denominator of Eq. (C16) and obtain

(0)
; q
e(0)† ; q |h |e±s± (q ; q) = (0)(0)
γ± (q) − γ (q )†∓ [1 + e± (qL )e± (qL )]s∓± (q L ; q L )

Kq

= √
.
†
4 gηq [γ±(0) (q) − γ(0) (q )][1 + e (q )e (q )](C16)E(q L ),

(C21)

√
√
where E(q L ) = Kq L /[4 gηqL γ+(0) (qL )] bKq L /[2 gηkqL3 ].
To make the physical meaning of the shifted eigenvectors clear,

012416-7

OKAMOTO, KOMURA, AND FOURNIER

we examine the time evolution of the vector[ĥ(q ,t)|h; q  + ρ(q ,t)|ρ; q ].
|(t) =PHYSICAL REVIEW E 96, 012416 (2017)

the expressions in (ii) of Table II, we find(C22)

We may suppose that the modes for ∀ q ∈ Qq L \ {q L } decay
much faster than the modes for q L , because in the zeroth-order
γ±(0) (q ) is much larger than γ±(0) (qL ). Then, after sufficiently
large time t satisfying t  γ±(0) (q )−1 for ∀ q ∈ Qq L \ {q L },the modes for ∀ q ∈ Qq L \ {q L } are in the quasiequilibrium state ĥqe (q ; ĥ(q L )) and ρqe (q ; ĥ(q L )) determined by
Fq [δF /δh] = Fq [δF /δρ ± ] = 0 for a given value of ĥ(q L ).
Then we may approximate Eq. (C22) as|(t)

ĥ(q L ,t)|h; q L  + ρ(q L ,t)|ρ; q L ĥqe (q ; ĥ(q L ,t))|h; q ++ ρqe (q ; ĥ(q L ,t))|ρ; q  .

(C23)This will be justified later [see discussion around Eq. (C30)].
Using Fq [δF /δρ] = 2k[ρ(q ) − dq 2 h(q )] and Eq. (A4), we
obtain Eq. (12) at linear order in ν. We can then show that thefollowing relation holds up to the first order in ν:A (ĥ(q L ,t),ρ(q L ,t)) |e(0) ; q L +ν|e(1) ; q L  ,

=+,−

(C24)

k 2 d 2 Kq S.
√
gbκ̃ 2 qS2

νγ+(1) (q S )γ+(0) (qS )νγ−(1) (q S )

νKq S
,
√
g κ̃qS4(C28)

νkd 2 Kq S
νKq
∼ √ S4 .
√4g κ̃κqSg κ̃qS(C29)

Therefore, √the effect of cytoskeleton is negligible if (aqS )4 
2
2a νKq S /( 3κ̃) ≈ 0.0722Kq S . Since Kq  10, the correction
to the rates due to the cytoskeleton is not very large. However,
the associated modes do not necessarily decay to zero fort  1/γ±(0) (qS ). In fact, if a large-scale Fourier mode exists
in Qq S , i.e., q L ∈ Qq S (or equivalently q S ∈ Qq L ), the Fouriermodes for q S are rather excited to the quasiequilibrium state by
the large-scale Fourier mode ĥ(q L ), as discussed above. We can
(0)†
(1)†
confirm this explicitly by operating e± ; q S | + νe± ; q S | toEq. (C22). Owing to Eq. (C17), the resultant quantity decays
with the rate γ±(0) (qS ) + νγ±(1) (q S ) ∼ γ±(0) (qS ). Hence, for t 
γ±(0) (qS )−1 , we have
(0)†

(1)†

[e± ; q S | + νe± ; q S |]|(t)

where
A+ (ĥ(q L ),ρ(q L )) =
1
ν Ē(q L )e− +ĥ(q L )e− − e+e− − e+†νe Ē(q L )ρ(q) ,(C25)− 1− −
e− − e+and†
νe+ Ē(q L )1
A− (ĥ(q L ),ρ(q L )) =1−
ρ(q L )
e− − e+e− − e+ν Ē(q L )ĥ(q L ) ,− e+ +(C26)
e− − e+with Ē(q L ) = 2ηd 2 qL3 E(q L )/c(qL ) and e± = e± (qL ). Therefore, after sufficiently large time t, A+ and A− decay with therates in Eq. (11). The above discussion indicates that the rates
in Eq. (11) and their associated eigenvectors correspond to therelaxation of the (long-wavelength) modes of q L accompanied
by instantaneous relaxation of the other (short-wavelength)
modes of ∀ q ∈ Qq L \ {q L } to the quasiequilibrium state. Notethat even though at initial time we set ĥ(q ) = ρ(q ) = 0 forq ∈ Qq L \ {q L }, after sufficiently large time (t  γ±(0) (q )−1 )they are excited by nonzero ĥ(q L ) to their quasiequilibrium
values in Eq. (12).

(C27)

We also examine the magnitude of γ±(1) compared with that of
γ±(0) . Using the approximate expressions in (ii) of Table II, we
obtainγ−(0) (qS )q ∈Qq L \{q L }|(t)Kq(1)√ S , γ− (q S )4 gηqSγ+(1) (q S )q ∈Qq L

0.

(C30)

Using the approximations in Table II, we can show thatEq. (C30) is equivalent to the quasiequilibrium condition
Eq. (12) with q = q S and q = q L .
APPENDIX D: EFFECTS OF FRICTION BETWEEN THE
BILAYER AND THE PROTEINS AT THE
CYTOSKELETON VERTICES

We can take into account the friction between the bilayer
and the proteins at the vertices of the cytoskeleton bymodifying the lateral force balance Eq. (6) toδF+ ∂j τij± ± Tiz± ∓ b(vi+ − vi− )0 = − ∂iδρ ±

−λvi± δ(x − R ).
(D1)
In the above, the friction coefficient λ is assumed to be common
forboth the upper and the lower monolayers. The new termλ  vi± δ(x − R ) also breaks the translational symmetry andhence leads to the coupling between different Fourier modes.
With this new term, Eq. (B11) is modified to
δF
2
δF
∂ρ(q)
− −
c(q)
+ iqFq
iq
∂t
δρ +
δρ

(v + − v − )δ(x − R ) .
+ λq̂ · Fq

0=−Case (ii): q  2ηk/(bκ̃)

Next we discuss how the rate γ±(0) (qS ) with qS = |q S | 2ηk/(bκ̃) is altered by the cytoskeleton. In this regime, using
012416-8

(D2)

DYNAMICS OF A BILAYER MEMBRANE COUPLED TO A . . .

PHYSICAL REVIEW E 96, 012416 (2017)

With the use of the identity Eq. (A3), the new term is rewrittenas
Fq(v + − v − )δ(x − R )The correction to the rate γ±(0) (q) due to the friction is thengiven by(0)†

(0)
D± (q) ≡ − λe± ; q|ρ SL |e±
; q
†



2  1
=− √
[q ρ̇(q ) + q ⊥ w⊥ (q )],i g q ∈Q q 2

(D3)

q

+
−
− v⊥
)/2.
w⊥ = −iq(v⊥To eliminate w⊥ fromwhere
Eq. (D2), we need the transverse part of Eq. (D1). As in asimilar way to derive Eq. (D2), we obtain2cw (q)w⊥ (q)iq
(v + − v − )δ(x − R ) ,+ λq̂ ⊥ · Fq

0=−

(D4)



with cw (q) = 2b + ηq + μq . Equations (B9), (D2), and (D4)
are the complete set of the relaxation equations for ρ and ĥ.
We can see from Eq. (D3) that the Fourier modes for q are
coupled to the Fourier modes for ∀ q ∈ Qq , as in the casewithout the friction at the vertices of the cytoskeleton.
The full equations are also represented in terms of operatorsand vectors in the Hilbert space Sq , as in the previous section.
To perform perturbation calculations, we regard both ν and λas small parameters. Then, to the first order in ν and λ, thegoverning equation is2∂|(t)q = −(SL + νh − λρ SL )|(t)q ,
∂t

(D5)

where the operator ρ is defined as
1   q ·qρ (q) = √|ρ; q ρ; q |.
g q ∈Q q ∈Q cw (q )q 2
q

=− √

(D6)

q

[1] D. Mizuno, C. Tardin, C. F. Schmidt, and F. C. MacKintosh,
Science 315, 370 (2007).
[2] M. Guo, A. J. Ehrlicher, M. H. Jensen, M. Renz, J. R. Moore,
R. D. Goldman, J. Lippincott-Schwartz, F. C. MacKintosh, and
D. A. Weitz, Cell 158, 822 (2014).
[3] A. S. Mikhailov and R. Kapral, Proc. Natl. Acad. Sci. USA 112,
E3639 (2015).
[4] F. Brochard and J. F. Lennon, J. Phys. (Paris) 36, 1035 (1975).
[5] A. Zilker, H. Engelhardt, and E. Sackmann, J. Phys. (Paris) 48,
2139 (1987).
[6] U. Seifert and S. A. Langer, Europhys. Lett. 23, 71 (1993).
[7] J. Prost and R. Bruinsma, Europhys. Lett. 33, 321 (1996).
[8] T. Betz, M. Lenz, J.-F. Joanny, and C. Sykes, Proc. Natl. Acad.
Sci. USA 106, 15320 (2009).
[9] 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. USA 107, 1289
(2010).
[10] R. Rodríguez-García, I. López-Montero, M. Mell, G. Egea,
N. S. Gov, and F. Monroy, Biophys. J. 108, 2794 (2015).

λγ±(0) (q)e± (q)e± (q)†gcw (q)[1 + e± (q)e± (q)]

.

(D7)

To measure the relevance of D± (q), we shall consider the ratioq|D± /γ±(0) |. Using the approximations in (i) bσ/(2kη)
2ηk/(bκ̃) of Table II, we obtain D+ (q) 
λ
λ


(D8)
 (0)  2√gb ∼ ba 2 ,γ+ (q) D− (q) 
λd 2 q
λkd 2
λ
λ


∼ √ ∼ 2 . (D9)√
 (0)  2√gη
g
κ̃b
gb
ba
γ− (q)
Similarly, with the approximations in (ii) q  2ηk/(bκ̃) in
Table II, the ratio |D± /γ±(0) | is comparable to or smaller
than λ/(ba 2 ). Therefore, the friction due to the network is
negligible as long as the friction coefficient per area λ/a 2 ismuch smaller than the coefficient b for the intermonolayer
friction. We use Saffman-Delbrück theory to estimate thevalue of λ [23]. In this theory two physical situations are
examined: (i) membrane of finite size with surrounding fluid
being neglected and (ii) membrane of infinite size withsurrounding fluid being taken into account. Since λ is a “bare”
friction constant between a protein and a monolayer withsurrounding fluid being neglected, we may use the resultfor (i). Then we can set λ 4π μ/ ln(L/r0 ), with L themembrane size and r0 the protein size. Using L/r0 102 ,b 2 × 108 J s/m4 , μ 2 × 10−9 J s/m2 , and a 10−7 m,
we estimate λ/(a 2 b) 5 × 10−3 . We can thus neglect theeffects of λ. Note that the estimation with (ii) also leads tothe same conclusion. For (ii), we set λ 4π μ/ ln[μ/(ηr0 )]and obtain λ/(a 2 b) 1.8 × 10−3 when r0 2 × 10−9 m.

[11] K. Yasuda, S. Komura, and R. Okamoto, Phys. Rev. E 93, 052407(2016).
[12] 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).
[13] J.-B. Fournier, N. Khalifat, N. Puff, and M. I. Angelova, Phys.
Rev. Lett. 102, 018102 (2009).
[14] N. Gov, A. G. Zilman, and S. Safran, Phys. Rev. Lett. 90, 228101(2003).
[15] J.-B. Fournier, D. Lacoste, and E. Raphaël, Phys. Rev. Lett. 92,
018102 (2004).
[16] C. Dubus and J.-B. Fournier, Europhys. Lett. 75, 181(2006).
[17] G. Popescu, T. Ikeda, K. Goda, C. A. Best-Popescu, M. Laposata,
S. Manley, R. R. Dasari, K. Badizadegan, and M. S. Feld, Phys.
Rev. Lett. 97, 218101 (2006).
[18] L. Miao, U. Seifert, M. Wortis, and H.-G. Döbereiner, Phys.
Rev. E 49, 5389 (1994).
[19] J.-B. Fournier, J. Nonlinear Mech. 75, 67 (2015).

012416-9

OKAMOTO, KOMURA, AND FOURNIERPHYSICAL REVIEW E 96, 012416 (2017)[20] R. Okamoto, Y. Kanemori, S. Komura, and J.-B. Fournier, Eur.
Phys. J. E 39, 52 (2016).
[21] T. V. Sachin Krishnan, R. Okamoto, and S. Komura, Phys. Rev.
E 94, 062414 (2016).

[22] A. A. Evans, B. Bhaduri, G. Popescu, and A. J. Levine, Proc.
Natl. Acad. Sci. USA 114, 2865 (2017).
[23] P. G. Saffman and M. Delbrück, Proc. Natl. Acad. Sci. USA 72,
3111 (1975).

012416-10