PHYSICAL REVIEW E 95, 032417 (2017)

Anomalous diffusion in viscoelastic media with active force dipoles
Kento Yasuda, Ryuichi Okamoto, and Shigeyuki Komura*
Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 26 December 2016; published 27 March 2017)
With the use of the “two-fluid model,” we discuss anomalous diffusion induced by active force dipoles inviscoelastic media. Active force dipoles, such as proteins and bacteria, generate nonthermal fluctuating flowsthat lead to a substantial increment of the diffusion. Using the partial Green’s function of the two-fluid model, we
first obtain passive (thermal) two-point correlation functions such as the displacement cross-correlation function
between the two-point particles separated by a finite distance. We then calculate active (nonthermal) one-point
and two-point correlation functions due to active force dipoles. The time correlation of a force dipole is assumed
to decay exponentially with a characteristic time scale. We show that the active component of the displacement
cross-correlation function exhibits various crossovers from super-diffusive to subdiffusive behaviors dependingon the characteristic time scales and the particle separation. Our theoretical results are intimately related to themicrorheology technique to detect fluctuations in nonequilibrium environment.
DOI: 10.1103/PhysRevE.95.032417
I. INTRODUCTION

The cytoplasm of living cells is full of proteins andorganelles that play important active roles with the aid ofchemical fuels such as adenosine triphosphate (ATP) [1]. Insuch a nonequilibrium environment, the transport propertiesof chemical species drastically deviate from those in static
equilibrium conditions. For example, there are several experimental works reporting the anomalous diffusion of a tagged
particle in biological cells due to protein activities [2–6]. In
other systems, a large enhancement of diffusion was also
observed for a passive particle immersed in a bacterial bath[7,8] or in a suspension of algae Chlamydomonas [9], and sucha phenomenon has been also studied theoretically [10,11].
The modified diffusion in cells was attributed to nonequilibrium forces generated by molecular motors walking on
cytoskeletal networks [12,13]. Recently, Mikhailov and Kapral
proposed a different mechanism caused by nonequilibrium
conformational changes of proteins or enzymes [14,15]. They
showed that, in addition to thermal fluctuations, active proteinsin living cells generate nonthermal fluctuating flows thatlead to a substantial increment of the diffusion constant. A
chemotaxis-like drift of a passive particle was also predicted
when a spatial gradient of active proteins is present [14,15].
In these previous works, however, the three-dimensional(3D) cytoplasm and two-dimensional (2D) biomembrane weretreated as purely viscous fluids characterized by constant shear
viscosities [16].
In general, biological cells behave as viscoelastic materials [17,18]. Hoffman et al. experimentally determined thefrequency-dependent shear modulus of cultured mammaliancells by using various methods to measure their viscoelastic
properties [19,20]. Interestingly, they found two universal(weak) power-law dependencies of the shear modulus at low
frequencies corresponding to the cortical and intracellular
networks. At high frequencies, on the other hand, theyobserved an exponent of 3/4, which was attributed to themechanical response of actin fibers. Such an universal behavior*

komura@tmu.ac.jp

2470-0045/2017/95(3)/032417(14)of mechanical responses in living cells was also reported inother work [21].
Among various methods, microrheology is one of themost useful techniques to measure the rheological propertiesof living cells [22–26]. In this method, the local and bulk
mechanical properties of a single cell can be extracted froma Brownian motion of probe particles, including both thermal
and nonthermal contributions [27]. Concerning its theoreticalbackground, the generalized Stokes-Einstein relation (GSER),
equivalent to the fluctuation dissipation theorem (FDT), hasbeen used to analyze thermal diffusive motions. In nonthermal
situations, the GSER has been further extended to relateparticle mean squared displacement (MSD) and nonthermal
force fluctuations [28,29]. It should be noted, however,that the GSER contains various assumptions that can beviolated in several situations [27]. Therefore, it is necessary
to discuss both thermal and nonthermal Brownian motions ina viscoelastic medium which is described by a well-foundedtheoretical model.
In this paper, we discuss diffusive motion of passive
particles embedded in viscoelastic media that is described by
the “two-fluid model” for gels [30–32]. We especially focus on
the effects of nonthermal fluctuations induced by active force
dipoles which undergo cyclic motions. We calculate displacement cross-correlation functions (CCF) of two-point particles
for the passive situation induced by thermal fluctuations andthe active situation driven by force dipole fluctuations. Our
calculation is closely related to the “two-point microrheology”
method, which has several technical advantages comparedto the “one-point microrheology” [33]. As for the stochastic
property of a force dipole, we consider the case when there isno correlation between different times and also the case when
it decays exponentially with a characteristic time scale. If thedipole time scale is much larger than the viscoelastic timescale, we show that the active contribution of the displacement
CCF exhibits all the possible crossover behaviors between
super-diffusive and subdiffusive motions. Our predictions canbe applied not only for cells but also for bacterial suspensions
and systems containing active colloids.
Since our theory is based on the standard two-fluid model,
it has some similarities to the works by Levine and Lubensky

032417-1©2017 American Physical Society

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURA[34,35] or MacKintosh and Levine [12,13]. In the former
studies [34,35], they investigated the dynamics of rigid spheres
embedded in viscoelastic media by using the two-fluid model,
but did not consider the effects of nonthermal fluctuations. In
the latter studies [12,13], on the other hand, they developed
a model for F-actin networks driven out of equilibrium bymolecular motors. The main difference in our work is that
active force dipoles are embedded in the fluid and exert forceson the fluid itself. In this regard, we use the partial Green’s
function that connects the force acting on the fluid and thefluid velocity as discussed in Refs. [36,37]. In these works,
they emphasized the role of the intermediate length scale in theanalysis of microrheology data. In our separate work, startingfrom the two-fluid model, we have derived effective equationsof motions for tracer particles displaying local deformations
and local fluid flows [38].
In Sec. II, we describe the two-fluid model and show
its partial Green’s function both in the Fourier space andthe real space. In Sec. III, we discuss the passive two-point
correlation functions. Using the coupling mobilities and the
FDT in thermal equilibrium, we calculate the power spectral
density of the velocity CCFs and the displacement CCFs.
In Sec. IV, we shall investigate active one-point correlation
functions due to active force dipoles. We calculate the active
velocity auto-correlation function of a passive point particle
by assuming different time correlations of force dipoles.
We then discuss in Sec. V the active two-point correlation
functions which are useful for two-point microrheology. The
summary of our work and some discussions related to therecent experiments are given in Sec. VI.

μFIG. 1. Schematic representation of the two-fluid model. The
system consists of an elastic network characterized by the Lamé
coefficient μ, and a viscous fluid characterized by the shear viscosity
η. The elastic and fluid components are coupled through the mutual
friction. The length scale  characterizes the typical internal structureof the elastic network. Orange objects represent stochastic force
dipoles which are immersed in the fluid component. Two passive
point particles separated by a distance d are embedded in thefluid component. These passive particles undergo correlated random
Brownian motion due to thermal fluctuations and active stochastic
fluctuations induced by active force dipoles.

conservation,

A. Model description(2)

In the above, ρu and ρv are the mass densities of the twocomponents, μ and λ are the Lamé coefficients of the elastic
network, respectively, η is the shear viscosity of the fluid,p(r,t) is the pressure field, while fu and fv are external force
densities. The elastic and the fluid components are coupledthrough the mutual friction characterized by the friction
coefficient . We note that μ, λ, η, and  are constantsand do not depend on frequency. When the volume fractionof the network is denoted by φ, the above equations arefurther supplemented by the condition of the total volume∂u+ (1 − φ)v = 0.
∇· φ
∂t

(3)

In the following, we employ several simplifications of themodel. (i) We neglect inertial effects, which is justified at
sufficiently low frequencies. Hence the left-hand sides ofEqs. (1) and (2) are both neglected. (ii) We assume thatthe volume fraction of the network is vanishingly small, i.e.,φ  1. In this limit, Eq. (3) can be approximated as∇ · v ≈ 0.

(4)This equation can be regarded as the incompressibility condition of the fluid component.

(1)


∂v
∂u
= η∇ 2 v − ∇p −  v −+ fv .
ρv
∂t
∂t

η

d

II. TWO-FLUID MODEL

To describe viscoelastic media from a general point ofview, we employ the two-fluid model that has been broadly
used to describe the dynamics of polymer gels [30–32,39].
As schematically shown in Fig. 1, there are two dynamical
fields in this model; the displacement field u(r,t) of the elastic
network and the velocity field v(r,t) of the permeating fluid.
Here r is the 3D position vector and t is the time. The coupled
dynamical equations for these two field variables are given by


∂ 2u
∂u
ρu 2 = μ∇ 2 u + (μ + λ)∇(∇ · u) − − v + fu ,∂t∂t

PHYSICAL REVIEW E 95, 032417 (2017)

B. Partial Green’s function

The above linearized equations can be solved by performing
the Fourier transform in space and the Laplace transform intime for any function f (r,t) as defined by
 ∞
 ∞
d 3r
dt f (r,t)e−iq·r−st .
(5)f [q,s] =−∞0Here q is the 3D wave vector and s is the frequency in the
Laplace domain. The general Green’s function (represented by
a 6 × 6 matrix) connecting u and v to fu and fv was calculated
by Levine and Lubensky [34,35].
In this paper, we particularly focus on the response of thefluid velocity v due to the point force fv , and use the partial

032417-2

ANOMALOUS DIFFUSION IN VISCOELASTIC MEDIA . . .

PHYSICAL REVIEW E 95, 032417 (2017)

Green’s function defined byvα [q,s] = Gαβ [q,s]fv,β [q,s].

(a)

(6)

1

Hereafter, the Einstein summation convention over repeatedindices is employed. According to Refs. [34–37], the Green’s
function is given by
1 + (ηb /η)ξ 2 q 2(δαβ − q̂α q̂β ),ηb q 2 (1 + ξ 2 q 2 )G1Gαβ [q,s] =

0.8

(7)

0.4

with q = |q| and q̂ = q/q (see Appendix A for a detailed
derivation). In the above, the frequency-dependent bulk viscosity and characteristic length scale are defined byμμη 1/2,(8)ηb = η + , ξ =
s
sηb

0.2
0
0

respectively. Notice that the above 3 × 3 matrix is nothing but
the part of the general 6 × 6 matrix [35]. In order to study
the effects of molecular motors that generate forces in thecytoskeleton, one needs to take into account fu as discussed
in Refs. [12,13] and recently by us [38]. We note here thatthe partial Green’s function in Eq. (7) does not depend on
the compressional Lamé coefficient λ, while it appears in thegeneral 6 × 6 matrix.
The Green’s function in Eq. (7) can be inverted back fromthe Fourier space to the real space (but remaining in the Laplacedomain). Following the calculation in Appendix A, we obtain
1 − ηb /η1
1+
G1 (r/ξ ) δαβGαβ [r,s] =
8π ηr
ηb /η



1 − ηb /η
(9)
G2 (r/ξ ) r̂α r̂β ,
+ 1+
ηb /η
where r = |r| and r̂ = r/r. Here we have defined the two
scaling functions by


1
1
2
−z
G1 (z) = 1 + 2 − 2e1+ + 2 ,
(10)
z
z z


3
3
6
(11)
G2 (z) = 1 − 2 + 2e−z 1 + + 2 .
zz z
In Fig. 2, we plot both G1 (z) and G2 (z) as a function of z = r/ξ .
When G1 (z) = G2 (z) = 0, the Green’s function Gαβ is purely
determined by η, whereas it is fully described by ηb when
G1 (z) = G2 (z) = 1.

0.6

1

2

4

z

6

8

10

6

8

10

(b)

G2

0.6
0.4
0.2
0
0

z

FIG. 2. Scaling functions (a) G1 and (b) G2 [see Eqs. (10) and (11),respectively] appearing in the partial Green’s function of the two-fluid
model. The scaling variable is z = r/ξ , where r is the distance and ξ
is the frequency-dependent characteristic length scale [see Eq. (8)].
The asymptotic behaviors of these scaling functions, as analyticallygiven by Eqs. (12) and (13), respectively, are plotted with dotted lines.

For our later purpose, we focus here on the large-scale
behavior of the Green’s function. For r  ξ , we obtainGαβ [r,s] ≈These asymptotic behaviors are also plotted in Fig. 2 by thedotted lines, which provide a good approximation especiallyfor z  1.

4

0.8

C. Asymptotic expressions

Next we discuss the asymptotic behaviors of the partial
Green’s function and the scaling functions. We first note thatthe asymptotic expressions of the scaling functions are givenby
4z/3 − 3z2 /4, z  1,
G1 (z) ≈(12)
1 + 2/z2 ,z  1, 2
z  1,z /4,(13)
G2 (z) ≈2
1 − 6/z , z  1.

2sτ1(δαβ + r̂α r̂β )8π ηr 1 + sτ−12(δαβ − 3r̂α r̂β ),4π ηr 3 (1 + sτ )2(14)

where we have introduced the characteristic length and timescales as = (η/ )1/2 ,τ = η/μ.

(15)

As argued in Ref. [36], the first and the second terms of
Eq. (14) are proportional to 1/r and 2 /r 3 , respectively. The
competition between these two terms is characterized by thecrossover length .

032417-3

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 95, 032417 (2017)

D. Coupling mobilities

(a)

0

10

V1x V2x (ω) dIn the following sections, we shall consider correlated
motions of two-point particles embedded in the fluid component. For this purpose, we shall introduce the coupling
mobility between the two points Mαβ [r,s] that is directly
related to the partial Green’s function in Eq. (9). Since
Gαβ is generally expressed as Gαβ = C1 δαβ + C2 r̂α r̂β , the“longitudinal” and the “transverse” coupling mobilities are
given by Mxx = C1 + C2 and Myy = C1 , respectively. Hence,
they are
G1 (z) + G2 (z)
1
1−
,
(16)
Mxx [r,s] =4π ηr2(1 + sτ )1G1 (z)Myy [r,s] =
1−
,
(17)
8π ηr
1 + sτ
√where z = r/ξ = (r/) 1 + sτ . We shall use these coupling
mobilities in order to calculate various correlation functions in
the next sections. Since Mxy = 0 by symmetry, it is sufficient
to consider only the above two coupling mobilities [36].

-1

10

10

2

-2

d/ = 1d/ = 10
= 100

-3

10

10

-4
-4

10
0

10

10

-2

ωτ

10

0

2

10

(b)

V1y V2y (ω) d |-1

III. PASSIVE TWO-POINT CORRELATION FUNCTIONS

Here we discuss the correlated dynamics of two distinctive
passive particles immersed in a viscoelastic gel that is inthermal equilibrium. This situation is relevant to the “twopoint microrheology” experiments as discussed before [33].
Compared to the “single-particle microrheology” (with theuse of a finite size particle), there are several advantages toperform multiparticle microrheology [27].

10
10

-2

-3

10
10

-4

-5

10 -4
10

A. Fluctuation dissipation theorem

Consider a pair of point particles undergoing Brownian
motion separated by a distance d as shown in Fig. 1 (butwithout force dipoles). We denote the positions of these
two-point particles by R1 (t) = R1 + R1 (t) and R2 (t) =R2 + R2 (t), where d = |R2 − R1 |. Then the velocities ofthese point particles are given by V1 (R1 ,t) = Ṙ1 (t) andV2 (R2 ,t) = Ṙ2 (t). The quantities of interest are the velocity
cross-correlation function (CCF) V1α V2α (t) d , and the displacement CCF R1α R2α (t) d . Without loss of generality,
we define the x axis to be along the line connecting the twoparticles, i.e., R2 − R1 = d êx .
According to the fluctuation dissipation theorem (FDT),
the velocity CCFs in thermal equilibrium are related to thecoupling mobility in the Laplace domain by [27,36]V1α V2α [s] d = kB T Mαα [r = d,s],(18)

where kB is the Boltzmann constant, and T is the temperature.
The power spectral density (PSD) of the passive velocity CCF
can be obtained by using the relation
V1α V2α (ω) d = 2 V1α V2α [s = iω] d ,(19)

where ω is the frequency in the Fourier domain, and indicatesthe real part.
In Fig. 3, we plot the scaled PSDs V1x V2x (ω) d andV1y V2y (ω) d as a function of ωτ using the longitudinal and thetransverse coupling mobilities obtained in Eqs. (16) and (17),respectively. Different colors represent different distances, d,

10

-2

ωτ

10

0

2

10

FIG. 3. The passive component of the power spectral density(PSD) of the two-point velocity cross-correlation functions (CCFs)(a) V1x V2x (ω) d and (b) |V1y V2y (ω) d | [see Eq. (19)] as a functionof ωτ for d/ = 1, 10, 100. Here, τ = η/μ is the viscoelastic timescale, and d is the distance between the two-point particles immersed
in viscoelastic media described by the two-fluid model. Both CCFs
are scaled by kB T /(2π ηd) in order to make them dimensionless.
Since V1y V2y (ω) d takes negative values for smaller ωτ (shown by
the dashed lines), we have plotted in (b) its absolute value. The
number “2” in (a) indicates the slope representing the exponent of thepower-law behaviors.

between√ the two points. In Fig. 3(a), V1x V2x (ω) d increases forωτ > 2(/d) as ∼ω2 [see later Eqs. (21) and (22)], while
it saturates for ωτ > 1. Since V1y V2y (ω) d takes negativevalues for smaller ωτ when d/ = 10 and 100, we have plotted
its absolute value in Fig. 3(b). Notice that in Figs. 3(a) and3(b), the PSDs are scaled by kB T /(2π ηd) and take the same
asymptotic value in the large ωτ limit.
The passive displacement CCF in thermal equilibrium as
a function of time can be directly obtained by the following
inverse Laplace transform of the velocity CCF [26,27]:
 c+i∞12R1α R2α (t) d =ds 2 V1α V2α [s] d est ,2π i c−i∞s(20)032417-4

ANOMALOUS DIFFUSION IN VISCOELASTIC MEDIA . . .

ΔR1x ΔR2x (t) d

10

10

4

(a)

d/ = 1d/ = 10= 100

PHYSICAL REVIEW E 95, 032417 (2017)

particles. This is because the elastic property of the medium,
representing the polymer network, is pronounced in these timescales. For d/ = 1, on the other hand, the CCF is almost
proportional to t during the entire time region.
In Fig. 4(b), the absolute value of R1y R2y (t) d isplotted because it takes negative values for larger t. Thismeans that the relative transverse motion of the two-point
particles is anticorrelated when their separation d/ is large
enough. Nevertheless, the general time-dependent behaviorof the transverse CCF is almost the same as that of thelongitudinal one in Fig. 4(a).
The crossover behaviors of the passive displacement CCF
for large d/ showing the successive scaling as t → t 0 → t
can explain some of the apparent power-law behaviors of soft
matter [40] or biological cells [20]. It should be noted, however,that the passive displacement CCF in thermal equilibrium
exhibits only a subdiffusive behavior.

1
1

2

1
10

0

1
-2

10

0

10

2

4

10

10

ΔR1y ΔR2y (t) d |

t/τ
10

10

4

(b)

B. Large distance behaviors

1

2

In the limit of large distances d   between the two points,
we can use Eq. (14) for the partial Green’s function to obtain
the PSDs in the Fourier domain as1V1x V2x (ω) d ≈

1
10

0

1

(ωτ )2kB T
2π ηd 1 + (ωτ )2
1
kB T 2(ωτ )2,+−π ηd 3 [1 + (ωτ )2 ]2[1 + (ωτ )2 ]2
(21)

-2

10

0

10

2

10

4

10

V1y V2y (ω) d ≈t/τ
FIG. 4. The passive component of the two-point displacement
cross-correlation functions (CCFs) (a) R1x R2x (t) d and (b)|R1y R2y (t) d | [see Eq. (20)] as a function of t/τ for d/ = 1, 10,100. Here d is the distance between the two-point particles immersed
in viscoelastic media. Both CCFs are scaled by kB T τ/(2π ηd) in orderto make them dimensionless. Since R1y R2y (t) d takes negativevalues for larger t/τ (shown by the dashed lines), we have plotted in(b) its absolute value. The numbers indicate the slope representing
the exponent of the power-law behaviors.

where c is a real number. Performing the numerical inverse Laplace transform of Eq. (20), we plot in Fig. 4 thelongitudinal and the transverse CCFs R1x R2x (t) d andR1y R2y (t) d as a function of t/τ for different distances dbetween the two points.
In Fig. 4(a), the longitudinal CCF is proportional to t fort/τ < 1 and t/τ > (d/)2 /2 for d/ = 10 and 100 [see later
Eq. (23)]. In the former time region, which is smaller thanthe viscoelastic time scale τ = η/μ, the two-point particles
interact through the fluid component of the two-fluid model.
In the latter long-time region, on the other hand, the CCF obeys
the normal diffusive behavior as expected for any viscoelastic
material with a characteristic relaxation time. Between thesecrossover time scales, the CCF remains almost constant due tothe elastic component that suppresses the motion of the point(ωτ )2kB T
4π ηd 1 + (ωτ )2
1
kB T 2(ωτ )2.
−−2π ηd 3 [1 + (ωτ )2 ]2[1 + (ωτ )2 ]2(22)

For the large distance behavior of the displacement CCFs,
we obtain
R1x R2x (t) d ≈kB T τ
(1 − e−t/τ )
2π ηd



kB T τ 2 t
−t/τ
−t/τ
) − 2(1 − e) ,+(1 + eπ ηd 3 τ
kB T τ
(1 − e−t/τ )R1y R2y (t) d ≈
4π ηd


kB T τ 2 t
−t/τ
−t/τ
(1
+
e
)
−
2(1
−
e
)
.
−
2π ηd 3 τ

(23)

(24)

In the above expressions, the first term is proportional to t in the
short-time regime, whereas it saturates in the long-time limit.
Whereas the second term in each expression is proportional
to t in the long-time limit, which dominates the large-scale
behavior. These properties of the displacement CCF can beclearly observed in Fig. 4 especially for larger d/. Although
not plotted, Eqs. (23) and (24) almost completely recover thenumerical plots for d/ = 10 and 100 in Fig. 4.

032417-5

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURA
IV. ACTIVE ONE-POINT CORRELATION FUNCTIONSIn this section, we shall consider the collective advection
effects due to active force dipoles on passive particles
immersed in viscoelastic media. A simple “dimer model”
for a stochastic hydrodynamic force dipole was previously
discussed in Refs. [14,15]. To investigate the hydrodynamic
effects of the force dipoles, we employ the partial Green’s
function representing the response of the fluid velocity v due tothe force fv acting on the fluid component [see Eq. (6)]. This isdifferent from Refs. [12,13] where they discussed the effects ofmolecular motors generating forces in the cytoskeleton, whichcorresponds to the elastic component of the two-fluid model.
Our aim is to focus on the role of active force dipoles that exist
in the fluid component. A more unified treatment of these twodifferent sources of active forces has been investigated in our
separate publication [38].

[cf. Eq. (5)]. Performing the bilateral Fourier transform ofEq. (27), we obtain∂Gαβ [r,ω]Vα (R0 ,ω) ≈ d 3 r
∂rγ
×

−∞

Consider an oscillating dimer of length a(t) and the force
magnitude fd (t) with its orientation given by the unit vector ê.
In this case, the induced velocity of a passive particle at R due
to the oscillating dimer is given by [14] t∂Gαβ (R − r,t − t )êβ êγ m(t ), (26)Vα (R,t) ≈
dt
∂rγ
−∞
where we have used the approximation a  |R − r|, andm(t) = a(t)fd (t) denotes the magnitude of the force dipole.
We further consider a collection of such active force dipoles,
located at positions {Ri } with orientations {êi }. By summing
up for all the dipoles, the velocity of the passive particle isthen given by [14] t∂Gαβ (r,t − t )dtVα (R0 ,t) ≈
d 3r
∂rγ
−∞
×

êi,γ êi,β mi (t )δ(Ri − R0 − r), (27)iwhere we have assumed that the displacement of the passive
particle is small, and kept only the lowest-order term. Thisequation describes the motion of a passive point particle due
to nonthermal active noise arising from the collective operationof active force dipoles.
Hereafter we introduce the bilateral Fourier transform intime for any function f (t) as
 ∞dt f (t)e−iωt ,(28)
f (ω) =−∞

êi,γ êi,β mi (ω)δ(Ri − R0 − r),

(29)

i

where Gαβ [r,ω] = Gαβ [r,s = iω]. We shall use this expression to calculate the velocity correlation functions and themean squared displacement (MSD) of the passive particle.
B. Active autocorrelation functions

We now calculate the velocity autocorrelation function(ACF) of a passive particle located on average at R0 . If
the correlation between different force dipoles vanishes, i.e.,
mi mj (ω) = m2 (ω) δij , we get from Eq. (29)Vα Vα (R0 ,ω)∂Gαβ [r,ω] ∂Gα β [r,−ω]= d 3r∂rγ∂rγ

A. Velocity induced by active force dipoles

When a point force fv is applied to the fluid at a point r,
it induces a fluid velocity at another position R that advects
a point particle located there. As in the previous section, we
denote the position of this passive point particle by R(t) =R0 + R(t), and its velocity by V(R,t) = Ṙ(t). Using theGreen’s function calculated in Sec. II, we obtain the relation
between V and fv as
 t
Vα (R,t) =dt Gαβ (R − r,t − t )fv,β (r,t ).
(25)

PHYSICAL REVIEW E 95, 032417 (2017)×êi,γ êi,β êi,γ êi,β m2i (ω) δ(Ri − R0 − r)i

= ββ γ γ m2 (ω)
∂Gαβ [r,ω] ∂Gα β [r,−ω]c(R0 + r),
× d 3r
∂rγ
∂rγ

(30)

where c(r) = i δ(Ri − r) is the local concentration of force
dipoles at a point r in the fluid component. In the above, asymbol,ββ γ γ = êβ êβ êγ êγ1(δββ δγ γ + δβγ δβ γ + δβγ δβ γ ), (31)15
has been defined, and we have assumed that the orientationsof active force dipoles are not correlated with their positions.
In other words, we do not consider any nematic ordering offorce dipoles [14].
When active force dipoles are uniformly distributed in spacewith a constant concentration, c(r) = c0 , the velocity ACF
Vα Vα (ω) is isotropic, i.e., Vx2 (ω) = Vy2 (ω) = Vz2 (ω) ,and vanishes otherwise. Hence it is enough to consider only
the x direction, and we obtain the active PSD as
c0
Vx2 (ω) = ββ γ γ m2 (ω)
3∂Gαβ [r,ω] ∂Gαβ [r,−ω]× d 3r
∂rγ
∂rγ
=

=

c0
1
m2 (ω) I(ω).
3 × 82 × 15π 2 η2 

Here we have introduced the scaled PSD defined by∂gαβ [r,ω] ∂gαβ [r,−ω],I(ω) = 15ββ γ γ
d 3 r̄
∂ r̄γ
∂ r̄γ

(32)

(33)

together with gαβ [r,ω] = 8π ηGαβ [r,s = iω] and r̄ = r/.
Since the above integral diverges for short length scales, weneed to introduce a small cutoff length δ. Physically, δ can

032417-6

ANOMALOUS DIFFUSION IN VISCOELASTIC MEDIA . . .

PHYSICAL REVIEW E 95, 032417 (2017)

10

2

6

10

(a)

δ/ = 1
δ/ = 10
δ/ = 100

5

10
0

10

10

10

(ΔRx )2 (t)Vx2 (ω)

10

2

-2

=1
δ/ = 10
δ/ = 100

-4

10

10

-4

10

-2

10

ωτ

0

10

10

2

1

1

2

10 -1
10

10

0

1

10

2

10

10

3

4

10

10

t/τ

FIG. 5. The active component of the power spectral density (PSD)Vx2 (ω) [see Eq. (32)] as a function of ωτ for δ/ = 1, 10, 100. Herea single point particle is immersed in viscoelastic media described by
the two-fluid model, and δ is the cutoff length corresponding to theparticle size. In the plot, the PSD is scaled by c0 m2 (ω) /(2880π 2 η2 δ)in order to make it dimensionless. The number indicates the slope
representing the exponent of the power-law behavior.

The detailed derivation of this expression is given in
Appendix B. Although not plotted, we have confirmed thatEq. (34) perfectly reproduces the curve of δ/ = 100 in Fig. 5.
We should keep in mind, however, that to regard the cutoff
length δ as the particle size is only an approximation, andhence the numerical prefactor should not be taken as accurate
when we compare with experiments.
C. Uncorrelated force dipoles

In order to further calculate the active PSD, the statisticalproperty for the time correlation of a force dipole needs to be
specified. First we assume that it is only δ-correlated in timeand is given by
(35)

6

10

(b)

1

4

10

(ΔRx )2 (t)be regarded as the size of the passive particle. In Fig. 5, we
numerically plot the scaled PSD, I(ω), as a function of ωτfor different cutoff lengths δ/ = 1, 10, and 100. These values
correspond to the situation when√the passive particle is larger
than the mesh size. For ωτ > 2 3(/δ)2 , the PSD increasesas ∼ω2 and saturates for ωτ > 1 [see later Eq. (34)]. Notethat the asymptotic value of the scaled PSD for large ωτ isindependent of δ.
If we use Eq. (14) for the partial Green’s function Gαβin the large distance limit r  , the scaled PSD can beapproximately calculated as  3(ωτ )24(ωτ )2
−
I(ω) ≈ 48π
1 + (ωτ )2 δ
[1 + (ωτ )2 ]2 δ 5 

12
+
.
(34)
[1 + (ωτ )2 ]2 δm(t)m(t ) = Sδ(t − t ),

1

3

10

-6
-6

1

4

1

2

10

2
0

10

1

=1
δ/ = 10
δ/ = 100

-2

10

2

-4

10

-1

10

0

10

10

1

2

10

3

10

4

10

t/τ
FIG. 6. The active component of the mean squared displacement(MSD) (Rx )2 (t) [see Eq. (36)] as a function of t/τ for δ/ = 1, 10,100. Here a single point particle is immersed in viscoelastic media,
and δ is the cutoff length. (a) The case when the time correlationof the force dipole is δ-correlated [see Eq. (35)]. (b) The case when
the time correlation of the force dipole decays exponentially with acharacteristic relaxation time τd [see Eq. (39)], and we set here τd /τ =100. In these plots, (Rx )2 (t) is scaled by c0 Sτ/(2880π 2 η2 δ) inorder to make it dimensionless. The numbers indicate the slope
representing the exponent of the power-law behavior.

where S fixes the fluctuation amplitude. In the Fourier
representation, this simply means that m2 (ω) = S. Once weknow the active PSD of the velocity ACF, the corresponding
MSD of a passive particle in the x direction can be obtainedby the inverse Fourier transform: ∞dω 2(Rx )2 (t) =V 2 (ω) eiωt .
(36)
2 x−∞ 2π (iω)
In contrast to the inverse Laplace transform in Eq. (20), we also
take into account the initial condition by including a constantterm in the above transformation. In Fig. 6(a), we numericallyplot the scaled (Rx )2 (t) as a function of t/τ for δ/ = 1,
10, and 100. Here the MSD is proportional to t both for shorttime scales t/τ < 1 and for long time scales t/τ > (δ/)4 /12032417-7

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURA[see later Eq. (38)]. For the intermediate time range 1 < t/τ <(δ/)4 /12, on the other hand, the MSD is strongly suppressed
due to the elastic component of the two-fluid model, and it
remains almost constant.
If we use the asymptotic expression Eq. (14) for the partial
Green’s function as before, the PSD can be obtained from
Eqs. (32) and (34) as
 
 3


4(ωτ )2
1 c0 S(ωτ )22
Vx (ω) ≈−222
260π η  1 + (ωτ ) δ[1 + (ωτ ) ] δ 5 

12
+
.
(37)
2
2
[1 + (ωτ ) ] δ
Then, with the use of Eq. (36), the asymptotic MSD can beanalytically obtained as 
1 c0 Sτ

−t/τ
)
(1
−
e
(Rx )2 (t) ≈2
60π η δ 3
+ 2[e−t/τ (t/τ + 1) − 1]δ 5 
−t/τ
+ 6[e(t/τ + 3) + 2t/τ − 3]
. (38)
δThe first term in the right-hand side of the above equation
indicates that the normal diffusion occurs for the short timescale t  τ , while it saturates in the longer times. In the longtime limit, on the other hand, we can set e−t/τ ≈ 0, and one
finds that MSD is proportional to t for t/τ  (δ/)4 /12 [seethe third line of Eq. (38)], as mentioned above.

Using the asymptotic expression Eq. (14), the active PSD
is now given by (ωτ )21

1 c0 SVx2 (ω) ≈2
2
2
60π η  1 + (ωτd ) 1 + (ωτ ) δ 3
 5 4(ωτ )2
12
−
+
.
2
2
2
2
[1 + (ωτ ) ] δ[1 + (ωτ ) ] δ(41)Then the corresponding active one-point MSD can be obtained
up to the lowest order in /δ as1
1 c0 Sτ(Rx )2 (t) ≈60π η2  1 + τd /τ
 

e−t/τ

(τd /τ )e−t/τd−. (42)× 1+
1 − τd /τ
1 − τd /τδThis equation reduces to the first line of Eq. (38) in the limit of
τd → 0. By Taylor expanding the above expression for small t,
one can indeed show that the linear term in t vanishes, and theactive MSD increases as ∼t 2 . The full expression of the active
one-point MSD including higher-order terms is provided in the
Appendix C. The analytic expressions in Eqs. (41) and (42)are the general and important results of this paper.
V. ACTIVE TWO-POINT CORRELATION FUNCTIONS
A. Velocity cross-correlation functions

In this section, we consider the active velocity CCF between
the two points at R1 and R2 that are separated by a distance d,as shown in Fig. 1 and also discussed in Sec. III. With the use of
Eq. (27), the active two-point velocity CCF can be evaluated by
V1α V2α (R1 ,R2 ,ω) d
3 ∂Gαβ [r,ω] ∂Gα β [r ,−ω]
= d r
∂rγ
∂rγ

D. Exponentially correlated force dipoles

Next we calculate the PSD and the MSD when the timecorrelation of a force dipole decays exponentially with acharacteristic relaxation time τd , i.e.,m(t)m(t ) =S −|t−t |/τde.
2τdr =r−(R2 −R1 )

×

(39)

1
c0 S
1
I(ω),
2
2
2880π η  1 + (ωτd )2(40)

where I(ω) was defined before in Eq. (33).
In Fig. 6(b), we numerically plot the scaled (Rx )2 (t) asa function of t/τ when τd /τ = 100 for δ/ = 1, 10, and 100,
i.e., the distance between the two points is larger than the meshsize. For δ/ = 1, we find that the active MSD is proportional
to t 2 and exhibits a super-diffusive behavior within the time
region t < τd . For δ/ = 100, such a super-diffusive behavior
is observed only up to the viscoelastic time scale t/τ < 1, andthe MSD exhibits a normal diffusive behavior for t/τ > 1.
The active MSD for δ/ = 100 is further suppressed for larger
time scales. In the very long time limit, the active MSD will
be again proportional to t [38].

êi,γ êi,β êi,γ êi,βm2i (ω) δ(Ri − R1 − r)i

= ββ γ γ m2 (ω)



∂G
[r,ω]
[r
,−ω]
∂G
αβ
α
β
× d 3r∂rγ∂rγ

In this case, we have m2 (ω) = S/[1 + (ωτd )2 ] in the Fourier
representation. Some justification of the above simple expression will be separately discussed in Sec. VI. Mathematically,
Eq. (39) reduces to Eq. (35) in the limit of τd → 0. Then theactive PSD is given by
Vx2 (ω) =PHYSICAL REVIEW E 95, 032417 (2017)

r =r−(R2 −R1 )× c(R1 + r).

(43)

As before, we can generally set R2 − R1 = d êx without lossof generality. We also assume that the active force dipoles are
uniformly distributed in space with a constant concentration,
c0 . Then one can further rewrite asc01m2 (ω) Iαα (d,ω), (44)V1α V2α (ω) d = 2
2
8 × 15π η2 
where
Iαα (d,ω)= 15ββ γ γ∂gαβ [r,ω] ∂gα β [r ,−ω]
d r̄
,
∂ r̄γ
∂ r̄γ
r =r−d3(45)
with gαβ = 8π ηGαβ as defined before.

032417-8

ANOMALOUS DIFFUSION IN VISCOELASTIC MEDIA . . .

PHYSICAL REVIEW E 95, 032417 (2017)

2

10

V1x V2x (ω) d

10

(a)

-2

10

ΔR1x ΔR2x (t) d

0

5

10

2

-4

10

d/ = 1.1= 10
d/ = 100

-6

10

(a)

= 1.1
d/ = 10
= 100

4

10

1

3

10

2

10

1

10

1

-8

10

0

-4

-2

10

10

0

ωτ

10

10 -1
10

2

10

4

10

ΔR1y ΔR2y (t) d |(b)0

V1y V2y (ω) d |

10

-2

10

2

-4

10

10

1

2

10

3

4

10

10

t/τ

2

10

0

10

-6

10

(b)
1

3

10

2

10

1

1

10

1

-8

10

0

-4

10

-2

10

0

ωτ

10

10 -1
10

2

10

In Fig. 7, we numerically plot the scaled active PSDs
V1x V2x (ω) d and V1y V2y (ω) d as a function of ωτ fordifferent distances d/ = 1.1, 10 and 100, as before. (Thereason that we chose here d/ = 1.1 is that there wasa numerical stability issue exactly at d/ = 1.) The PSD
increases as ∼ω2 for the intermediate frequency range. Within
the lowest-order term in Eq. (14), the asymptotic expressionsof the active PSDs can be obtained as
16π  (ωτ )2,
d 1 + (ωτ )2

Iyy (d,ω) = Izz (d,ω) ≈8π  (ωτ )2.
d 1 + (ωτ )2

10

1

2

10

3

10

4

10

t/τ

FIG. 7. The active component of the scaled power spectral density(PSD) (a) V1x V2x (ω) d and (b) |V1y V2y (ω) d | [see Eq. (44)] as
a function of ωτ for d/ = 1.1, 10, 100. Here d is the distance
between the two-point particles immersed in viscoelastic media. Both
PSDs are scaled by c0 m2 (ω) /(960π 2 η2 d) in order to make themdimensionless. Since V1y V2y (ω) d takes negative values for smallerωτ (shown by the dashed lines), we have plotted in (b) its absolutevalue. The numbers indicate the slope representing the exponent of
the power-law behaviors.

Ixx (d,ω) ≈

0

10

(46)

(47)

FIG. 8. The active component of the two-point displacement
cross-correlation functions (CCFs) (a) R1x R2x (t) d and (b)|R1y R2y (t) d | as a function of t/τ for d/ = 1.1, 10, 100.
Here d is the distance between the two-point particles immersed inviscoelastic media, and the time correlation of the force dipole is δcorrelated [see Eq. (35)]. Both CCFs are scaled by c0 Sτ/(960π 2 η2 d)in order to make them dimensionless. Since R1y R2y (t) d takes
negative values for larger t/τ (shown by the dashed lines), we haveplotted in (b) its absolute value. The numbers indicate the slope
representing the exponent of the power-law behaviors.
B. Displacement cross-correlation functions

Performing the inverse Fourier transform of the active
two-point PSDs as before [see Eq. (36)], we obtain thecorresponding longitudinal and transverse displacement CCFs
R1x R2x (t) d and R1y R2y (t) d for the distances d/ =1.1, 10, and 100. In Figs. 8(a) and 8(b), we plot these quantities
when the time correlation of a force dipole is δ-correlated
as assumed in Eq. (35). Figure 8 should be compared with
Fig. 6(a) where we have shown the MSD for the active
one-point case. Both longitudinal and transverse displacement
CCFs are proportional to t for short time scales t/τ < 1 andalso for longer time scales. For the intermediate time range,
however, these CCFs are strongly suppressed and become
constant due to the elastic component of the two-fluid model.

032417-9

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURA(a)4

ΔR1x ΔR2x (t) d

102

VI. SUMMARY AND DISCUSSION

1

= 1.1
d/ = 10
= 100

10

1

0

2

10

1

-2

10

-4

2

10

-1

0

10

1

10

10

2

10

3

4

10

10

t/τ
4

ΔR1y ΔR2y (t) d |

10

(b)

1

2

10

1

0

10

2

1

-2

10

-4

10

2
-1

10

0

10

1

10

t/τ

2

10

PHYSICAL REVIEW E 95, 032417 (2017)

3

10

4

10

FIG. 9. The active component of the two-point displacement
cross-correlation functions (CCFs) (a) R1x R2x (t) d and (b)|R1y R2y (t) d | as a function of t/τ for d/ = 1.1, 10, 100.
Here d is the distance between the two-point particles immersed
in viscoelastic media. The time correlation of the force dipole
decays exponentially with a characteristic relaxation time τd [seeEq. (39)], and we set here τd /τ = 100. Both CCFs are scaled
by c0 Sτ/(960π 2 η2 d) in order to make them dimensionless. Since
R1y R2y (t) d takes negative values for larger t/τ (shown by thedashed lines), we have plotted in (b) its absolute value. The numbers
indicate the slope representing the exponent of the power-law
behaviors.

In Figs. 9(a) and 9(b), on the other hand, we consider thecase when the time correlation of a force dipole is characterized
by a relaxation time τd /τ = 100 [see Eq. (39)]. These figures
should be compared with Fig. 6(b) because the overallbehavior is similar. For δ/ = 1.1, the active displacement
CCFs are proportional to t 2 when t < τd , showing a strongsuper-diffusive behavior. For δ/ = 100, however, this superdiffusive behavior is observed only within the time region
smaller than the viscoelastic time scale, t/τ < 1, and the CCFs
increase as ∼t for t/τ > 1. For much longer time scales, theactive CCFs are further suppressed because of the elasticity.
In the long time limit, the active CCFs are both proportional
to t.

In this paper, we have discussed anomalous diffusion
induced by active force dipoles in viscoelastic media that is
described by the standard two-fluid model for gels. We firstreviewed the two-fluid model and showed its partial Green’s
function both in the Fourier and the real spaces. With the use of
the coupling mobilities and the FDT in thermal equilibrium,
we have calculated the PSD of the velocity CCFs and thedisplacement CCFs between the two-point particles both forthe longitudinal and the transverse directions. The obtained
results are useful to interpret the data obtained by two-point
microrheology experiments. The passive (thermal) two-point
CCF increases linearly with time at shorter and longer timescales, while it is suppressed and remains almost constant atintermediate time scales (see Fig. 4).
Moreover, we have calculated active (nonthermal) onepoint and two-point correlation functions due to active force
dipoles. We have used the relation between the velocity andthe dipole strength, as given by Eq. (26), and the formulation
in Ref. [14] in order to further calculate the active PSD of thevelocity CCFs. For the one-point case, one needs to introduce
a cutoff length scale, δ, in evaluating the integrals, whereas afinite distance, d, between the two-point particles plays the roleof the cutoff length in the two-point case. As for the statisticalproperty of force dipoles, we considered the case when theirmagnitude is uncorrelated in time [see Eq. (35)] and the casewhen it decays exponentially with a characteristic time τd [seeEq. (39)].
For the active case, the important results can be summarized
as follows. As shown in Fig. 6(b) (one-point case) or Fig. 9(two-point case), we have found that the active MSD or
the displacement CCFs exhibits various crossovers fromsuper-diffusive to subdiffusive behaviors depending on thecharacteristic time scales (τ = η/μ and τd ) and the particle
separation d (or the cutoff length δ for the one-point case). We
emphasize that the active displacement CCF is proportional tot 2 for time scales shorter than the viscoelastic time scale, t < τ ,and it is proportional to t for the intermediate time scales,τ < t < τd . Within the present model, the passive contribution
only describes subdiffusion, whereas the active contribution
is responsible for both subdiffusion and super-diffusion. Our
results are useful in understanding active properties of thecytoplasm using force spectrum microscopy combined withthe microrheology experiment [6], as further discussed below.
In Ref. [6], Guo et al. measured the MSD of microinjectedtracer particles in mellanoma cells. They showed that the MSD
was nearly constant at shorter time scales (t < 0.1 s), while
it exhibited a slightly super-diffusive behavior at longer timescales (t > 0.1 s), i.e., (R)2 ∼ t β with β ≈ 1.2. However,when they inhibited motor and polymerization activity bydepleting cells of ATP, the MSD was almost constant in
time, i.e., β ≈ 0. Such an ATP-dependent Brownian motion
was also observed in prokaryotic cells and yeast [4,5]. In
addition to the MSD measurement, Guo et al. performed
active microrheology experiment [24,25], and found that thefrequency-dependent elastic modulus follows a power-law
form, i.e., |G(ω)| ∼ ωα with α ≈ 0.15 [6].
For simplicity, one may assume that PSD of the active
force also obeys a power-law behavior, i.e., m2 (ω) ∼ ω−γ

032417-10

ANOMALOUS DIFFUSION IN VISCOELASTIC MEDIA . . .

PHYSICAL REVIEW E 95, 032417 (2017)with a different exponent γ . Among these three exponents, thefollowing scaling relation should hold [6,28]:
β = 2α + γ − 1.

(48)

In thermal equilibrium, γ = −α + 1 holds according to the
FDT and hence β = α. In this case, the anomalous diffusion
purely reflects the viscoelasticity of the surrounding media.
The exponentially correlated force dipoles in Eq. (39) leads
to the active PSD in Eq. (40), and hence γ = 2 for ωτd  1.
Experimentally, the value γ ≈ 2 was observed by Lau et al.
[28] and later reconfirmed by Guo et al. [6]. When α ≈ 0.15
and γ ≈ 2, Eq. (48) gives β ≈ 1.3 which is almost consistent
with the MSD measurement mentioned above (β ≈ 1.2). In
the older experiment [28], on the other hand, the measured
exponents were α ≈ 0.25, γ ≈ 2 and hence β ≈ 1.5. In bothof these experiments, they claimed that active forces dominate
the low-frequency regime, whereas thermal forces dominate
the high-frequency regime [6,28].
It should be reminded, however, that different values ofγ were reported by different groups [41–43]. For example,
a combination of active and passive microrheology measurements using PC3 tumor cells resulted in α ≈ 0.4, β ≈ 1.3, andγ ≈ 1.5, satisfying also the scaling relation Eq. (48). They
argued that such a difference can arise because active andpassive measurements were done in Ref. [28] with different
probes and at very different locations in the cell. In Refs. [41–43], they performed dual passive-active measurements witha unique probe. Given these situations, we consider that thepower-law behavior of the force fluctuations and its exponent
require further experimental and theoretical investigations.
We also point out that Eq. (48) cannot be always true
because the exponents can take values only 0  β  2 and0  α  1. Hence, if γ = 2 holds, one should always observea super-diffusive behavior because β = 2α + 1  1. However,subdiffusive behaviors (β < 1) in cells have been observed inmany cases [44–46]. Moreover, the above relation also restricts
the value of α to 0 < α  0.5 because β < 2, which is notalways the case [19,20].
In our work, we have assumed that the time correlation of
a force dipole is an exponentially decaying function with acharacteristic time τd , as given in Eq. (39). Hence its Fouriertransform has a Lorentzian form, and decays as ω−2 forωτd  1. A similar Lorentzian form of force fluctuations was
discussed by Levine and MacKintosh [12,13]. While someof the experiments that reported the exponent γ = 2 [6,28]
justify our assumption, different values of γ found in theother experiments [41–43] indicate that the dipole correlation
cannot be a simple exponentially decaying function. Hencea more detailed investigation for the statistical property of afluctuating force dipole is required. Currently, we are analyzing
the stochastic properties of a simple model of a catalytically
active bidomain protein [14]. In this model, the two protein
domains are represented by beads connected by an elastic
spring, and the two internal states, namely, free protein andligand-protein complex, are assumed.
Although our theory is general and can be applied not
only for cells but also for other macroscopic systems, it isuseful to give some typical parameter values corresponding
to a cell. Since the characteristic length scale  = (η/ )1/2
roughly corresponds to the mesh size of a polymer gel, it isroughly given by  ∼ 10−7 m for a typical cell. Hence thedistance between the two-point particles such as d/ = 100
means d ∼ 10−5 m. According to Ref. [6], we also haveη ∼ 10−3 Pa · s and μ ∼ 1 Pa so that the viscoelastic timescale can be estimated as τ = η/μ ∼ 10−3 s. This meansthat the dipole time scale τd /τ = 100 used in Figs. 6(b) and9 corresponds to τd ∼ 10−1 s. Although this time scale issomewhat larger than the cycle time of a single protein machine[14], it still gives a good estimate to characterize the collective
dynamics of a protein complex.
Recently, Fodor et al. [47] have made an attempt to theoretically reproduce the MSD data measured in the cytoplasm ofliving A7 [6]. They used one-dimensional Langevin equationin the presence of a random active force to calculate both thethermal and nonthermal MSD. Their theory has a similarity
to the present work because they also introduce two timescales that are analogous to τ = η/μ and τd in our theory.
An important new aspect in the present paper is that theinternal structure of the viscoelastic medium is properly taken
into account. Both thermal and nonthermal MSDs exhibit
complicated time sequences depending on the length-scaleof the observation relative to the mesh size . In Ref. [47], thesize of tracer particles was assumed to be always larger thanthe mesh size of the cytoskeletal network.
In our separate work, we have considered the two-fluid
model where active macromolecules, described as force
dipoles, cyclically operate both in the elastic and the fluid
components [38]. Through coarse-graining, we have derived
effective equations of motions for tracer particles displaying
local deformations and local fluid flows. The equation for
deformation tracers coincides with the phenomenologicalmodel by Fodor et al. [47] (see also the related publication[48]). Our analysis reveals that localization and diffusion
phenomena are generally involved. The motion of tracers
immobilized within the elastic subsystem is localized in thelong-time limit, but it can show a diffusion-like behavior
at the intermediate time scales shorter than the cooperative
correlation times of molecular motor aggregates operating in
the active gels [38].
Recently, Bruinsma et al. [49] investigated a large-scalecorrelated motion of chromatin inside the nuclei of living cells
by using another “two-fluid model” for polymer solutions[50] (but not for gels). They derived the response functions
that connect the chromatin density and velocity correlation
function to the correlation functions of the active sources thatare either scalar or vector quantities. One of the differences intheir theory is that the form of the complex viscoelastic moduli
needs to be specified in order to compare with experiments,
whereas the viscoelasticity naturally arises from the present
two-fluid model. It is interesting to note that their active PSD
also contains Lorentzian-type frequency dependence as wehave obtained such as in Eq. (40). It would be interesting
to calculate the active MSD based on this different two-fluid
model.
Finally, we mention that anomalous diffusion observed incolloidal gels has been also explained in terms of force dipoles
due to structural inhomogeneities [51,52]. Assuming that such
inhomogeneities are randomly distributed, it was shown thatthe relaxation time of the dynamic structure factor is inversely
proportional to the wave number. In Ref. [53], the MSD

032417-11

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 95, 032417 (2017)exhibits diffusive motion at short times and super-diffusive
motion at long times.

the Laplace transform in time, Eqs. (1), (2), and (4) can be
represented in the steady state as

ACKNOWLEDGMENTS

0 = −μq 2 u[q,s] − (μ + λ)q(q · u[q,s])

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

− (su[q,s] − v[q,s]),(A1)

0 = −ηq 2 v[q,s] − iqp[q,s]− (v[q,s] − su[q,s]) + fv [q,s],(A2)

q · v[q,s] = 0.

(A3)

Taking the inner products of both Eqs. (A1) and (A2) with q,
and using Eq. (A3), we obtain

APPENDIX A: PARTIAL GREEN’S FUNCTION

In this Appendix, we show the derivation of Eqs. (7)and (9) [34–37]. By using the Fourier transform in space andp[q,s] = −iq · fv [q,s]
.
q2

(A4)

From Eq. (A1), we can solve for u as(μ + λ)q 2
δαβ −q̂α q̂β vβ [q,s].
uα [q,s] =μq 2 + s
(μq 2 + s)(2μq 2 + λq 2 + s)(A5)

Substituting Eqs. (A4) and (A5) into Eq. (A2), we obtain the following equation:
 2q [ημq 2 + (sη + μ)]s 2 (μ + λ)q 2
δ
q̂
q̂
−
αβ
α β vβ [q,s] = (δαβ − q̂α q̂β )fv,β [q,s].
μq 2 + s
(μq 2 + s)(2μq 2 + λq 2 + s)(A6)Then we can solve for v asvα [q,s] =μq 2 + s(δαβ − q̂α q̂β )fv,β [q,s].
q 2 [ημq 2 + (sη + μ)](A7)

In terms of ηb and ξ defined in Eq. (8), we finally obtain Eq. (7).
Next, we derive the real-space representation of the partial Green’s function [54]. We first assume that it has the form ofGαβ [r,s] = C1 δαβ + C2 r̂α r̂β ,(A8)Gαα [r,s] = 3C1 + C2 ,(A9)Gαβ [r,s]r̂α r̂β = C1 + C2 .

(A10)

so thatHence we have
1 − ηb /η
d 3 q 1 + (ηb /η)ξ 2 q 2 iq·r
1
−r/ξ
e
1+
(1 − e3C1 + C2 = 2
=
) ,
(A11)
(2π )3 ηb q 2 (1 + ξ 2 q 2 )
2π ηr
ηb /η




d 3 q 1 + (ηb /η)ξ 2 q 2
1
1 − ηb /η 
2 iq·r2−r/ξ2
C1 + C2 =1−
2(ξ/r)]e=+2e[(ξ/r)
+
(ξ/r)
]
[1
−
(
q̂
·
r̂)
1
+
.
(2π )3 ηb q 2 (1 + ξ 2 q 2 )
4π ηr
ηb /η


(A12)
Solving for C1 and C2 , we finally arrive at Eq. (9) with Eqs. (10) and (11).

APPENDIX B: DERIVATION OF EQ. (34)

In this Appendix, we show the derivation of Eq. (34). Here we use the dimensionless form of the Green’s function gαβ =8π ηGαβ , and consider its asymptotic expression [see Eq. (14)]gαβ [r,ω] =

iωτ 
23
(δαβ − 3r̂α r̂β ) ≡ Aαβ (r,ω) − Bαβ (r,ω),(δαβ + r̂α r̂β ) − 3r(1 + iωτ )r (1 + iωτ )2032417-12(B1)

ANOMALOUS DIFFUSION IN VISCOELASTIC MEDIA . . .

PHYSICAL REVIEW E 95, 032417 (2017)

where we have defined the functions Aαβ and Bαβ in the last equation. The spatial derivatives of these functions with respect tor̄ = r/ arer̄α δβγ + r̄β δαγ − r̄γ δαβr̄α r̄β r̄γ
iωτ
∂
,
(B2)
Aαβ (r,ω) =
−3 5
r̄ 3
r̄
∂ r̄γ
1 + iωτr̄α r̄β r̄γr̄α δβγ + r̄β δαγ + r̄γ δαβ∂2Bαβ (r,ω) =+ 15 7
−3
.
(B3)
r̄ 5
r̄
∂ r̄γ
(1 + iωτ )2Using these results, we can further calculate the following quantities:∂∂(ωτ )2r̄ −4 ,Aαβ (r,ω)Aαβ (r,−ω) = 12
∂ r̄γ
∂ r̄γ
1 + (ωτ )2∂∂iωτ − (ωτ )2 −6r̄ ,15ββ γ γAαβ (r,ω)Bαβ (r,−ω) = −72
∂ r̄γ
∂ r̄γ
[1 + (ωτ )2 ]2∂∂−iωτ − (ωτ )2 −6r̄ ,15ββ γ γBαβ (r,ω)Aαβ (r,−ω) = −72
∂ r̄γ
∂ r̄γ
[1 + (ωτ )2 ]2∂∂720r̄ −8 .
15ββ γ γBαβ (r,ω)Bαβ (r,−ω) =
∂ r̄γ
∂ r̄γ
[1 + (ωτ )2 ]215ββ γ γ

Hence the dimensionless PSD in Eq. (33) can be obtained as∂gαβ [r,ω] ∂gαβ [r,−ω]d 3 r̄
I(ω) = 15ββ γ γ
∂ r̄γ
∂ r̄γ


 ∞
2
(ωτ )
(ωτ )2
720
−4
−6
−8
r̄
r̄
r̄
= 4π
d r̄ r̄ 2 12−144+
1 + (ωτ )2[1 + (ωτ )2 ]2[1 + (ωτ )2 ]2
δ/
 
 3
 5 (ωτ )2(ωτ )2
12
−4
= 48π
+
.
1 + (ωτ )2 δ
[1 + (ωτ )2 ]2 δ
[1 + (ωτ )2 ]2 δ

(B4)
(B5)
(B6)
(B7)

(B8)

Hence we finally obtain Eq. (34).

APPENDIX C: FULL EXPRESSION OF EQ. (42)The full expression of Eq. (42), including higher-order terms in /δ, is given as follows:

 

1
(τd /τ )e−t/τd
1 c0 Sτe−t/τ
2+−
(Rx ) (t) ≈22260π η 1 + τd /τ
1 − (τd /τ )1 − (τd /τ )δ 
[−2(t/τ )(τd /τ )2 + 2t/τ − 6(τd /τ )2 + 2]e−t/τ
4(τd /τ )3 e−t/τd2(2τd /τ + 1)  3++−[1 − (τd /τ )2 ]2[1 − (τd /τ )2 ]2[1 + τd /τ ]2δ
6[−(t/τ )(τd /τ )2 + t/τ − 5(τd /τ )2 + 3]e−t/τ12(τd /τ )5 e−t/τd++22[1 − (τd /τ ) ][1 − (τd /τ )2 ]2
 5 
6[2(τd /τ )3 + 4(τd /τ )2 + 6τd /τ + 3] 12t−+.
[1 + τd /τ ]2τδ[1] B. Alberts, A. Johnson, P. Walter, J. Lewis, and M. Raff,
Molecular Biology of the Cell (Garland Science, New York,
2008).
[2] C. P. Brangwynne, G. H. Koenderink, F. C. MacKintosh, and D.
A. Weitz, J. Cell Biol. 183, 583 (2008).
[3] F. C. MacKintosh, Proc. Natl. Acad. Sci. USA 109, 7138
(2012).
[4] S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Proc. Natl.
Acad. Sci. USA 109, 7338 (2012).

(C1)[5] B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O’Hern, E. R.
Dufresne, and C. Jacobs-Wagner, Cell 156, 183 (2014).
[6] 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).
[7] X.-L. Wu and A. Libchaber, Phys. Rev. Lett. 84, 3017 (2000).
[8] D. T. N. Chen, A. W. C. Lau, L. A. Hough, M. F. Islam, M.
Goulian, T. C. Lubensky, and A. G. Yodh, Phys. Rev. Lett. 99,
148302 (2007).

032417-13

KENTO YASUDA, RYUICHI OKAMOTO, AND SHIGEYUKI KOMURA[9] K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci, and R. E.
Goldstein, Phys. Rev. Lett. 103, 198103 (2009).
[10] Y. Hatwalne, S. Ramaswamy, M. Rao, and R. A. Simha,
Phys. Rev. Lett. 92, 118101 (2004).
[11] P. T. Underhill and M. D. Graham, Phys. Fluids 23, 121902(2011).
[12] F. C. MacKintosh and A. J. Levine, Phys. Rev. Lett. 100, 018104(2008).
[13] A. J. Levine and F. C. MacKintosh, J. Phys. Chem. B 113, 3820
(2009).
[14] A. S. Mikhailov and R. Kapral, Proc. Natl. Acad. Sci. USA 112,
E3639 (2015).
[15] R. Kapral and A. S. Mikhailov, Physica D 318, 100(2016).
[16] Y. Koyano, H. Kitahata, and A. S. Mikhailov, Phys. Rev. E 94,
022416 (2016).
[17] B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas,
and J. J. Fredberg, Phys. Rev. Lett. 87, 148102 (2001).
[18] P. Kollmannsberger and B. Fabry, Annu. Rev. Mater. Res. 41,
75 (2011).
[19] B. D. Hoffman, G. Massiera, K. M. Van Citters, and J. C.
Crocker, Proc. Natl. Acad. Sci. USA 103, 10259 (2006).
[20] B. D. Hoffman and J. C. Crocker, Annu. Rev. Biomed. Eng. 11,
259 (2009).
[21] X. Trepat, L. Deng, S. S. An, D. Navajas, D. J. Tschumperlin,
W. T. Gerthoffer, J. P. Butler, and J. J. Fredberg, Nature 447,
592 (2007).
[22] T. G. Mason and D. A. Weitz, Phys. Rev. Lett. 74, 1250
(1995).
[23] T. G. Mason, K. Ganesan, J. H. van Zanten, D. Wirtz, and S. C.
Kuo, Phys. Rev. Lett. 79, 3282 (1997).
[24] F. Gittes, B. Schnurr, P. D. Olmsted, F. C. MacKintosh, andC. F. Schmidt, Phys. Rev. Lett. 79, 3286 (1997).
[25] B. Schnurr, F. Gittes, F. C. MacKintosh, and C. F. Schmidt,
Macromolecules 30, 7781 (1997).
[26] T. G. Mason, Rheol. Acta 39, 371 (2000).
[27] T. M. Squires and T. G. Mason, Annu. Rev. Fluid Mech. 42, 413(2010).
[28] A. W. C. Lau, B. D. Hoffmann, A. Davies, J. C. Crocker, andT. C. Lubensky, Phys. Rev. Lett. 91, 198101 (2003).
[29] A. W. C. Lau and T. C. Lubensky, Phys. Rev. E 80, 011917(2009).
[30] P. G. de Gennes, Macromolecules 9, 587 (1976).
[31] P. G. de Gennes, Macromolecules 9, 594 (1976).

PHYSICAL REVIEW E 95, 032417 (2017)[32] F. Brochard and P. G. de Gennes, Macromolecules 10, 1157
(1977).
[33] J. C. Crocker, M. T. Valentine, E. R. Weeks, T. Gisler, P. D.
Kaplan, A. G. Yodh, and D. A. Weitz, Phys. Rev. Lett. 85, 888(2000).
[34] A. J. Levine and T. C. Lubensky, Phys. Rev. Lett. 85, 1774
(2000).
[35] A. J. Levine and T. C. Lubensky, Phys. Rev. E 63, 041510(2001).
[36] A. Sonn-Segev, A. Bernheim-Groswasser, H. Diamant, and
Y. Roichman, Phys. Rev. Lett. 112, 088301 (2014).
[37] A. Sonn-Segev, A. Bernheim-Groswasser, and Y. Roichman,
Soft Matter 10, 8324 (2014).
[38] K. Yasuda, R. Okamoto, S. Komura, and A. S. Mikhailov,
Europhys. Lett. 117, 38001 (2017).
[39] H. Diamant, Eur. Phys. J. E 38, 32 (2015).
[40] D. T. N. Chen, Q. Wen, P. A. Janmey, J. C. Crocker, and A. G.
Yodh, Annu. Rev. Condens. Matter Phys. 1, 301 (2010).
[41] C. Wilhelm, Phys. Rev. Lett. 101, 028101 (2008).
[42] F. Gallet, D. Arcizet, P. Bohec, and A. Richert, Soft Matter 5,2947 (2009).
[43] D. Robert, T.-H. Nguyen, F. Gallet, and C. Wilhelm, PLoS ONE
5, e10046 (2010).
[44] I. M. Tolić-Nørrelykke, E.-L. Munteanu, G. Thon, L.
Oddershede, and K. Berg-Sørensen, Phys. Rev. Lett. 93, 078102(2004).
[45] I. Golding and E. C. Cox, Phys. Rev. Lett. 96, 098102 (2006).
[46] S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Phys. Rev. Lett.
104, 238102 (2010).
[47] É. Fodor, M. Guo, N. S. Gov, P. Visco, D. A. Weitz, and F. van
Wijland, Europhys. Lett. 110, 48005 (2015).
[48] É. Fodor, W. W. Ahmed, M. Almonacid, M. Bussonnier, N.
S. Gov, M.-H. Verlhac, T. Betz, P. Visco, and F. van Wijland,
Europhys. Lett. 116, 30008 (2016).
[49] R. Bruinsma, A. Y. Grosberg, Y. Rabin, and A. Zidovska,
Biophys. J. 106, 1871 (2014).
[50] M. Doi and A. Onuki, J. Phys. II (France) 2, 1631 (1992).
[51] L. Cipelletti, S. Manley, R. C. Ball, and D. A. Weitz, Phys. Rev.
Lett. 84, 2275 (2000).
[52] A. Duri and L. Cipelletti, Europhys. Lett. 76, 972 (2006).
[53] Y. Gao, J. Kima, and M. E. Helgeson, Soft Matter 11, 6360
(2015).
[54] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics(Clarendon Press, Oxford, 1986).

032417-14