PHYSICAL REVIEW E 95, 052407 (2017)

Lateral diffusion induced by active proteins in a biomembrane
Yuto Hosaka, Kento Yasuda, Ryuichi Okamoto, and Shigeyuki Komura*
Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 14 March 2017; published 16 May 2017)
We discuss the hydrodynamic collective effects due to active protein molecules that are immersed in lipid
bilayer membranes and modeled as stochastic force dipoles. We specifically take into account the presence ofthe bulk solvent that surrounds the two-dimensional fluid membrane. Two membrane geometries are considered:
the free membrane case and the confined membrane case. Using the generalized membrane mobility tensors, weestimate the active diffusion coefficient and the drift velocity as a function of the size of a diffusing object. The
hydrodynamic screening lengths distinguish the two asymptotic regimes of these quantities. Furthermore, thecompetition between the thermal and nonthermal contributions in the total diffusion coefficient is characterizedby two length scales corresponding to the two membrane geometries. These characteristic lengths describe thecrossover between different asymptotic behaviors when they are larger than the hydrodynamic screening lengths.
DOI: 10.1103/PhysRevE.95.052407
I. INTRODUCTION

Biomembranes consisting of lipid bilayers can be regarded
as thin two-dimensional (2D) fluids, and membrane protein
molecules as well as lipid molecules are allowed to move laterally [1,2]. These membrane inclusions are subject to the thermal motion of lipid molecules, leading to random positional
fluctuations. Such Brownian motion plays important roles invarious life processes, such as the transportation of materials or
the reaction between chemical species [3]. To describe lateral
diffusion of membrane proteins, a drag coefficient of a cylindrical disk moving in a 2D fluid sheet has been studied theoretically in various membrane environments [4–10]. The obtaineddrag coefficient was used to estimate the diffusion coefficientsof membrane proteins through Einstein’s relation under theassumption that the system is in thermal equilibrium [11].
In recent experiments, however, it has been shown that
motions of particles inside cells are dominantly driven byrandom nonthermal forces rather than thermal fluctuations[12,13]. In these experimental works, they found that nonthermal forces in biological cells are generated by active
proteins undergoing conformational changes with a supplyof adenosine triphosphate (ATP). These active fluctuations
lead to enhanced diffusion of molecules in the cytoplasm[14,15]. Biomembranes also contain various active proteinsthat act, for example, as ion pumps by changing their shapes
to exert forces to the adjacent membrane and solvent [2].
Lipid bilayers containing such active proteins have beencalled “active membranes,” and their out-of-plane fluctuations(deformations) have already been investigated both experimentally and theoretically [16–18]. However, lateral motionsof inclusions in membranes that are induced by active proteins
have not yet been considered. Since such active forces giverise to enhanced diffusion, one needs to take into account bothactive nonthermal fluctuations as well as passive thermal ones
to calculate diffusion in membranes.
Recently, Mikhailov and Kapral discussed enhanced diffusion due to nonthermal fluctuating hydrodynamic flows, which
are induced by oscillating active force dipoles [see Fig. 1(a)]*komura@tmu.ac.jp

2470-0045/2017/95(5)/052407(10)[19,20]. They calculated the active diffusion coefficient ofa passive particle immersed either in a three-dimensional(3D) cytoplasm or in a 2D membrane, and they showed thatit exhibits a logarithmic size dependence for the 2D case.
Moreover, a chemotaxis-like drift of a passive particle waspredicted when gradients of active proteins or ATP are present[19]. Later, Koyano et al. showed that lipid membrane rafts, inwhich active proteins are concentrated, can induce a directed
drift velocity near the interface of a domain [21]. In theseworks, they considered membranes that are smaller in size thanthe hydrodynamic screening length. Huang et al. performedcoarse-grained simulations of active protein inclusions inlipid bilayers [22,23]. In Ref. [23], they showed that active
proteins undergoing conformational motions not only affectthe membrane shape but also laterally stir the lipid bilayer
so that lipid flows are induced. Importantly, the flow pattern
induced by an immobilized protein resembles the 2D fluid
velocity fields that are created by a force dipole.
Following Refs. [19,20], we assume that an active protein
behaves like an oscillating force dipole, which acts on thesurroundings to generate hydrodynamic flows that can induce
motions of passive particles in the fluid. In this paper, we investigate active diffusion and drift velocity of a particle in “free”
and “confined” membranes, which are completely flat andinfinitely large. In the free membrane case, a thin 2D fluid sheet
is embedded in a 3D solvent having typically a lower viscosity
than that of the membrane. In the confined case, which mimics
a supported membrane [24], a membrane is sandwiched by two
rigid walls separated by a finite but small distance from it. For
both the free and confined membrane cases, we employ generalmobility tensors that take into account the hydrodynamic
effects mediated by the surrounding 3D solvent [25–28]. Using
the general mobility tensors, we numerically calculate theactive diffusion coefficient and the drift velocity as a functionof the diffusing particle size for the entire length scales.
Furthermore, several asymptotic expressions are also derived
in order to compare with numerical estimates and thermal
contributions. Importantly, our result leads to characteristic
length scales describing a crossover from nonthermal tothermal diffusive behaviors for large scales.
In the next section, we present the expressions for the active
diffusion coefficient and the drift velocity in 2D membranes

052407-1©2017 American Physical Society

PHYSICAL REVIEW E 95, 052407 (2017)

HOSAKA, YASUDA, OKAMOTO, AND KOMURA

homogeneously and isotropically distributed in the membrane,
and they exert only in-plane lateral forces. The total diffusion
coefficient is given by D = DT + DA , where DT is the thermal
contribution and is determined by Einstein’s relation (which
will be discussed in Sec. IV), and DA is the active nonthermal
contribution given by [19]!
∂Gαβ (r) ∂Gαβ ′ (r)
Sc0
,
(1)
DA =
d 2 r "ββ ′ γ γ ′
2
∂rγ
∂rγ ′
where r = (x,y) denotes a 2D vector and we have introduced
a notation"ββ ′ γ γ ′ = 18 (δββ ′ δγ γ ′ + δβγ δβ ′ γ ′ + δβγ ′ δβ ′ γ ).

FIG. 1. (a) The conformational change of an oscillating force
dipole representing an active protein. Within a turnover cycle of theforce dipole separated by a distance x(t), it exerts two oppositely
directed forces ±F(t) at time t. The integral intensity of a force
dipole is S (see the text). (b) Schematic picture showing a flatand infinitely large membrane of 2D viscosity ηm that is locatedat z = 0. The membrane is surrounded by a bulk solvent of 3D
viscosity ηs , and the two flat walls are located at z = ±h. The solvent
velocity is assumed to vanish at the surfaces of these walls. The“free membrane” and the “confined membrane” cases correspond to
the limits of h → ∞ and h → 0, respectively. The yellow passive
particle undergoes Brownian motion due to thermal and nonthermal
fluctuations. The latter contribution is induced by active force dipolesthat are distributed homogeneously in the membrane with a 2D
concentration c0 .

[19]. We also review the general mobility tensors for the freeand confined membrane cases [25–28]. Using these expressions, we calculate in Sec. III the active diffusion coefficientfor the two geometries. In Sec. IV, we compare the thermal
diffusion coefficient with the obtained nonthermal diffusion
coefficient, and we discuss the characteristic crossover lengths.
In Sec. V, we obtain the drift velocities as a function of theparticle size. A summary of our work and some numericalestimates for the obtained quantities are given in Sec. VI.
II. ACTIVE TRANSPORT AND MOBILITY TENSORSIN MEMBRANES
A. Active diffusion coefficient

Active proteins in a 2D biological membrane, modeled as
oscillating force dipoles, produce nonequilibrium fluctuations
and cause an enhancement of the lateral diffusion of a passive
particle. We assume that the spatially fixed force dipoles are(2)

Throughout this paper, the summation over repeated greek
indices is assumed. In Eq. (1), S is the integral intensity ofa force dipole, c0 is the constant 2D concentration of active
proteins, and Gαβ (r) is the membrane mobility tensor, which
will be discussed later separately.
Within a fluctuating “dimer model” as presented inFig. 1(a), the magnitude of a force dipole is given by m(t) =x(t)F (t), where x(t) is the distance between the two spheres
and F (t) is the magnitude of the oppositely directed forces.
The statistical average of the dipole magnitude vanishes, i.e.,⟨m(t)⟩ = 0, whereas" ∞the integral intensity S of a force dipole
is given by S = 0 dt ⟨m(t)m(0)⟩ [19]. Since we assume
that active proteins are homogeneously distributed in themembrane as shown in Fig. 1(b), it is sufficient to consider
only the isotropic diffusion as given by Eq. (1).
In deriving Eq. (1), the size of a dipole is assumed to bemuch smaller than the distance between the passive particle
and active force dipoles [19]. At large distances, almost any
object that changes its shape would create a flow field that
corresponds to some force dipole. It should be noted, however,that the above expression is not accurate when the distance
between them becomes smaller. As for the mobility tensor
in 3D fluids, it is known that the Rotne-Prager mobility
tensor takes into account higher-order corrections to the Oseen
mobility tensor and gives a more accurate approximation atshort distances [20]. Such a better approximation has not been
worked out so far for 2D fluid membranes, and we shall only
consider the lowest-order contribution (see later calculations).
In the above, we have also assumed that force dipoles are
spatially fixed in the membrane. Since no forces are applied
to fix the dipoles, such an approximation is justified when
the dynamics of force dipoles is much slower than that of thepassive particle.
B. Drift velocity

Although we have assumed above that c0 is constant, active
proteins are often distributed inhomogeneously in the membrane due to heterogeneous structures such as sphingolipidenriched domains [29,30]. According to the “lipid raft”
hypothesis, these domains act as platforms for membrane
signaling and trafficking [31]. Hence it is also important toconsider the effects of nonuniform spatial distribution of active
proteins and to see how it affects the lateral dynamics inmembranes.

052407-2

PHYSICAL REVIEW E 95, 052407 (2017)

LATERAL DIFFUSION INDUCED BY ACTIVE PROTEINS . . .

When a spatial concentration gradient ∇c of active protein
is present, it gives rise to an unbalanced induced force between
two points in the membrane. Hence passive particles are
subjected to a drift toward either lower or higher concentrationof active proteins, and a chemotaxis-like drift can occur.
When the absolute value of the concentration gradient |∇c| isassumed to be constant, the induced drift velocity of a passive
particle in the direction ∇c is given by [19]!
∂ 2 Gαβ (r) ∂Gδβ ′ (r)V = −S|∇c| d 2 r "ββ ′ γ γ ′ n̂α(r · n̂). (3)∂rγ ∂rδ∂rγ ′Here, the unit vector n̂ = ∇c/|∇c| denotes the direction of
the concentration gradient of active proteins. We shall employ
the above expression to obtain the lateral drift velocity in amembrane by using the membrane mobility tensor as discussed
below.
C. Membrane mobility tensors

Since we discuss active diffusion in an infinitely large
flat membrane, we use the 2D membrane mobility tensor,
which also takes into account the hydrodynamic effects ofthe surrounding 3D solvent. We consider a general situation
as depicted in Fig. 1(b), where a fluid membrane of 2D shear
viscosity ηm is surrounded by a solvent of 3D shear viscosity
ηs . Furthermore, we consider the case in which there are twowalls located symmetrically at an arbitrary distance h from theflat membrane [25–28].
We denote the in-plane velocity vector of the fluid membrane by v(r) and the lateral pressure by p(r). Assuming thatthe incompressibility condition holds for the fluid membrane,
we write its hydrodynamic equations as∇ · v = 0,(4)

ηm ∇ 2 v − ∇p + fs + F = 0.

(5)

# · v± = 0,∇

(6)

The second equation is the 2D Stokes equation, where fs is theforce exerted on the membrane by the surrounding solvent,and F is any external force acting on the membrane. If wedenote the upper and lower solvents with the superscripts ±,the two solvent velocities v± (r,z) and pressures p± (r,z) obey
the following hydrodynamic equations, respectively:# ± = 0,#2 v± − ∇pηs ∇vα (r) = vα± (r,0),Gαβ (k) =(8)

where α = x,y. Furthermore, the solvent velocity vanishes at
the walls located at z = ±h, i.e., vα± (r, ± h) = 0.
By solving the above coupled hydrodynamic equations inFourier space with k = (kx ,ky ) being the 2D wave vector,
the 2D mobility tensor Gαβ (k) defined through vα (k) =ηm$δαβ − k̂α k̂β%,2k + νk coth(kh)(9)

where k = |k| and k̂α = kα /k, and the ratio of the two viscosities ν −1 = ηm /(2ηs ) defines the Saffman-Delbrück (SD)hydrodynamic screening length [4,5]. Notice that ηm and ηshave different dimensions, and ν −1 has a dimension of length.
To perform analytical calculations, the two limiting casesof Eq. (9) are considered, i.e., the “free membrane” case andthe “confined membrane” case corresponding to the limitsof h → ∞ and h → 0, respectively [26–28]. For the freemembrane case, we take the limit kh ≫ 1 in Eq. (9) and obtainthe following asymptotic expression:
GFαβ (k) =δαβ − k̂α k̂β.
ηm (k 2 + νk)(10)δαβ − k̂α k̂β,ηm (k 2 + κ 2 )(11)

Hereafter, we shall denote the quantities for the free membranecase with the superscript “F.” For the confined membrane case,on the other hand, we take the opposite limit kh ≪ 1 and obtain
GCαβ (k) =where κ −1 = (h/ν)1/2 is the Evans-Sackmann (ES) screening
length [7], and we use the superscript “C” for the quantities
related to the confined membrane case. We note that the ES
screening length κ −1 is the geometric mean of ν −1 and h sothat we typically have κ −1 < ν −1 .
Taking the inverse Fourier transform of Eqs. (10) and (11),
we obtain the mobility tensors in real space for the two limitingcases as [26–28]&
1
2F
H0 (νr) − Y0 (νr) +Gαβ (r) =4ηm
π ν2r 2'H1 (νr) Y1 (νr)
+
δαβ
−
νr
νr
&
4
1
2H1 (νr)− 2 2+
+
4ηm
πν r
νr
'
2Y1 (νr)− H0 (νr) + Y0 (νr) r̂α r̂β
(12)
−
νr
and

(7)

# stands for the 3D differential operator.
where ∇
We assume that the surrounding solvent cannot permeate
the membrane, and we impose the no-slip boundary condition
between the membrane and the surrounding solvent at z = 0[4,5,25–28]. Hence we require the conditions
vz± (r,0) = 0,Gαβ (k)Fβ (k) can be obtained as [25–28]

GCαβ (r) =&'1
K1 (κr)1
K0 (κr) +− 2 2 δαβ
2π ηm
κr
κ r
&
'
2
1
2K1 (κr)−K0 (κr) −+ 2 2 r̂α r̂β ,
+
2π ηm
κr
κ r
(13)

respectively, where we have used the notations r = |r| andr̂α = rα /r. In the above, Hn (z) are the Struve functions,
Yn (z) are the Bessel functions of the second kind, and Kn (z)are the modified Bessel functions of the second kind. The
physical meaning of the above expressions was also discussed
in Refs. [32–34]. We note that if there is only one wall insteadof two, the definition of the ES length needs to be modified
as κ −1 → (2h/ν)1/2 [34]. In the next sections, we shall use

052407-3

PHYSICAL REVIEW E 95, 052407 (2017)

HOSAKA, YASUDA, OKAMOTO, AND KOMURATo understand the above behaviors, we next discuss theasymptotic behaviors of DAF for both small and large δ values.
Expanding the mobility tensor in Eq. (12) for νr ≪ 1 andνr ≫ 1, we have [34]
)
(
1
2
F
δαβ + ẑα ẑβgαβ(z) ≈ ln − γ −
(15)
z
2

0

10

-3

10

and

-6

F
gαβ
(z) ≈

10

-2

0

10

(16)

respectively, where γ = 0.5722 . . . is Euler’s constant. By
substituting Eqs. (15) and (16) into Eq. (14), we can analytically obtain the asymptotic forms of the active diffusion
coefficient as a function of δ = νℓc .
As obtained in Ref. [19], we find for δ ≪ 1

-9

10

10

2
ẑα ẑβ ,z

2

10FIG. 2. The plot of the scaled active diffusion coefficient DA as afunction of the scaled cutoff length δ = νℓc and ϵ = κℓc for the freemembrane case [solid line, see Eq. (14)] and the confined membranecase [dashed line, see Eq. (19)], respectively. Here DA is scaled by
Sc0 /(256π ηm2 ). The numbers in this plot indicate the slope of thecurves and represent the powers of the algebraic dependencies.

Eqs. (12) and (13) to calculate the active diffusion coefficientsand the drift velocity.
III. ACTIVE DIFFUSION COEFFICIENT

DAF ≈

F
where δ = νℓc is the dimensionless cutoff, and gαβ(z) =F
4π ηm Gαβ is the corresponding dimensionless mobility tensor[see Eq. (12)]. We have first evaluated the above integral
numerically. In Fig. 2, we plot the obtained DAF as a functionof δ = νℓc by the solid line. We see that the active diffusion
coefficient depends only weakly on the particle size at small
scales, whereas it shows a stronger size dependence describedby a power-law behavior at large scales. The crossover between
these two behaviors is set by the condition δ ≈ 1.

(17)

where a large cutoff length L is introduced because theintegral in Eq. (14) also diverges logarithmically at large
distances. To match with the numerical estimation, we obtainL ≈ 0.682ν −1 . The above logarithmic dependence on ℓcmeans that DAF depends only weakly on the particle size. We
also note that the above expression contains only the membrane
viscosity ηm , and does not depend on the solvent viscosity ηs .
This is because the hydrodynamics at small scales is primarily
dominated by the 2D membrane property.
In the opposite limit of δ ≫ 1, on the other hand, we show
in Appendix A that the active diffusion coefficient becomes

A. Free membranes

We first calculate the active diffusion coefficient for the freemembrane case by substituting Eq. (12) into Eq. (1). Since theintegrand in Eq. (1) diverges logarithmically at short distances,
we need to introduce a small cutoff length ℓc . Physically, ℓc is
given by the sum of the size of a passive particle (undergoinglateral Brownian motion) and that of a force dipole [20]. In thefollowing, we generally assume that force dipoles are smaller
than the diffusing object, whose size is represented by ℓc .
This is further justified when we consider lateral diffusion ofa passive object that is larger than the SD or ES screening
lengths.
Introducing a dimensionless vector z = νr scaled by the
SD length, we can write the active diffusion coefficient for thefree membrane case as
! ∞
F
F
∂gαβ
(z) ∂gαβ
′ (z)
Sc0
F
2
′γ γ ′
d
z
"
,
(14)
DA =
ββ
2
32π 2 ηm
∂z
∂z
γ
γ′
δ

Sc0
L
ln ,
2
32π ηmℓc

DAF ≈5Sc0 1,256π ηs2 ℓ2c(18)which is an important result of this paper. This asymptotic
expression decays as 1/ℓ2c and depends now only on ηs ,indicating that the membrane lateral dynamics is governedby the surrounding 3D fluid at large scales. From the obtained
asymptotic expressions in Eqs. (17) and (18), the behavior of
DAF in Fig. 2 is explained as a crossover from a logarithmic
dependence to an algebraic dependence with a power of −2.
B. Confined membranes

Next we consider the confined membrane case. With theuse of Eq. (13), the active diffusion coefficient can be written
as
! ∞
C
C
∂gαβ
(w) ∂gαβ
′ (w)
Sc0
2
′
′
DAC =
d
w
"
, (19)ββ γ γ22
32π ηm ϵ∂wγ∂wγ ′where w = κr is a different dimensionless variable, ϵ = κℓc
Cis a differently scaled cutoff, and gαβ(w) = 4π ηm GCαβ is thecorresponding dimensionless mobility tensor [see Eq. (13)].
Performing the numerical integration of Eq. (19), we plot in
Fig. 2 the active diffusion coefficient DAC as a function ofϵ = κℓc by the dashed line. For small ϵ values, the behavior
of DAC is similar to that of DAF , while DAC decays much fasterthan DAF for large ϵ values.

052407-4

PHYSICAL REVIEW E 95, 052407 (2017)

LATERAL DIFFUSION INDUCED BY ACTIVE PROTEINS . . .

To discuss these size dependencies, we use the asymptotic
expressions of Eq. (13) for κr ≪ 1 and κr ≫ 1 given by [34]
)
(
1
2
C
δαβ + ŵα ŵβ(w) ≈ ln − γ −
(20)
gαβ
w
2

0

10

and

-3

C
gαβ
(w) ≈ −

2(δαβ − 2ŵα ŵβ ),
w2

10
(21)

-6

10

respectively. Note that Eq. (20) is identical to Eq. (15) when
w is replaced by z. Hence, in the limit of ϵ ≪ 1, the active
diffusion coefficient for the confined membrane case should
be identical to Eq. (17) and is given by [19]
DAC ≈

Sc0
L
ln .
2
32π ηm ℓc

-9

10

(22)
-2

10

−1

The large cutoff length should be taken here as L ≈ 1.12κ .
As mentioned before, the 2D hydrodynamic effect is more
important at small scales, and DAC is logarithmically dependent
on the particle size.
In the large size limit of ϵ ≫ 1, on the other hand, we also
show in Appendix A that DAC behaves asymptotically as
DAC ≈Sc0 h2,16π ηs2 ℓ4c(23)which is another important result. The obtained expression
decays as 1/ℓ4c , which is much stronger than Eq. (18) for thefree membrane case. According to Eqs. (22) and (23), thebehavior of DAC in Fig. 2 can be understood as a crossoverfrom a logarithmic dependence to an algebraic dependencewith a power of −4.
IV. TOTAL DIFFUSION COEFFICIENT

Having obtained the active diffusion coefficients for thefree and the confined membrane cases, we now discuss the totallateral diffusion coefficients in membranes by considering both
thermal and nonthermal contributions. Concerning the thermal
diffusion coefficient DTF for the free membrane case, we use anempirical expression obtained by Petrov and Schwille [35,36],
&
'
2
4δ
δ2 2
kB T
ln − γ +− ln
DTF (δ) =
4π ηm
δ
π
2
δ
&
'
−1
δ3 2
c1 δ b1
× 1 − ln +
,
(24)
π
δ
1 + c2 δ b2

where kB is the Boltzmann constant, T is the temperature,and the four numerical constants are chosen as c1 = 0.737 61,
b1 = 2.748 19, c2 = 0.521 19, and b2 = 0.514 65 [36]. For thefree membrane case, there is no exact analytical expression of
the thermal diffusion coefficient that covers the entire size
range, except for the case in which a 2D polymer chain is
confined in a fluid membrane [26]. Equation (24) is known torecover the correct asymptotic limits of the thermal diffusion
coefficients both for δ ≪ 1 [4,5] and δ ≫ 1 [6].
On the other hand, the thermal diffusion coefficient DTC
for the confined membrane case was explicitly calculated by
Evans et al. [7] and also by Ramachandran et al. [8–11]. In

0

10

2

10

FIG. 3. The plot of the scaled thermal diffusion coefficient DT asa function of the scaled cutoff length δ = νℓc and ϵ = κℓc for the freemembrane case [solid line, see Eq. (24)] and the confined membranecase [dashed line, see Eq. (25)], respectively. Here DT is scaled by
kB T /(4π ηm ). The numbers in this plot indicate the slope of the curves
and represent the powers of the algebraic dependencies.

this case, the resulting expression is given by&'−1
ϵK1 (ϵ)kB T ϵ 2+
DTC (ϵ) =.
4π ηm 4K0 (ϵ)(25)

In Fig. 3, we plot DTF as a function of the particle size δ by thesolid line, and DTC as a function of ϵ by the dashed line for thewhole size range. Their asymptotic behaviors are discussed
separately below.
When we consider the total diffusion coefficient D =DT + DA , we shall neglect the contribution from thermal
fluctuations of force dipoles. These fluctuations can arise whenforce dipoles contain structural internal degrees of freedom.
However, such a contribution to the diffusion coefficient issmall compared to DT because it should be proportional to the
product of kB T and the concentration of force dipoles c0 .
A. Free membranes

For the free membrane case, the total diffusion coefficientis given by D F = DTF + DAF , where the active nonthermal
contribution DAF was discussed in the previous section. Using
Eqs. (24) and (17) in the limit of δ ≪ 1, we asymptoticallyhave [4,5]
(
)
2
kB T
Sc0
L
ln
DF ≈
−γ +
ln ,
(26)
2
4π ηm
νℓc
32π ηmℓcwhere both contributions are proportional to ln(1/ℓc ).
For δ ≫ 1, on the other hand, we obtain from Eqs. (24) and(18) [6]052407-5

DF ≈kB T 1
5Sc0 1+.
16ηs ℓc256π ηs2 ℓ2c(27)

PHYSICAL REVIEW E 95, 052407 (2017)

HOSAKA, YASUDA, OKAMOTO, AND KOMURA

Since the ℓc -dependencies in Eq. (27) are different between
the thermal and nonthermal contributions, we can introduce a
new crossover length defined by
ℓ∗ =

5Sc0
.
16π kB T ηs

0

10

(28)
-3

10

This length scale characterizes a crossover from the 1/ℓ2c dependence to the 1/ℓc -dependence. When ℓc ≪ ℓ∗ (but stillν −1 ≪ ℓc ), the nonthermal contribution dominates over thethermal one, while in the opposite limit of ℓc ≫ ℓ∗ the thermal
contribution is of primary importance.

-5

10

B. Confined membranes-8

In the case of confined membranes, the total diffusion
coefficient now becomes D C = DTC + DAC . In the limit ofϵ ≪ 1, we have from Eqs. (25) and (22) [7,8]
(
)
2
kB T
Sc0
L
C
ln
D ≈
−γ +
ln ,
(29)
2
4π ηm
κℓc
32π ηm ℓcwhere both contributions exhibit a logarithmic dependence onℓc as in the free membrane case.
In the opposite limit of ϵ ≫ 1, we find from Eqs. (25) and(23) [7,8]DC ≈kB T hSc0 h2+.
2π ηs ℓ2c
16π ηs2 ℓ4c(30)

Similar to the free membrane case, we can consider anothercharacteristic length defined by()
Sc0 h 1/2
ℓ∗∗ =
.
(31)
8kB T ηs

This length scale characterizes a crossover from the 1/ℓ4c dependence to the 1/ℓ2c -dependence. We note that ℓ∗∗ isessentially the geometric mean of ℓ∗ and h. Numericalestimates of these two characteristic length scales will be
discussed in Sec. VI.
V. DRIFT VELOCITY
A. Free membranes

In this section, we calculate the drift velocity V of a passive
particle due to a concentration gradient of active force dipoles.
For the free membrane case, we substitute Eq. (12) into Eq. (3)and obtain
! ∞
S|∇c|
F
V =−
d 2 z "ββ ′ γ γ ′2
16π 2 ηm
δ
× n̂α

F
F
∂ 2 gαβ(z) ∂gδβ′ (z)∂zγ ∂zδ∂zγ ′(z · n̂),(32)

F
where δ = νℓc and gαβ(z) = 4π ηm GFαβ as before. Performing
the numerical integration of Eq. (32), we plot in Fig. 4 the drift
velocity V F as a function of δ by the solid line. Similar to theactive diffusion coefficient DAF , the drift velocity V F depends
weakly on the particle size at small scales, while it exhibits astronger size dependence at large scales. Such a crossover also
occurs around δ ≈ 1.
We next discuss the asymptotic behaviors of V F for small
and large δ values. With the use of Eqs. (15) and (16), we show

10

-2

0

10

2

10

10

FIG. 4. The plot of the scaled drift velocity V as a function of thescaled cutoff length δ = νℓc and ϵ = κℓc for the free membrane case[solid line, see Eq. (32)] and the confined membrane case [dashed line,
see Eq. (35)], respectively. Here V is scaled by S|∇c|/(128π ηm2 ). The
numbers in this plot indicate the slope of the curves and represent thepowers of the algebraic dependencies.

in Appendix B that the asymptotic behaviors of V for δ ≪ 1
and δ ≫ 1 are
VF ≈

S|∇c|
L
ln
2
32π ηm
ℓc

(33)

VF ≈

13S|∇c| 1,256π ηs2 ℓ2c

(34)

and

respectively, where we choose L ≈ 1.85ν −1 . Note that Eq. (33)
was previously derived in Ref. [19] for a 2D membrane, while
Eq. (34) is a new result. As we see in Eqs. (33) and (34), there
is a crossover from a logarithmic to an algebraic dependencewith a power of −2 when δ is increased. These behaviors
are consistent with the numerical plot in Fig. 4 for the freemembrane case.
B. Confined membranes

Finally, we calculate the drift velocity for the confined
membrane case. Substituting Eq. (13) into Eq. (3), we now
obtain
! ∞
S|∇c|
VC =−
d 2 w "ββ ′ γ γ ′2
16π 2 ηm
ϵ
× n̂α

C
C
∂ 2 gαβ(w) ∂gδβ′ (w)∂wγ ∂wδ

∂wγ ′(w · n̂),(35)

C
where ϵ = κℓc and gαβ(w) = 4π ηm GCαβ as before. In Fig. 4,
we present numerically calculated V C as a function of ϵ bythe dashed line. As ϵ is increased, we see a crossover from alogarithmic to an algebraic dependence, although V C decays
faster than V F at large scales.

052407-6

PHYSICAL REVIEW E 95, 052407 (2017)

LATERAL DIFFUSION INDUCED BY ACTIVE PROTEINS . . .

TABLE I. Summary of the asymptotic dependencies of the thermal diffusion coefficient DT , the active diffusion coefficient DA , and thedrift velocity V on the passive particle size ℓc . The numbers after the asymptotic expressions correspond to the equation numbers in this paper.
Cases

LimitsDT

DA

V

free membrane(hk ≫ 1)
confined membrane(hk ≪ 1)νℓc ≪ 1νℓc ≫ 1κℓc ≪ 1κℓc ≫ 1ln(1/ℓc ) (26)
1/ℓc (27)ln(1/ℓc ) (29)
1/ℓ2c (30)ln(1/ℓc ) (17)1/ℓ2c (18)ln(1/ℓc ) (22)
1/ℓ4c (23)ln(1/ℓc ) (33)
1/ℓ2c (34)ln(1/ℓc ) (36)1/ℓ4c (37)The asymptotic behaviors of V C for small and large ϵ values
can be discussed similarly. Using Eqs. (20) and (21), we obtainin Appendix B the asymptotic expressions of V C for ϵ ≪ 1and ϵ ≫ 1 as
VC ≈

S|∇c|
L
ln
2
32π ηm
ℓc

(36)

VC ≈

3S|∇c| h2,
16π ηs2 ℓ4c

(37)

and

respectively, and we choose L ≈ 3.05κ −1 to coincide with thenumerical integration. We note that Eqs. (33) and (36) are
identical and depend only on ηm for small sizes [19].
From Fig. 4 and Eqs. (33), (34), (36), and (37), we seethat the drift velocity V is always positive. This meansthat passive particles move toward higher concentrations ofactive proteins, and a chemotaxis-like drift takes place inthe presence of protein concentration gradients [19–21]. The
dominant viscosity dependence of V switches from ηm to ηs
as the particle size exceeds the corresponding hydrodynamic
screening length, namely ν −1 or κ −1 .
VI. SUMMARY AND DISCUSSION

In this paper, we have investigated lateral diffusion inducedby active force dipoles embedded in a biomembrane. In
particular, we have calculated the active diffusion coefficientand the drift velocity for the free and the confined membranecases by taking into account the hydrodynamic coupling
between the membrane and the surrounding bulk solvent. The
force dipole model in Refs. [19,20] and the general membrane
mobility tensors obtained in Refs. [25–28] have been employedin our work. When the size of a passive diffusing particle
is small, the active diffusion coefficients for the free andthe confined membranes represent the same logarithmic size
dependence, as shown in Eqs. (17) and (22), respectively [19].
In the opposite large size limit, we find algebraic dependencies
with powers −2 and −4 for the two cases, as given by Eqs. (18)and (23), respectively. These are the important outcomes of thispaper, and they are also summarized in Table I together withother asymptotic expressions.
In our work, we have assumed that the total diffusion
coefficient is provided by the sum of thermal and nonthermal
contributions. For small particle sizes, we have shown that boththe total D F and D C exhibit a logarithmic size dependence [19],whereas different contributions have different size dependencies for large particle sizes. From this result, we have obtained

two characteristic length scales that describe the crossoverfrom nonthermal to thermal behaviors when the particle sizeis larger than the hydrodynamic screening length. The drift
velocity in the presence of a concentration gradient of active
proteins exhibits the same size dependencies as the active
diffusion coefficient for the two membrane geometries.
Here we give some numerical estimates of the obtained
crossover length scales. Using typical values such as kB T ≈4 × 10−21 J, ηs ≈ 10−3 Pa s, h ≈ 10−9 m, S ≈ 10−42 J2 s, andc0 ≈ 1014 m−2 [19], we obtain ℓ∗ ≈ 2 × 10−6 m [see Eq. (28)]and ℓ∗∗ ≈ 6 × 10−8 m [see Eq. (31)]. On the other hand, the SD
and the ES screening lengths are typically ν −1 ≈ 5 × 10−7 m
and κ −1 ≈ 2 × 10−8 m, respectively [4,5,7,8]. Hence ℓ∗ andℓ∗∗ are typically larger than ν −1 and κ −1 , respectively.
Moreover, the values of S and c0 can vary significantly in
one membrane to another, as pointed out in Ref. [19]. For
example, when active proteins are confined in raft domains[29–31], the 2D concentration c0 can be much larger. When, forexample, c0 ≈ 1015 m−2 (while S is the same as above) [21],the crossover length can be estimated as ℓ∗ ≈ 2 × 10−5 m
and ℓ∗∗ ≈ 2 × 10−7 m. If ℓ∗ and ℓ∗∗ are much larger than
the screening lengths ν −1 and κ −1 , respectively, as in thiscase, the three different scaling regimes of the total diffusion
coefficient are expected as the particle size is increased, i.e.,ln(1/ℓc ) → 1/ℓ2c → 1/ℓc for the free membrane case, andln(1/ℓc ) → 1/ℓ4c → 1/ℓ2c for the confined membrane case.
Momentum in a membrane is conserved over distances
smaller than the hydrodynamic screening length (either ν −1or κ −1 ), whereas it leaks to the surrounding fluid beyond that
length scale [32–34]. Within a membrane, the velocity decays
as ln(1/r) at short distances, as shown in Eqs. (15) and (20),due to the momentum conservation in two dimensions. These2D behaviors also lead to the logarithmic dependence of theactive diffusion coefficients in Eqs. (17) and (22). For the freemembrane case, the velocity decays as 1/r at large scales
as shown in Eq. (16) due to the momentum conservation in
the 3D bulk. This behavior is reflected in the first term ofEq. (27) for the thermal diffusion coefficient [6]. As shown in
Eq. (21), however, the velocity decays as 1/r 2 at large scales
for the confined membrane case. This behavior essentiallyarises from the mass conservation in two dimensions while thetotal momentum is not conserved due to the presence of thewalls, which break the translational symmetry of the system[32–34]. The corresponding contribution is the first term ofEq. (30) for the thermal diffusion coefficient [7,8].
The active diffusion coefficient DAF obtained in Eq. (18) for
the free membrane case essentially reflects the hydrodynamics
of the surrounding bulk 3D solvent. Hence our result can becompared with that in Ref. [19] obtained for a purely 3D fluid

052407-7

PHYSICAL REVIEW E 95, 052407 (2017)

HOSAKA, YASUDA, OKAMOTO, AND KOMURA

system:APPENDIX A: DERIVATION OF EQS. (18) AND (23)

DA3D ≈

Sc03D 1,60π ηs2 ℓc(38)which decays as 1/ℓc and is different from Eq. (18). In fact,
such a difference arises from the different dimensions of thedipole concentrations, i.e., c0 is the 2D concentration in ourcase, while c03D is the 3D concentration in Ref. [19]. A similar
comparison can also be made for the drift velocity of free membranes in Eq. (34) and that in Ref. [19] for a 3D fluid system:
V 3D ≈

S|∇c3D | 1.
30π ηs2 ℓc(39)

The same reason holds for the different ℓc -dependence.
At this stage, we also comment that both the active
diffusion coefficient DA and the drift velocity V exhibit thesame ℓc -dependence. Although the integrands in Eqs. (1) and(3) look apparently different, their physical dimensions are
identical because the first derivative of the mobility tensor in
Eq. (1) corresponds to the product of the second derivative and(r · n̂) in Eq. (3). This is the simple reason that they exhibit
the same ℓc -dependence. One can also easily confirm that
V is positive when we make use of the membrane mobility
tensor, because the integrand in Eq. (3) is the product of
the first and the second derivatives of the mobility tensor,
which have opposite signs. This leads to V > 0 indicating achemotaxis-like drift, as mentioned before.
In this work, we have assumed that active proteins generate
forces only in the lateral directions. On the other hand, actual
active motors such as bacteriorhodopsin can also exert forces tothe surrounding solvent [16–18]. Although we did not take into
account such normal forces that induce membrane undulation,
consideration of normal forces as well as lateral ones will
provide us with a general understanding of active diffusion inbiomembranes [37].
We have also assumed that the force dipoles are fixed ina membrane and are distributed homogeneously. It would beinteresting to consider the case when active proteins can also
move laterally in the membrane and even interact with eachother through a nematic-like interaction [38]. The full equationof motion now involves potential-of-mean-force interactionsin the multiparticle diffusion equations that describe the
combined motions of the passive particle and active proteinsin the membrane. Although the dynamics of the active
protein concentration is essentially determined by a diffusion
equation, it is a complicated problem because not only thermal
diffusion but also active nonthermal diffusion should be taken
into account. Our work is a step toward such a full descriptionof very rich biomembrane dynamics.

Since Eqs. (17) and (22) have been obtained in Ref. [19],we show here the derivation of Eqs. (18) and (23). Substituting
Eq. (16) into Eq. (14), we get! ∞Sc0
d 2 z "ββ ′ γ γ ′DAF = 2 28π ηm δ(())
ẑα ẑβ ′∂∂ ẑα ẑβ
,
(A1)
×
∂zγ
z
∂zγ ′
z
where z = νr. Since()∂ ẑα ẑβ13= 3 (δαγ zβ + δβγ zα ) − 5 zα zβ zγ ,
∂zγ
z
z
z
the integrand in Eq. (A1) becomes(())
ẑα ẑβ ′∂∂ ẑα ẑβ
∂zγ
z
∂zγ ′
z

1
1
δβγ δβ ′ γ ′ + 6 [δγ γ ′ zβ zβ ′ − 2(δβγ zβ ′ zγ ′ + δβ ′ γ ′ zβ zγ )]4z
z
3+ 8 zβ zβ ′ zγ zγ ′ .
(A3)z
By operating "ββ ′ γ γ ′ , we have(())
ẑα ẑβ ′∂5∂ ẑα ẑβ"ββ ′ γ γ ′
= 4.
(A4)
∂zγ
z
∂zγ ′
z
8z
=

After the integration, we obtain Eq. (18).
Similarly, we substitute Eq. (21) into Eq. (19) and obtain
! ∞
Sc0
C
DA = 2 2
d 2 w "ββ ′ γ γ ′8π ηm ϵ()()δαβ − 2ŵα ŵβ∂δαβ ′ − 2ŵα ŵβ ′
∂
,
×
∂wγ
w2
∂wγ ′
w2
(A5)
where w = κr. Since()∂δαβ − 2ŵα ŵβ
∂wγ
w2
=−

8
2
(δαβ wγ + δβγ wα + δαγ wβ ) + 6 wα wβ wγ ,
4
w
w
(A6)

we obtain()()∂δαβ − 2ŵα ŵβ∂δαβ ′ − 2ŵα ŵβ ′
∂wγ
w2
∂wγ ′
w2
=

ACKNOWLEDGMENTS

We thank A. S. Mikhailov and T. Kato for useful discussions. S.K. and R.O. acknowledge support from a Grant-in-Aid
for Scientific Research on Innovative Areas “Fluctuationand Structure” (Grant No. 25103010) from the Ministry of
Education, Culture, Sports, Science, and Technology of Japan,
and a Grant-in-Aid for Scientific Research (C) (Grant No.
15K05250) from the Japan Society for the Promotion of
Science (JSPS).

(A2)

4
4
δβγ δβ ′ γ ′ + 8 [δββ ′ wγ wγ ′ + δβ ′ γ wβ wγ ′6ww
+ δβγ ′ wβ ′ wγ + δγ γ ′ wβ wβ ′ − 2(δβγ wβ ′ wγ ′+ δβ ′ γ ′ wβ wγ )].

(A7)

By operating "ββ ′ γ γ ′ , we have()()
δαβ − 2ŵα ŵβ∂δαβ ′ − 2ŵα ŵβ ′∂′′"ββ γ γ
∂wγ
w2
∂wγ ′
w2
4
.
w6
After the integration, we obtain Eq. (23).

052407-8=(A8)

PHYSICAL REVIEW E 95, 052407 (2017)

LATERAL DIFFUSION INDUCED BY ACTIVE PROTEINS . . .

APPENDIX B: DERIVATION OF EQS. (34) AND (37)In this appendix, we show the derivation of Eqs. (34) and (37). Substituting Eq. (16) into Eq. (32), we obtain
(
)
(
)
!
ẑα ẑβ∂ẑδ ẑβ ′∂2S|∇c| ∞ 2F(z · n̂).
d z "ββ ′ γ γ ′ n̂αV =−2
4π 2 ηm∂zγ ∂zδ
z
∂zγ ′
z
δ

In the above, the derivatives are()
ẑα ẑβ1∂2
3
15= 3 (δαδ δβγ + δαγ δβδ ) − 5 (δαδ zβ zγ + δβδ zα zγ + δαγ zβ zδ + δβγ zα zδ + δγ δ zα zβ ) + 7 zα zβ zγ zδ∂zγ ∂zδ
z
z
z
z

and

(B1)

(B2)

(
(
)
)
ẑα ẑβẑδ ẑβ ′∂1∂2= − 6 [2δβ ′ γ ′ (δαγ zβ + δβγ zα ) − (δαγ ′ δβγ + δαγ δβγ ′ )zβ ′ ]∂zγ ∂zδ
z
∂zγ ′
z
z

3
[(δαγ ′ zβ zγ + δβγ ′ zα zγ − δαγ zβ zγ ′ − δβγ zα zγ ′ + δγ γ ′ zα zβ )zβ ′ − 2δβ ′ γ ′ zα zβ zγ ]z8
3
− 10 zα zβ zβ ′ zγ zγ ′ .
(B3)
z

−

By operating "ββ ′ γ γ ′ , we have()()
ẑα ẑβ∂ẑδ ẑβ ′13zα∂2"ββ ′ γ γ ′=− 6 .
∂zγ ∂zδ
z
∂zγ ′
z
8z

After the integration, we obtain Eq. (34).
Next we substitute Eq. (21) into Eq. (35) and find
(
)
(
)
!
δαβ − 2ŵα ŵβ∂δδβ ′ − 2ŵδ ŵβ ′∂2S|∇c| ∞ 2
C
d w "ββ ′ γ γ ′ n̂αV =−(w · n̂).
2
4π 2 ηm∂wγ ∂wδ
w2
∂wγ ′
w2
ϵ

(B4)

(B5)

Here the derivatives are()δαβ − 2ŵα ŵβ28∂2= − 4 (δαβ δγ δ + δαγ δβδ + δαδ δβγ ) + 6 (δαβ wγ wδ + δβδ wα wγ + δαδ wβ wγ∂wγ ∂wδ
w2w
w+ δαγ wβ wδ + δβγ wα wδ + δγ δ wα wβ ) −48
wα wβ wγ wδ
w8

(B6)

and

(
)
(
)
δαβ − 2ŵα ŵβ∂δδβ ′ − 2ŵδ ŵβ ′4∂2= − 8 [3δβ ′ γ ′ (δαβ wγ + δαγ wβ + δβγ wα )∂wγ ∂wδ
w2
∂wγ ′
w2
w

− (δαβ δγβ ′ + δαγ δββ ′ + δαβ ′ δβγ )wγ ′ − (δαβ δγ γ ′ + δαγ δβγ ′ + δαγ ′ δβγ )wβ ′ ]16[(δαβ wβ ′ wγ − δββ ′ wα wγ − δαβ ′ wβ wγ
w 10+ δαγ wβ wβ ′ + δβγ wα wβ ′ − δβ ′ γ wα wβ )wγ ′+By operating "ββ ′ γ γ ′ , we find"ββ ′ γ γ ′− (δβγ ′ wα wγ + δαγ ′ wβ wγ + δγ γ ′ wα wβ )wβ ′ + 3δβ ′ γ ′ wα wβ wγ ].
()()δαβ − 2ŵα ŵβ∂δδβ ′ − 2ŵδ ŵβ ′12wα∂2=− 8 .
∂wγ ∂wδ
w2
∂wγ ′
w2
w

(B7)

(B8)

After the integration, we obtain Eq. (37).

[1] S. J. Singer and G. L. Nicolson, Science 175, 720 (1972).
[2] B. Alberts, A. Johnson, P. Walter, J. Lewis, and M. Raff,
Molecular Biology of the Cell (Garland Science, New York,
2008).
[3] R. Lipowsky and E. Sackmann, Structure and Dynamics ofMembranes (Elsevier, Amsterdam, 1995).

[4] P. G. Saffman and M. Delbrück, Proc. Natl. Acad. Sci. (USA)72, 3111 (1975).
[5] P. G. Saffman, J. Fluid Mech. 73, 593 (1976).
[6] B. D. Hughes, B. A. Pailthorpe, and L. R. White, J. Fluid Mech.
110, 349 (1981).
[7] E. Evans and E. Sackmann, J. Fluid Mech. 194, 553 (1988).

052407-9

PHYSICAL REVIEW E 95, 052407 (2017)

HOSAKA, YASUDA, OKAMOTO, AND KOMURA[8] S. Ramachandran, S. Komura, M. Imai, and K. Seki, Eur. Phys.
J. E 31, 303 (2010).
[9] K. Seki, S. Ramachandran, and S. Komura, Phys. Rev. E 84,
021905 (2011).
[10] K. Seki, S. Mogre, and S. Komura, Phys. Rev. E 89, 022713(2014).
[11] S. Komura, S. Ramachandran, and M. Imai, in Non-Equilibrium
Soft Matter Physics, edited by S. Komura and T. Ohta (World
Scientific, Singapore, 2012), p. 197.
[12] B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O’Hern, E. R.
Dufresne, and C. Jacobs-Wagner, Cell 156, 183 (2014).
[13] 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).
[14] K. Yasuda, R. Okamoto, S. Komura, and A. S. Mikhailov,
Europhys. Lett. 117, 38001 (2017).
[15] K. Yasuda, R. Okamoto, and S. Komura, Phys. Rev. E 95, 032417(2017).
[16] J.-B. Manneville, P. Bassereau, D. Lévy, and J. Prost, Phys. Rev.
Lett. 82, 4356 (1999).
[17] J.-B. Manneville, P. Bassereau, S. Ramaswamy, and J. Prost,
Phys. Rev. E 64, 021908 (2001).
[18] S. Ramaswamy, J. Toner, and J. Prost, Phys. Rev. Lett. 84, 3494
(2000).
[19] A. S. Mikhailov and R. Kapral, Proc. Natl. Acad. Sci. (USA)112, E3639 (2015).
[20] R. Kapral and A. S. Mikhailov, Physica D 318, 100 (2016).
[21] Y. Koyano, H. Kitahata, and A. S. Mikhailov, Phys. Rev. E 94,
022416 (2016).

[22] M.-J. Huang, A. S. Mikhailov, and H.-Y. Chen, Eur. Phys. J. E
35, 119 (2012).
[23] M.-J. Huang, R. Kapral, A. S. Mikhailov, and H.-Y. Chen,
J. Chem. Phys. 138, 195101 (2013).
[24] M. Tanaka and E. Sackmann, Nature (London) 437, 656(2005).
[25] K. Inaura and Y. Fujitani, J. Phys. Soc. Jpn. 77, 114603(2008).
[26] S. Ramachandran, S. Komura, K. Seki, and G. Gompper, Eur.
Phys. J. E 34, 46 (2011).
[27] S. Ramachandran, S. Komura, K. Seki, and M. Imai, Soft Matter
7, 1524 (2011).
[28] S. Komura, S. Ramachandran, K. Seki, and M. Imai, Adv. Planar
Lipid Bilayers Liposomes 16, 129 (2012).
[29] K. Simons and E. Ikonen, Nature (London) 387, 569 (1997).
[30] S. Komura and D. Andelman, Adv. Colloid Interface Sci. 208,
34 (2014).
[31] D. Lingwood and K. Simons, Science 327, 46 (2010).
[32] H. Diamant, J. Phys. Soc. Jpn. 78, 041002 (2009).
[33] N. Oppenheimer and H. Diamant, Biophys. J. 96, 3041
(2009).
[34] N. Oppenheimer and H. Diamant, Phys. Rev. E 82, 041912(2010).
[35] E. P. Petrov and P. Schwille, Biophys. J. 94, L41 (2008).
[36] E. P. Petrov and P. Schwille, Soft Matter 8, 7552 (2012).
[37] S. Komura, K. Yasuda, and R. Okamoto, J. Phys.: Condens.
Matter 27, 432001 (2015).
[38] A. W. C. Lau and T. C. Lubensky, Phys. Rev. E 80, 011917(2009).

052407-10