PHYSICAL REVIEW E 89, 022713 (2014)

Diffusion coefficients in leaflets of bilayer membranes
Kazuhiko Seki
NRI, National Institute of Advanced Industrial Science and Technology (AIST), AIST Central 5,
Higashi 1-1-1, Tsukuba, Ibaraki 305-8565, Japan

Saurabh Mogre
Indian Institute of Technology Bombay, Powai, Mumbai, Maharashtra 400 076, India

Shigeyuki Komura
Department of Chemistry, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 2 December 2013; published 19 February 2014)
We study diffusion coefficients of liquid domains by explicitly taking into account the two-layered structure
called leaflets of the bilayer membrane. In general, the velocity fields associated with each leaflet are different andthe layers sliding past each other cause frictional coupling. We obtain analytical results of diffusion coefficientsfor a circular liquid domain in a leaflet, and quantitatively study their dependence on the interleaflet friction. We
also show that the diffusion coefficients diverge in the absence of coupling between the bilayer and solvents, evenwhen the interleaflet friction is taken into account. In order to corroborate our theory, the effect of the interleaflet
friction on the correlated diffusion is examined.
DOI: 10.1103/PhysRevE.89.022713

PACS number(s): 87.16.D−, 82.35.Lr, 68.05.−n

I. INTRODUCTION

Diffusion of inclusions in membranes occurs in twodimensional (2D) media embedded in solvents, and is different
from that in a homogeneous three-dimensional (3D) fluid.
In the pioneering work by Saffmann and Delbrück [1,2],
the membrane was regarded as a thin plane sheet boundedon both sides by solvents. Solvents are dragged by theflow in the membrane through the coupling to the solvent
flow at the solvent-membrane interfaces. The transfer ofmomentum between the membrane and solvents necessitates
the introduction of a length called the Saffmann-Delbrück (SD)hydrodynamic screening length [1–3]. The SD expression of
the diffusion coefficient has been elaborated to apply to thelarger size of objects [4]. Although some experiments were notexplained by the SD theory [5], the SD theory and its extended
expressions have been supported by certain experiments [6–8].
When the membrane is supported on a solid substrate, avery thin layer of solvent exists between the membrane and thesolid support [9]. For supported membranes, the assumptionof infinite thickness of solvents in SD theory is irrelevant; thediffusion coefficient decreases when the thickness of solvent
layer is decreased [9,10]. The length scale of momentum
transfer is characterized by a new length scale, called theEvans-Sackmann (ES) hydrodynamic screening length, in thelimit of thin layer of solvent. When the size of the diffusing
object exceeds the ES hydrodynamic screening length, thediffusion coefficient is predicted to show the inverse squaresize dependence [9]. The inverse square size dependence has
been observed in supported lipid bilayers [11].
In the earlier works, diffusion of solid molecules in homogeneous membranes has been studied. Recently, microheterogeneity of the membrane has attracted great interest [12–15].
The microdomains in multicomponent membranes have been
expected to be involved in signal transduction and control
intracellular transport [12–15]. By phase separation of ternary
mixtures of cholesterol, saturated and unsaturated lipids,
1539-3755/2014/89(2)/022713(12)regions rich in saturated lipids and cholesterol, coexist withthose rich in unsaturated lipids. The smaller regions form
circular viscous domains and diffuse laterally. The diffusion
coefficients of liquid domains have been obtained usinghydrodynamic calculations by ignoring a shear flow acrossthe membrane [16–20]. However, it is possible to form
membranes with bilayer consisting of two leaflets having
different compositions [14,15]. The liquid domains in eachleaflet do not necessarily lie on top of one another; theliquid domains can be formed in one of the leaflets without
overlapping to those in the other leaflet. In general, the velocity
fields associated with each leaflet are different and the layers
sliding past each other cause frictional coupling [21–25].
When relative motion of the leaflets does not play asignificant role, the bilayer can be considered as a single
fluid medium. However, the interleaflet friction would induce different transport in two leaflets of the bilayer forsupported membranes, because each leaflet faces different
environment [24]. Furthermore, when the solvent thickness issmall, typical for supported membranes, the relative velocity
difference between the wall and the membrane induces a shear
flow in the solvent through the stick boundary condition. The
resultant drag force from the solvent can be sufficiently large
to decouple the two leaflets.
The diffusion of solid particles in a leaflet under theinterleaflet friction was previously studied numerically [24].
The complex interplay between the interleaflet friction andthe drag from solvents was shown. In this paper, we derive
analytical results of diffusion coefficients for circular liquiddomains. Although we use the same mobility tensor derived
by Brown et al. [24], we present the analytical expressions
showing the effect of the interleaflet friction on the diffusion
coefficients for the general thickness of solvent. By examining
the asymmetry in diffusion coefficients associated with theupper and lower leaflet for supported membranes, we find thatthe value of interleaflet friction can be estimated. In addition,
we show that the momentum is dissipated by the interleaflet

022713-1©2014 American Physical Society

KAZUHIKO SEKI, SAURABH MOGRE, AND SHIGEYUKI KOMURA

friction, but this dissipation mechanism is not sufficient toovercome logarithmic long-range correlations in the flow fields
leading to a divergence of the diffusion coefficient. In order
to corroborate our theory, the correlated diffusion is studied
under the influence of the interleaflet friction.
In experiments, the frictional interaction between theleaflets has been studied in detail by applying the interlayer
shear stress induced by a large local deformation [26–29].
Zhang and Granick concluded that the difference in diffusion
coefficient of lipids between two leaflets is small for supported
bilayers [30]. Detailed investigation on diffusion of lipids inleaflets on supported membranes suggests that the diffusion
coefficient of lipids in the lower leaflet is reduced from that of
the upper leaflet depending on the type of substrate [31–33].
Simulations have been also performed to estimate the valueof interleaflet friction [34–36]. From coarse-grained molecularsimulations, the lateral diffusion of lipids in the lower leaflet tothe solid substrate was shown to be slower than that of the upper
leaflet by one order of magnitude [36]. We discuss the resultson the asymmetry in diffusion coefficients for the supported
membranes using the values of the interleaflet friction thus
estimated.
In Sec. II, the diffusion coefficient of a circular domain
in a bilayer membrane is obtained by taking into accountthe interleaflet friction for arbitrary thickness and viscosityof solvents. In Sec. III, we study the effect of the interleaflet
friction on the SD length and ES length. Symmetric environments are considered for simplicity. Diffusion in a supported
bilayer is considered in Sec. IV. In Sec. V, we study correlated
dynamics of two point particles. In Sec. VI, we show that thefinite diffusion coefficient cannot result from the interleaflet
friction alone without solvents. The final section is devoted forsummary and discussion.
II. DIFFUSION COEFFICIENT OF A CIRCULAR
DOMAIN IN A LEAFLET
A. Hydrodynamic model

We explicitly take into account the two-layered structure of
the membrane as shown in Fig. 1. The structure is sometimes
referred to as the leaflet and its hydrodynamics has beenstudied previously [21,24]. These leaflets are denoted by theindices “+” and “−” for the upper and lower monolayer,
respectively. The monolayers of the membrane, called leaflets,
are coupled through the drag force between the upper andlower leaflets.
We consider a circular liquid domain of radius R which isimmersed in either leaflet. For simplicity, we assume that theviscosity of the liquid inside the domain is the same as thatof the leaflet in which the domain is located. We also assume
that the domain maintains a circular shape by the line tension
and does not undergo deformation as it moves. The system is
illustrated in Fig. 1. We will derive the general expression forthe diffusion coefficient in either case.
The 2D viscosity of the leaflets are denoted by η± ,depending on the upper or lower leaflet. The bilayer membrane
is sustained in a solvent at distance H + and H − below andabove the walls, respectively. The viscosity of the solvent
above and below the membrane is represented by ηf+ and ηf− ,respectively.

PHYSICAL REVIEW E 89, 022713 (2014)

Wa ll

H+

η+
f
η+

H-

2R

ηf

η-

Wa ll

FIG. 1. (Color online) The geometry of the bilayer in a generalenvironment. Liquid domains of size 2R diffuse in leaflets.

The in-plane flow in each leaflet can be described using
Stokes equations with additional terms associated with thecoupling between the leaflet and the solvent through theinterface and the friction between leaflets [21–24],
η± ∇ 2 v ± − ∇p± + K ± ∗ v ± ∓ (v + − v − ) + f ± = 0, (1)with the incompressibility condition,∇ · v ± = 0,(2)where p± and v ± indicate the hydrodynamic pressure andfluid velocities of the leaflets, respectively. In the above,
the position vector r is abbreviated. The term (v + − v − )accounts for the frictional coupling between
the two leaflets.
The convolution term K ± ∗ v ± = d 2 r K ± (r − r )v ± (r )represents the drag from the solvent. Using Fourier transform K̂(k) = d 2 r exp (−ik · r) K(r), the convolution canbe expressed as a product. The drag force from the solvent
depends on the geometry of the system. For the system underconsideration, we have K̂ ± = −ηf± k coth(kH ± ) for arbitrary
solvent thickness [37–40]. The external in-plane force isdenoted by f ± .
B. Diffusion coefficient

We consider the case when the circular domain moves witha velocity U in leaflet α. The origin of the coordinates istaken to be the center of the domain and the x coordinate
is chosen in the direction of U. The external in-plane force
applied at the periphery of the circular domain in Eq. (1)is given by f ± = f (n) cos θ nδ(r − R)/(2π R) [16], where nis the outward normal unit vector on the circumference of thedomain, θ is the angle between U and r, and cos θ term reflects
the symmetry of the system. f (n) is a constant that should givethe external force after the integration over the circular domainboundary.

022713-2

DIFFUSION COEFFICIENTS IN LEAFLETS OF BILAYER . . .

PHYSICAL REVIEW E 89, 022713 (2014)

We can derive the diffusion coefficient of the circulardomain in either leaflet using a similar method as in ourprevious publication [18]. By introducing a component of amobility tensor in Fourier space given in Appendix A [24], thediffusion coefficient of the liquid domain is obtained as, ∞[J1 (kR)]2αDdom= kB Tdk
π kR 20η(−α) k 2 + ηf(−α) k coth(kH (−α) ) + 
×  
,

(α)
(α) 2
(α) ) +  − 2α η k + ηf k coth(kH

Wall

H

η

H

η

f

η

(3)
where α denotes the leaflet, −α denotes the leaflet opposite toα and J1 (z) denotes the Bessel functions of the first kind oforder 1 [41].
The expression for the diffusion coefficient contains various
physical quantities with varying magnitudes and different
units. In order to systematically compare various cases andto have a concise understanding of the results, we make thevariables dimensionless mainly using the SD hydrodynamic
screening length given by 1/ν ± ≡ η± /ηf± [1–3]. The following dimensionless variables are defined:
ρ ± = ν ± R,

h± = ν ± H ± , λ± = 1,
η± (ν ± )2

μ=

η
f

Wall

FIG. 2. (Color online) The geometry of the bilayer in a symmetric environment. Liquid domains of size 2R shown by red (dark gray)diffuse in leaflets.

η−
.
η+
(4)

Here, ρ ± and h± represent the dimensionless radius of
the liquid domain and the dimensionless solvent thickness,
respectively. The ES√hydrodynamic screening length is givenby the smaller of h± /ν. A new length, η− /, is associated with the interleaflet friction. This quantity can becalled the interleaflet sliding length, and contributes to thehydrodynamic screening factor in addition to the contribution
from the solvent. The variable μ is the ratio between the 3D
viscosities of the leaflets. In the later sections, the index of thedimensionless variables is dropped when it is redundant, e.g.,h = h± whenever it is clear from the context or the variables
do not differ between the indices.

We can see from Eq. (5) and Fig. 3 that the diffusion
coefficient decreases from D0 to D0 /2 as the interleaflet
friction is changed from 0 to ∞. The qualitatively similar
behavior is obtained by Camley and Brown for a solid circular
object in a freely suspended membrane when the symmetricleaflets are faced to the same solvent medium [24]. We also
note in Fig. 3 that the diffusion coefficient as a function of λbehaves very differently when ρ
1. The reason is separately
discussed in Appendix C.

1

D/ D 0III. DIFFUSION FOR A LIQUID DOMAIN
IN A SYMMETRIC ENVIRONMENT

For simplicity we consider the symmetric system as
illustrated in Fig. 2. We assume a solvent layer of equal
thickness on both sides of the membrane, i.e., H + = H − = H .
The viscosity of the solvents and two leaflets are also taken tobe equal, i.e., ηf+ = ηf− = ηf and η+ = η− = η.
For this case, the diffusion coefficient is simplified to,
Dsym = 12 [D0 + D1 (λ)],
where
D1 (λ) =kB Tπ ηρ 2

 ∞
dκ
0

(5)

1
[J1 (κρ)]2.
2κκ + κ coth(κh) + 2λ
(6)

In the above, D0 represents D1 (λ = 0), which is the case in
the absence of the interleaflet friction [18].

0.75

0.5
−2
10

10

−1

0

10

λ

1

10

10

2

FIG. 3. Normalized diffusion coefficient D/D0 as a functionof λ for h = 1 is shown by using Eq. (5). Here, D0 = D(λ = 0).
The values of π ηD0 /(kB T ) are 0.5882(ρ = 0.1), 0.1404(ρ = 1),
and 0.0044(ρ = 10). The solid, dashed, and dotted lines represent
ρ = 0.1,1,10, respectively.

022713-3

KAZUHIKO SEKI, SAURABH MOGRE, AND SHIGEYUKI KOMURA
A. Infinite solvent medium

In the limit of H → ∞, we have coth(κh) ≈ 1 in Eq. (6).
The expression for the diffusion coefficient becomes
 ∞
1
1
kB T
[J1 (κρ)]2+ 2.
dκ
Dsym =
22π ηρ 0κ
κ2 + κ
κ + κ + 2λ
(7)
The above integral can be calculated by using MATHEMATICA
and is expressed in terms of Meijer G functions [42,43],
 1 1 

1
1
kB T
,
32
2
G2 4 (c1 ρ)  12 2
Dsym =
√
2
3/2
4π ηρ
0, 2 ,1,−1
1 − 8λ π c1
 1 1 
1,
− 3/2 G32 24 (c2 ρ)2  12 2
π c2
0, 2 ,1,−1
 1 3 

1
1
,
32
2
,
(8)
+
− 1 − 3/2 G2 4 ρ  2 3 22λπ0,1, 2 ,−1where c1 and c2 are the roots of the quadratic equation, κ +κ + 2λ = 0. The limits for the domain size can now be taken
to this expression.
By taking the limit of small domain size, i.e., ρ1, c1 ρ1, and c2 ρ
1, in Eq. (8), the diffusion coefficient is shownto depend logarithmically on the domain size,

2η
1
kB T
1
Dsym ≈ln−γ + +
4π η
ηf R
4 2(c1 − c2 )PHYSICAL REVIEW E 89, 022713 (2014)

We can analytically integrate Eq. (11) and the result is given
by,




ηf
ηf
kB T H 1− I1 R
K1 R
Dsym =
2π R 2 ηf 2
ηH
ηH


(ηf /H ) + 21
1
− I1 R+(ηf /H ) + 2 2η 
(ηf /H ) + 2.
(12)
× K1 Rη√
By taking the limit of ρ
1, ρ/ h1 and√
1 in Eq. (12), we obtain,
ρ/ (1/ h) + 2λ 
2 ηf H1
kB T
Dsym ≈ln
4π ηRη [1 + (2H /ηf )]1/4−γ +2× (c2 ln c2 − c1 ln c1 ) ,(9)where γ = 0.5772 . . . is Euler’s constant [41].
By taking the limit of large domains, i.e., c1 ρ
1, ρ1and c2 ρ
1, in Eq. (8), we obtain,


1
2kB T
kB T
2η
Dsym ≈.
(10)+−
3π 2 ηf R
8π R 2  ηf2
In the asymptotic limit of large domain size, the diffusion
coefficient is independent of the membrane viscosity η. The
asymptotic domain size dependence is given by the first term,
1/R, which is independent of the interleaflet friction. The
second term showing 1/R 2 -dependence is associated with theinterleaflet friction.
Equation (10) implies that the diffusion coefficient dependson λ if it is smaller than λ∗sym = 3π/(16ρ). By plotting the
diffusion coefficient as a function of λ as shown in Fig. 3,
we confirm that the inflection point where Dsym depends moststrongly on λ, can be estimated by using λ∗sym when ρ  0.1.

Dsym ≈1
κ 2 + 1/ h + 2λ+1
κ 2 + 1/ h

.

(11)

kB T H2π ηf R 2

1 + H /ηf
1 + 2H /ηf

.

(14)

As in the case of an infinite solvent medium, the diffusion
coefficient is independent of the membrane viscosity. The
diffusion coefficient shows 1/R 2 dependence in this limit,
which is stronger than that in Eq. (10). The size dependenceof the diffusion coefficient is qualitatively the same as thatobtained by regarding the membrane as a sheet of 2D fluid
without any leaflet structure, though there are quantitative
differences given by H /ηf dependence.
C. Arbitrary thickness of solvent medium

The integrand in Eq. (6) can be expressed by a sum using
the partial fraction expansion as shown in Appendix B. The
result after the integration is given as
Dsym =∞
1
kB T 
Ws
− I1 (ωs ρ)K1 (ωs ρ)22π ηρ s=12
1
− I1 (ωs ρ)K1 (ωs ρ) ,
+ Ws
2

(15)

where the weight factor is given by

B. Thin solvent medium

×

(13)

As in the case of an infinite solvent medium, Eq. (9), thediffusion coefficient is shown to depend logarithmically on
the domain size.
√
By taking the limit of ρ
1, ρ/ h1 and√ρ/ (1/ h) + 2λ
1, Eq. (12) reduces to,Ws =In the limit of, H → 0, we have coth(κh) ≈ 1/κh in Eq. (6).
By taking this limit, the expression for the diffusion coefficient
becomes ∞kB T[J1 (κρ)]2dκ
Dsym =
22π ηρ 0κ

1.
42hωs4 + (1 + h − 4λh)ωs2 + 2λ + 4hλ2

(16)and ωs satisfies the characteristic equation,ωs2 = ωs cot(ωs h) + 2λ.

(17)

In the above, Ws = 2/[ωs2 (1 + h + ωs2 h)] and ωs = cot(ωs h)
are Ws and ωs obtained by setting λ = 0, respectively.
Moreover, I1 (z) and K1 (z) are the modified Bessel functionsof order 1 [41]. We can show that the rigorous condition for
H → 0 to obtain Eq. (12) is h = νH < π . This is the condition

022713-4

DIFFUSION COEFFICIENTS IN LEAFLETS OF BILAYER . . .

PHYSICAL REVIEW E 89, 022713 (2014)

0

10

−1

πηD/ kB T

10

η

f

−2

10

η

−3

10

−4

10

−5

10 −2
10

0

10

ρ

2

10

η

H

FIG. 4. π ηD/kB T as a function of ρ for h = 1 is shown by using
Eq. (5). The solid, dashed, and dotted lines represent λ = 0.1,1,10,
respectively. The dash-dotted line indicates the asymptotic 1/ρ 2
dependence.

to obtain Eq. (12) from Eq. (15) by retaining the first term in
the expansion and ignore the rest.
In the above, 1/(ωs ν) represents the hydrodynamic screening length under the interleaflet friction. The diffusion
coefficient is given by the infinite sum of the expression
associated with different hydrodynamic screening lengths.
In this sense, hydrodynamic screening is related to multiple
length scales given by 1/(ωs ν) with s = 1,2,3, . . .. The weightapproximately satisfies the relation Ws ∼ 1/ωs4 . Below, we
assume that 1/(ωs ν) is sorted from small to large values. The
difference between the successive screening length, 1/(ωs ν)
and 1/(ωs+1 ν), is approximately equal to π/ h. The weightfactor Ws rapidly decreases if h = νH < π . In this case,1/(ω1 ν) can be regarded as the virtual hydrodynamic screening
length. The difference due to retaining only the first term in
the diffusion coefficient is numerically investigated, and it is at
most 10% when h = 1. The difference increases by increasing
the value of h.
In Fig. 4 we show the diffusion coefficient as a functionof the domain size. The diffusion coefficient decreases by
increasing λ. When ρ is small, the variation of the diffusion
coefficients for the different values of the interleaflet friction
decreases. This can be understood by noticing that the diffusion
coefficient depends logarithmically on the strength of theinterleaflet friction.

IV. DIFFUSION IN A SUPPORTED BILAYER

In this section, we consider the case where asymmetry inthe system is introduced by changing the thickness of thesolvent layer above and below the membranes. The valuesof the viscosities are kept the same as in the symmetric casefor simplicity, although a similar calculation is possible withdifferent values as well. Figure 5 illustrates the system being
considered.

η

f

Wa ll

FIG. 5. (Color online) The geometry of the bilayer in an asymmetric environment. Liquid domains of size 2R shown by red (darkgray) diffuse in leaflets.
A. Formally exact results

We note that the domain under consideration can be locatedin either of the two leaflets and diffuses differently in each case.
We denote the diffusion coefficient when the domain is in theupper and lower leaflets by D + and D − , respectively.
The situation of supported bilayer can be mathematically
expressed as H + → ∞ and ν − H −1. In this limit, wehave coth(κh+ ) ≈ 1 and coth(κh− ) ≈ 1/κh− . By taking above
limits and using Eq. (3), the expression for the diffusion
coefficient becomes, ∞g α (κ)kB T[J1 (κρ)]2
α
,
Dsup
=
dκ
2
+
π ηρ 0κg (κ)g − (κ) − λ2(18a)where g α is given by,
1
+ λ,
h

(18b)

g − (κ) = κ 2 + κ + λ,(18c)g + (κ) = κ 2 +andrespectively.
As before, by applying a partial fraction expansion to
Eq. (18), we obtain [42]4kB T  g α (ωi ) −1αDsup=π ηρ 2 i=1 f  (ωi ) 2ωi 1 1 

1
,
32
2
,
(19)
× 1 − 3/2 G2 4 (ωi ρ)  12 2
0,
,1,−1
π
2
where ωi are the roots of the denominator, f (κ) =g + (κ)g − (κ) − λ2 = 0 and their real parts are negative, whilef  (κ) is the derivative of the denominator with respect to κ.
We note that ωi are independent of ρ, and use this fact inunderstanding the size dependence.

022713-5

KAZUHIKO SEKI, SAURABH MOGRE, AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 89, 022713 (2014)

As in the previous section, by taking limits of ρ in Eq. (19),1,
we obtain for ωi ρ4
2
1
kB T  ωi g α (ωi )
α
ln
−γ +
, (20)Dsup≈
4π η i=1 f (ωi )ωi ρ41
kB T  g α (ωi )2π η i=1 f  (ωi )
4

α
Dsup
≈

−1
ωi ρ 2πηDα/kB T

and for ωi ρ

0

10


.

(21)

−1

10

As in the previous case, the diffusion coefficient depends
logarithmically on the domain size for small domains, or the
inverse square of the domain size for large domains.
−2

10 −2
10

B. Approximate results

Although we have obtained the formal solution, Eq. (19),in terms of the roots of a fourth order polynomial equation,
finding the λ or h dependence of the roots is a formidable
task. In order to avoid this difficulty, we simplify the modelby taking the limit of ηf+ → 0. When the solvent is water, thesolvent viscosity is about two orders of magnitude smaller thanthat of lipids. The limit of ηf+ → 0 could be justified for a wide
range of the values of h, λ, and ρ. We thoroughly investigate
the mathematical structure of the diffusion constant in the limitof ηf+ → 0. Then, we compare the analytical solution in thelimit of ηf+ → 0 and that with a finite ηf+ .
By taking ηf+ = 0 in Eq. (3) and assuming that the solvent
layer below the membrane is very thin, i.e., ν − H −
1, we
get,
 ∞
g α (κ)kB T[J1 (κρ)]2
α
,
Dsup
=
dκ
−
2
+
−
πη ρ 0κg (κ)g (κ) − λ− λ+(22a)where we denoteg + (k) = κ 2 + 1/ h− + λ− ,g − (k) = κ 2 + λ+ . (22b)

For simplicity, we consider the case when the viscosity of theupper leaflet is the same as that of the lower leaflet.
The integral of Eq. (22a) can be evaluated analytically to
give,

dα
1
kB T d α
H
α
1
−
Dsup
=
+
π R 2 2b
ηf c + − c −
c−




νc−
νc−
K1 R
× I1 R
H
H





νc+
νc+
dα
K1 R
, (23a)− 1 − + I1 Rc
H
H
where we have introduced the notations,d+ = 1 +

H
,
ηf

and

d− =

H
,
ηf


H
1
1+2
c =
±
2
ηf
±

H
1+ 2
ηf

(23b)

2
,

(23c)

−1

10

10

0

λ

10

1

10

2

3

10

FIG. 6. π ηD α /kB T as a function of λ is shown for h = 0.1 andρ = 1. The circular and triangular markers denote D + and D − forthe supported case, while the solid and the dashed lines denote D +and D − in the limit of ηf+ → 0, respectively.

where the double sign corresponds. The first term shows 1/R 2
size dependence. It is independent of the membrane viscosity,and is a linear function of the inverse of the interleaflet friction
constant. It can be seen from Eq. (23) that when 1 ≈ 2H /ηf(1/ h ≈ 2λ), D changes rapidly by changing the interleaflet
friction coefficient. The diffusion coefficient as a functionof λ for a supported bilayer is shown in Fig. 6. Though
the argument is not rigorous, we see from Fig. 6 that the
diffusion coefficient indeed changes when λ is close to theabove mentioned value. This value is equal to λ∗2 obtained forthe domain in a symmetric environment (see Appendix C fora detailed discussion on λ∗2 ). By decreasing ρ, we note that the
inflection point is not characterized by λ∗2 but is close to λ∗1 asin the case of a symmetric environment. (see Appendix C fora detailed√discussion on λ∗1 ).
For R νcα /H
1, Eq. (23a) can be expanded to yield,  
2 ηf H  ηf 1/4
kB T
1αDsup ≈
−γ +
ln
4π η
R
η
H
4

ln(cα /c(−α) )+,(24)
4 1 + (2H /ηf )2
where (−α) represents the leaflet opposite to α as before. The
diffusion coefficient shows √a logarithmic size dependence in
1, cannot be satisfied
this limit. The condition, R νcα /Hin the  → ∞ limit because cα given by Eq. (23) diverges.
The order√of the limit is important in this respect.
For R νcα /H
1, the diffusion coefficient simplifies
into
kB T d α
α
.
(25)
≈
Dsup
2π R 2 This diffusion coefficient scales with the domain size as 1/R 2 .
−In this limit, Dsup
is independent of the value of the interleaflet
friction coefficient.

022713-6

DIFFUSION COEFFICIENTS IN LEAFLETS OF BILAYER . . .

As seen in Fig. 6, the λ dependence of the diffusion
coefficient for the supported case is well captured by thatobtained in the limit of ηf+ → 0. Both of these cases have alarge asymmetry between D + and D − as the friction between
the leaflets changes. This implies that contribution to friction
from the fluid above the bilayer is small, and the friction fromthe lower fluid is dominant. D − has a very low friction
dependence and does not change much by changing λ.
The λ dependence can be rationalized by taking the limits.
In the limit of λ → ∞, Eq. (22a) becomes,
 ∞
1
kB T
[J1 (κρ)]2
α
(26)
=
dκ
Dsup
2
2
2π ηρ 0κ
κ + 1/(2h)kB T H 1
νν=− I1 R
K1 R. (27)
2π ηf R 2 22H2HFigure 6 shows that for very large values of λ, D + , andD − converge to the same value for both the cases. Thisis because the leaflets are moving together and losing theirbilayer nature due to the high interleaflet friction. Whereas in+
the limit of λ → 0, Dsup
in Eq. (22a) shows a divergence, and−Dsup in Eq. (22a) reduces to the same expression as Eq. (27)with ν/2 replaced by ν. This replacement implies that themonolayer thickness should be used instead of the bilayer
thickness to obtain the membrane viscosity in the absence of
the interleaflet friction. Note that the 2D membrane viscosity
is its 3D viscosity multiplied by the membrane thickness. The−
slight change of Dsup
in Fig. 6 is the result of the change in thehydrodynamic screening length from 2/ν to 1/ν.
The diffusion coefficient as a function of ρ for a supported
bilayer is shown in Fig. 7. We note that the difference between
D + and D − increases for both cases with increasing ρ asshown in Fig. 7. This implies that smaller sized domains havelower asymmetry between the upper and the lower leaflets,
and it increases with increasing the size of the domain.

PHYSICAL REVIEW E 89, 022713 (2014)
V. CORRELATED DYNAMICS BETWEEN TWO LEAFLETS

So far, we have studied the diffusion of a liquid domain ineither side of the leaflet of bilayer. According to the strengthof the interleaflet coupling, the diffusion in a leaflet influences
the diffusion in another leaflet. In this section, we study thecorrelated diffusion of two point particles in different leaflets.
We define the x axis using the line connecting theinitial positions of two particles. By introducing the mobility
tensor [24,40]

ki kj
αβ
αβ
(28)
Gij (k) = M (k) δij − 2 ,kwhere M αβ (k) is given in Eq. (A2), we can calculate
αβ
the correlation between displacements, ri(1) rj(2) r =αβ
2kB T tGij (r), by the inverse Fourier transformation.
We consider the longitudinal and transverse coupling ofdisplacements for two particles. The longitudinal coupling difαβαβfusion coefficient is defined by DL (r) ≡ rx(1) rx(2) r /(2t)and is expressed as, ∞dk J1 (kr) αβαβ
M (k).
DL (r) = kB T
(29)
2π r
0
On the other hand, the transverse coupling diffusion coefficientαβ
αβ
is defined by DT (r) ≡ ry(1) ry(2) r /(2t) and is expressedas, ∞dk ∂J1 (kr) αβαβ
M (k),DT (r) = kB T
(30)
2π
∂r
0
where ∂J1 (z)/∂z = J0 (z) − J1 (z)/z should be noted [41].
Using appropriate function for M αβ (k), the above expressions
can be evaluated for all the cases discussed in the previous
sections.
As a relevant example, we calculate the coupling displacements for the supported case mentioned in Sec. IV. Using
the appropriate mobility tensor, we obtain the longitudinal
diffusion coefficient by using MATHEMATICA as [42]

0

10

αβ

πηDα /kB T

DL (r) =+ H−1 (−ωi ρ)} ,

−2

10

(31)

where hαβ (x) = g α (x)δαβ + λ(1 − δαβ ), and ωi , g α (x) andf (x) were defined in Eq. (18). Y1 (z) is the Bessel function of
the second kind of order 1, and H−1 (z) is the Struve H functionof order −1 [41]. Here, ρ denotes the distance between two
particles normalized by ν, ρ = νr. The transverse coupling
diffusion coefficient is calculated from Eq. (30) by changing
the order of integration and differentiation [42]

−4

10

−6

10 −2
10

4
kB T  hαβ (ωi ) −2
+ π {Y1 (−ωi ρ)
4π ηρ i=1 f  (ωi ) ωi ρ

0

10

ρ

kB T  ωi hαβ (ωi )8π η i=1 f  (ωi )
4

2

10

αβ

DT (r) =FIG. 7. π ηD α /kB T as a function of ρ is shown for h = 0.1 andλ = 1 for the supported case and in the limit of ηf+ → 0. D + and D −are shown by circular and triangular markers for the supported case,whereas they are shown by the solid and dashed lines in the limit ofηf+ → 0, respectively. 1/ρ 2 is shown by a dash-dotted line.
022713-7×2(2 + ωi ρ)− π {Y0 (−ωi ρ) − Y2 (−ωi ρ)(ωi ρ)2+ H−2 (−ωi ρ) − H0 (−ωi ρ)} .

(32)

KAZUHIKO SEKI, SAURABH MOGRE, AND SHIGEYUKI KOMURA(a)PHYSICAL REVIEW E 89, 022713 (2014)
(b)

0.25

0.08

0.2

0.06

πηD T /kB T

πηD L /kB T

0.1

0.15
0.1

0.04
0.02
0

0.05
−0.02
0
−2

0

10

10

10

λ

−0.04 −2
10

2

10

0

2

λ

10

FIG. 8. The (a) longitudinal and (b) transverse correlation as a function of λ calculated from Eq. (31) and Eq. (32), respectively. The solidline indicates the case that both particles are embedded in the upper leaflet, whereas the dashed line indicates the case that both particles areembedded in the lower leaflet. The dotted line indicates that two particles are embedded in different leaflets. Parameters used are h = 0.1,
ρ = 1.

Equations (31)–(32) are obtained for point particles and valid
without any restrictions on the separation distance. We study
the above equations by taking the limits for ρ. Equations (31)and (32) becomekB T
αβDL,T (r) ≈4ωi hαβ (ωi )4π η i=1

in the limit of ωi ρ
αβDL,T (r) ≈

−2
ln
ωi ρ


−γ ±

1
, (33)
2


4
kB T  ∓hαβ (ωi ),2π ηρ 2 i=1 ωi f  (ωi )(34)f  (ωi )1, andin the limit of ωi ρ
1. The upper and lower sign corresponds
to the longitudinal and transverse component, respectively.
The correlated diffusion coefficients as a function of thefrictional coupling are shown in Fig. 8. From Fig. 8(a),we note that the longitudinal correlated diffusion exhibits
a similar behavior as the diffusion coefficient if α = β. As
the value of the frictional coupling increases, the correlation
decreases when the particles are in the same leaflet, whereas
it increases when the particles are in the different leaflets.
It is seen that the correlated diffusion between two particles
in the different leaflets is small when λ  1. However, as
the value of the frictional coupling increases, the correlated
diffusion coefficient in the different leaflets tends to converge
to the smaller value of the correlated diffusion in the sameleaflet.
As shown in Fig. 8(b), the transverse components show
some interesting features. For particles in the same leaflet, thecorrelation reduces from a positive value to a negative value
for particles in the upper leaflet on increasing intermembrane
friction, while the correlation remains negative and approaches
the value of the upper leaflet when particles are in the lower
leaflet.

The correlated diffusion coefficients as a function of theseparation between two inclusions are shown in Fig. 9. The
longitudinal correlated diffusion exhibits a similar behavior as
the diffusion coefficient. The transverse component decreases
to a negative value and turns to increase as the separation
increases if α = α  . In the description of the bilayer as a
uniform 2D fluid, a similar nonmonotonic dependence ofthe coupling diffusion coefficients on the separation between
the inclusions has been reported [3]. The correlated diffusion
coefficient in different leaflets monotonically decreases to zeroby increasing the separation.

VI. SUMMARY AND DISCUSSION

We have investigated the effect of interleaflet friction on
the domain diffusion in a leaflet of a bilayer. We consideredboth symmetric and asymmetric environments of solvents. For
both cases, we have shown that the effect of the interleaflet
friction on the leaflet diffusion is important if λ is smaller
than λ∗2 = 1/(2h); the value is taken from the inflection
point of diffusion coefficient as a function of the interleaflet
friction for large domains (see Appendix C for details). For
small domains, the inflection point can be characterized byλ∗1 = [1/ρ − cot(h/ρ)]/(2ρ). The critical domain size canbe obtained from λ∗1 = λ∗2 . We have also studied the size
dependence of the leaflet-diffusion coefficient.
According to recent experiments, the value of interleaflet
friction is estimated to be  ∼ 108 –109 Pa s/m [26,29,33,44].
The relation, λ∗2 = 1/(2h), can be rewritten as ∗ = ηf /(2H ).
By substituting the typical values, ηf ∼ 1 × 10−3 Pa s andH ∼ 1 nm, we find ∗ ≈ 5 × 105 Pa s/m. However, if solventsare confined in molecular scales, the solvent viscosity isreported to be significantly increased, ηf ∼ 0.1 Pa s [24,30].
In the extremely confined situation, ∗ can be increased to 5 ×
107 Pa s/m. Experimental values of  ∼ 108 − 109 Pa s/m is022713-8

PHYSICAL REVIEW E 89, 022713 (2014)

(a)

(b)

1.2

1

1

0.8

0.8

0.6

πηD T /kB T

πηD L /kB T

DIFFUSION COEFFICIENTS IN LEAFLETS OF BILAYER . . .

0.6
0.4

0.4

0.2
0.2
0
0
−2
10

−1

10

10

ρ

0

1

10

10

2

10

−2

−1

10

10

0

1

10

ρ

2

10

FIG. 9. The (a) longitudinal and (b) transverse correlation as a function of ρ calculated from Eq. (31) and Eq. (32), respectively. The solidline indicates the case that both particles are embedded in the upper leaflet, whereas the dashed line indicates that both particles are embeddedin the lower leaflet. The dotted line indicates that two particles are embedded in different leaflets. Parameters used are h = 0.1, λ = 1.

not much different from the value of ∗ for the extremely
confined case. The value of λ∗1 relevant to the small size ofdiffusing objects can be order of magnitude larger than that ofλ∗2 . For small domains, by using λ∗1 we estimate that ∗ canbe 109 Pa s/m for the extremely confined case and 107 Pa s/mwithout molecular scale confinements.
According to simulations, the value of interleaflet friction
is reported to be  ∼ 2.4 × 106 Pa s/m [34,35]. There, theleaflet diffusion at small sizes can be affected by the interleaflet
friction. Indeed, the lateral diffusion of the lower leaflet was
shown to be slower than that of the upper leaflet by theattractive interaction of the support within a coarse-grained
molecular simulation [36].
In this paper, we have studied asymmetry in leaflet diffusion
coefficients induced by the asymmetry in environments. Byanchoring lipopolymer, the fluidity of a leaflet composing
a bilayer can be obstructed and the fluidity of leaflets canbe asymmetric [45,46]. When the fluidity of leaflets isasymmetric, the interleaflet friction should be in principle
taken into account. Although the model studied in this paper issimpler than the situation of lipopolymer-grafted bilayers, ourstudy could be useful to understand the effect of the asymmetry
in the fluidity between two leaflets on the leaflet diffusion.
Finally, we discuss the diffusion of a circular liquid domain
in an isolated membrane without any solvent to study the actual
contribution of the momentum dissipation by the interleaflet
friction. In the absence of the coupling to 3D solvent, the SD
length given by the ratio between the 2D membrane viscosity
and the 3D solvent viscosity diverges. This divergence reflects
the so-called Stokes paradox; the diffusion coefficient cannot
be obtained by solving hydrodynamic equation in the limit
of low Reynolds number in pure 2D fluid since the influenceof translational motion of a circular body extends to large
distances due to lack of enough momentum dissipation in two
dimensions. Although the coupling between solvents and thebilayer is absent, the momentum can be dissipated through theinterleaflet friction. In order to see whether the Stokes paradox
can be resolved by interleaflet friction alone, we set ηf± = 0 in
Eq. (3) resulting in,
kB T
D=
π R2

 ∞
dk
0

[J1 (kR)]2kη− k 2 + b.
k 2 (η+ η− k 2 + η+ η− )(35)

The above expression can be rewritten as
D=

 ∞
kB T
1
η−
[J1 (kR)]2
dk+
π R 2 η+ + η− 0
k3
η+

 ∞
1
[J1 (kR)]2.
×dkkk 2 + (η+ + η− )/(η+ η− )0(36)

In Eq. (36), we note that the first term diverges at k = 0,
because J1 (z) ≈ z when z
1 and J1 (kR)2 /k 3 ≈ 1/k inthis limit. This shows that frictional coupling of leaflets is
not enough to provide for a finite diffusion coefficient, andcontribution from solvents is essential. We get an interesting
result if we allow η− → ∞, which implies that the lower leafletof the membrane behaves as a solid. In this limit, we obtain
kB T
D=
π R2

 ∞
dk
0

[J1 (kR)]2k1.
η+ k 2 + (37)Equation (37) does give a finite result for the diffusion
coefficient, which can be explained by the lower leaflet
providing required friction for diffusion, instead of the solvent.
In the above, we show that the frictional coupling between
two leaflets alone is not enough to resolve the Stokes paradox.
Only in the case that one of the leaflets of the membrane isimmobile, a finite diffusion coefficient can be obtained without

022713-9

KAZUHIKO SEKI, SAURABH MOGRE, AND SHIGEYUKI KOMURA
(b)

2

2

1

1

0

0

y

y

(a)

PHYSICAL REVIEW E 89, 022713 (2014)

−1

−1

−2

−2
−2

−1

0

1

2

−2

x

−1

0

1

2

x

FIG. 10. (Color online) Velocity field for the supported case calculated from Eq. (18) when the domain is in the (a) upper and (b) lower
leaflets. The center of the domain moves with the velocity U in the positive direction of the x-axis. Parameters used are h = 0.1, ρ = 1, λ = 1.

solvents. The particular model was already applied to study theinterleaflet friction from the diffusion in a leaflet coupled to an
immobile leaflet [47], and the diffusion in a polymer-supportedmonolayer where the polymer plays the roll of an immobile
leaflet [48]. However, if both leaflets are viscous, the model
should include the coupling to solvents to have finite diffusion
coefficients.
ACKNOWLEDGMENTS

K.S. and S.K. are supported by Grant-in-Aid for Scientific
Research (Grant No. 24540439) from the MEXT of Japan. This
research was conducted under Agreement for Collaboration
between IIT Bombay and Nanosystem Research Institute,
AIST. S.K. would like to acknowledge support from theGrant-in-Aid for Scientific Research on Innovative Areas“Fluctuation & Structure” (Grant No. 25103010) from the
MEXT of Japan, and the JSPS Core-to-Core Program “International research network for nonequilibrium dynamics of soft
matter”.

space and is given by [24] (−α) 2
k + ηf(−α) k coth(kH (−α) ) δαβ + η,M (k) =  

(α)
(α) 2
(α) ) +  − 2α η k + ηf k coth(kH
αβ(A2)
where (−α) represents the leaflet opposite to α.
By noticing that the velocity at (R,0) is equal to U, thefriction coefficient, ζ , is obtained from the linear relation,
U = −(1/ζ )F. Using the Einstein relation D = kB T /ζ [50],
we obtain the diffusion coefficient as in Eq. (3).
The velocity fields can be drawn by using Eq. (A1). As an
example, we draw the velocity fields for the supported case ofSec. IV. The velocity field in the lower leaflet is shown when aliquid domain is located in the upper (lower) leaflet in Fig. 10.
By comparing the velocity fields, we could see that a part ofthe velocity field in the lower leaflet flows opposite to the flow
velocity of the liquid domain compared to that in the upper
leaflet.

APPENDIX A: VELOCITY FIELD

Denoting the total force on the circular domain in the steadystate by F ± , the velocity field is expressed as [49]v α (r) = − ∞dk M αα (k)

0

+

 ∞
0

dk M αα (k)J1 (kR)J1 (kr)
êx
π kRr

APPENDIX B: PARTIAL FRACTION EXPANSION TOSOLVE REQUIRED INTEGRALS

We now derive the result in Eq. (15) from Eqs. (5) and (6)using the method of partial fractions. A general relation used
to obtain partial fraction expansion for g(x)/f (x) is [51]
 g(xi )
1
g(x)
=
,

f (x) all roots f (xi ) (x − xi )J1 (kR)J2 (kr) sin(θ )êθ F α ,
πR

(B1)

(A1)
where êx is the unit vector in the direction of the x axis, êθ is theunit vector tangential to the circle pointed in the direction ofrotation, indices α, β denote the leaflets. In the above equation,
M αβ (k) represents a component of a mobility tensor in Fourier

where xi are the roots of the denominator, f (x), and the sum iscarried out over all roots. We find the partial fraction expansion
for the D1 (λ) term in Eq. (5). The other term can be obtained byputting λ = 0 and following the same procedure. The fraction
under consideration is 1/[κ 2 + κ coth(κh) + 2λ]. Equating the022713-10

DIFFUSION COEFFICIENTS IN LEAFLETS OF BILAYER . . .

PHYSICAL REVIEW E 89, 022713 (2014)

denominator to zero, we obtain,
κs2 + κs coth(κs h) + 2λ = 0,

(B2)

where s is the index of the root. It is convenient to introduce a transformation κs = iωs when an infinite number of roots exists.
This gives,
ωs2 = ωs cot(ωs h) + 2λ.

(B3)

We substitute cot(ωs h) = (ωs2 − 2λ)/ωs in Eq. (B1) to obtain,∞1
ωs / i
1
,
=
2
4
2
2
κ + κ coth(κh) + 2λ s=−∞ hωs + (1 + h − 4λh)ωs + 2λ + 4hλ κ − iωs(B4)where we can assume ω−s = −ωs without loss of generality. Since the expression within the parenthesis is even, we obtain

∞

2ωs2
1
1
=
.
(B5)
4 + (1 + h − 4λh)ω2 + 2λ + 4hλ2
2 + ω2κ 2 + κ coth(κh) + 2λ
hω
κ
s
s
s
s=1
If the first term in the expansion dominates, Eq. (6) becomesD1 (λ) =

2
kB T
4
2
π ηρ hω1 + (1 + h − 4λh)ω12 + 2λ + 4hλ2

1
− I1 (ω1 ρ)K1 (ω1 ρ) .
2

(B6)

The inflection region for the diffusion coefficient as afunction of λ is very different when ρ
1 in Fig. 3. The
inflection originates from the variation of the Bessel functions.
When ω1 ρ
1, the term including Bessel functions canbe approximated by a logarithmic function of the variable
by using the expansion, I1 (z)K1 (z)/z2 ∼ (1/4) ln z for z 
1. The logarithmic dependence is weak and the condition,
ω1 ρ
1, should be satisfied to obtain sharp variable dependence from the Bessel functions. By substituting ω1 ρ ∼ 1 into
the characteristic equation, Eq. (17), we obtain λ∗1 = [1/ρ −cot(h/ρ)]/(2ρ). At the inflection point, the term including
Bessel functions should depend sharply on the variables andthe variable itself should depend sharply on the value of λ.
We now discuss the λ dependence of the variable ω1 ρ. In thelimit√of √h → 0, we have cot(h/ρ) ∼ ρ/ h and obtain, ω1 ρ ∼(ρ/ h) 1 + 2λh. Therefore, the variable, ω1 ρ, depends
largely on λ if λ value exceeds λ∗2 = 1/(2h). By comparison,

we note that λ∗2 can be used even when h is as large as 10.
In summary, the larger of λ∗1 or λ∗2 characterizes the inflection
point. When ρ is small, the inflection point is characterized byλ∗1 . By increasing ρ, λ∗1 decreases and the inflection point ischaracterized by λ∗2 . The switching between λ∗1 and λ∗2 takes
place at ρ satisfying λ∗1 = λ∗2 . As seen in Fig. 3, if the inflection
point is characterized by λ∗1 , the inflection point can become
order of magnitude larger than λ∗2 by decreasing ρ.
The value characterizing the inflection point, √λ∗2 , can beobtained by equating the interleaflet sliding√length, η/b, with
the ES hydrodynamic screening length, h/ν. The velocity
field in the membrane perturbed by a diffusing particle is
2D-like within the ES length and becomes 3D-like if the
distance from the diffusing particle exceeds the ES length.
If the interleaflet sliding length is smaller than the ES length,
the 2D-like feature of the velocity field is not affected byinterleaflet sliding. On the other hand, if the interleaflet sliding
length is larger than the ES length, the 2D-like feature of thevelocity field is largely influenced by the interleaflet sliding
and the diffusion coefficient is increased by the sliding effect.

[1] P. G. Saffman and M. Delbrück, Proc. Natl. Acad. Sci. USA 72,
3111 (1975).
[2] P. G. Saffman, J. Fluid Mech. 73, 593 (1976).
[3] N. Oppenheimer and H. Diamant, Phys. Rev. E 82, 041912(2010).
[4] B. D. Hughes, B. A. Pailthorpe, and L. R. White, J. Fluid Mech.
110, 349 (1981).
[5] Y. Gambin, R. Lopez-Esparza, M. Reffay, E. Sierecki, N. S. Gov,
M. Genest, R. S. Hodges, and W. Urbach, Proc. Natl. Acad. Sci.
USA 103, 2098 (2006).
[6] P. Cicuta, S. L. Keller, and S. L. Veatch, J. Phys. Chem. B. 111,
3328 (2007).
[7] E. P. Petrov and P. Schwille, Biophys J. 94, L41 (2008).

[8] S. Ramadurai, A. Holt, V. Krasnikov, G. van den Bogaart, J. A.
Killian, and B. Poolman, J. Am. Chem. Soc. 131, 12650 (2009).
[9] E. Evans and E. Sackmann, J. Fluid. Mech. 194, 553 (1988).
[10] H. A. Stone and A. Ajdari, J. Fluid Mech. 369, 151 (1998).
[11] F. Harb, J. Sarkis, N. Ferte, and B. Tinland, Eur. Phys. J. E 35,
118 (2012).
[12] K. Simons and E. Ikonen, Nature (London) 387, 569 (1997).
[13] D. Lingwood and K. Simons, Science 327, 46 (2010).
[14] D. W. Allender and M. Schick, Biophys. J. 91, 2928 (2006).
[15] M. D. Collins and S. L. Keller, Proc. Natl. Acad. Sci. USA 105,124 (2008).
[16] R. De Koker, Ph.D. thesis, Stanford University, 1996
(unpublished).

APPENDIX C: INFLECTION POINT OF THE DIFFUSION
COEFFICIENT AS A FUNCTION OF λ022713-11

KAZUHIKO SEKI, SAURABH MOGRE, AND SHIGEYUKI KOMURA[17] S. Ramachandran, S. Komura, M. Imai, and K. Seki, Eur. Phys.
J. E 31, 303 (2010).
[18] K. Seki, S. Ramachandran, and S. Komura, Phys. Rev. E 84,
021905 (2011).
[19] Y. Fujitani, J. Phys. Soc. Jpn. 81, 084601 (2012).
[20] Y. Fujitani, J. Phys. Soc. Jpn. 82, 084403 (2013).
[21] U. Seifert and S. Langer, Europhys. Lett. 23, 71 (1993).
[22] H. A. Stone and H. M. McConnell, Proc. R. Soc. London A 448,
97 (1995).
[23] A. J. Levine and F. C. MacKintosh, Phys. Rev. E 66, 061606(2002).
[24] B. A. Camley and F. L. H. Brown, Soft Matter 9, 4767 (2013).
[25] T. Han and M. Haataja, Soft Matter 9, 2120 (2013).
[26] E. Evans and A. Yeung, Chem. Phys. Lipids 73, 39 (1994).
[27] R. M. Raphael and R. E. Waugh, Biophys. J. 71, 1374 (1996).
[28] T. Pott and P. Méléard, Europhys. Lett. 59, 87 (2002).
[29] A. F. Bitbol, J. B. Fournier, M. I. Angelova, and N. Puff, J. Phys.:
Condens. Matter 23, 284102 (2011).
[30] L. Zhang and S. Granick, J. Chem. Phys. 123, 211104 (2005).
[31] C. Scomparin, S. Lecuyer, M. Ferreira, T. Charitat, andB. Tinland, Eur. Phys. J. E 28, 211 (2009).
[32] M. Hetzer, S. Heinz, S. Grage, and T. M. Bayerl, Langmuir 14,
982 (1998).
[33] R. Macháň and M. Hof, Biochim. Biophys. Acta 1798, 1377
(2010).
[34] S. A. Shkulipa, W. K. den Otter, and W. J. Briels, Phys. Rev.
Lett. 96, 178302 (2006).
[35] W. K. den Otter and S. A. Shkulipa, Biophys. J. 93, 423 (2007).
[36] C. Xing and R. Faller, J. Phys. Chem. B 112, 7086 (2008).

PHYSICAL REVIEW E 89, 022713 (2014)[37] K. Inaura and Y. Fujitani, J. Phys. Soc. Jpn. 77, 114603 (2008).
[38] D. K. Lubensky and R. E. Goldstein, Phys. Fluids 8, 843(1996).
[39] Th. M. Fischer, J. Fluid Mech. 498, 123 (2004).
[40] S. Ramachandran, S. Komura, K. Seki, and G. Gompper, Eur.
Phys. J. E 34, 46 (2011).
[41] M. Abramowitz and I. A. Stegun, Handbook of Mathematical
Functions (Dover, New York, 1972).
[42] Wolfram Research Inc., MATHEMATICA 9 (Wolfram Research,
Champaign, 1988).
[43] In Eq. (39) of Ref. [18], the term proportional to 1/R 4 should be
eliminated. This is due to an error in MATHEMATICA 7. The otherparts including Eqs. (40)–(44) are not affected by the correction.
[44] J. B. Fournier, N. Khalifat, N. Puff, and M. I. Angelova, Phys.
Rev. Lett. 102, 018102 (2009).
[45] M. Tanaka and E. Sackmann, Nature (London) 437, 656 (2005).
[46] H.-Y. Zhang and R. J. Hill, J. R. Soc. Interface 8, 127 (2011).
[47] R. Merkel, E. Sackmann, and E. Evans, J. Phys. France 50, 1535
(1989).
[48] M. Tanaka, J. Hermann, I. Haase, M. Fischer, and S. G. Boxer,
Langmuir 23, 5638 (2007).
[49] This result is the correction of a similar result from a previous
publication [18]. The velocity at (R,0) was correct and the otherparts are not affected by the correction in Ref. [18].
[50] R. Kubo, M. Toda, and N. Hashitsume, Statistical Physics II:
Nonequilibrium Statistical Mechanics (Springer-Verlag, Berlin,1992).
[51] G. B. Arfken and H. J. Weber, Mathematical Methods for
Physicists, 5th ed. (Academic, San Diego, 2001).

022713-12