Soft Matter
View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

PAPER

Cite this: Soft Matter, 2020,
16, 10734

View Journal | View Issue

Mechanochemical enzymes and protein machinesas hydrodynamic force dipoles: the active
dimer model
Yuto Hosaka,a

Shigeyuki Komura*a and Alexander S. Mikhailovbc

Mechanochemically active enzymes change their shapes within every turnover cycle. Therefore, they induce
circulating flows in the solvent around them and behave as oscillating hydrodynamic force dipoles. Because of
Received 22nd June 2020,
Accepted 14th September 2020
DOI: 10.1039/d0sm01138j

non-equilibrium fluctuating flows collectively generated by the enzymes, mixing in the solution and diﬀusion ofpassive particles within it are expected to get enhanced. Here, we investigate the intensity and statisticalproperties of such force dipoles in the minimal active dimer model of a mechanochemical enzyme. In theframework of this model, novel estimates for hydrodynamic collective eﬀects in solution and in lipid bilayers

rsc.li/soft-matter-journal

under rapid rotational diﬀusion are derived, and available experimental and computational data is examined.

1 Introduction
Ligand-induced mechanochemical motions are typical for
enzymes. Binding or dissociation of a ligand (i.e., substrate or
product) to such proteins, as well as chemical reactions within theligand-bound state, are often accompanied by conformationaltransitions in them. Thus, these macromolecules would repeatedly
change their shapes in each next turnover cycle. The primary roleof mechanochemical motions is to enable and facilitate catalytic
reaction events. In the enzymes that operate as protein machines
or molecular motors and catalytically convert ATP or GTP, suchmotions are moreover employed to bring about the required
machine function or to generate work.
Since enzymes are in solution, their active conformationalchanges are accompanied by flows in the fluid around them.
Such non-equilibrium flows can aﬀect internal mechanicalmotions in the enzymes and also influence translational androtational diﬀusion of such proteins, as demonstrated by MD
simulations for a model protein1 and adenylate kinase.2 It has
been discussed whether hydrodynamic self-propulsion of enzymes
could furthermore occur, in the models where either instantaneous transitions3,4 or ligand-induced continuous conformationalmotions take place.5–7
Lipid bilayers, forming biological membranes, behave as twodimensional (2D) fluids on submicrometer scales.8,9 Biomembranes
often include many active protein inclusions, such as ion pumps ora

Department of Chemistry, Graduate School of Science, Tokyo Metropolitan
University, Tokyo 192-0397, Japan. E-mail: komura@tmu.ac.jpbComputational Molecular Biophysics, WPI Nano Life Science Institute, Kanazawa
University, Kakuma-machi, Kanazawa, 920-1192, Japanc
Department of Physical Chemistry, Fritz Haber Institute of the Max Planck Society,
Faradayweg 4-6, 14195 Berlin, Germany

10734 | Soft Matter, 2020, 16, 10734--10749

transporters. Essentially, they represent protein machines poweredby ATP hydrolysis or other catalytic reactions in them. Within each
operation cycle, the shapes of their membrane domains typicallychange, inducing 2D fluid flows in the lipid bilayer around them.10
As a result, active protein inclusions might even propel themselves
through biomembranes.11
Collective conformational activity of enzymes and protein
machines leads to the development of non-thermal fluctuatingflows in solution or a lipid bilayer. Other particles (i.e., passive
tracers) are advected by these non-equilibrium flows, and, as
previously shown,12 increased mixing in such systems anddiﬀusion enhancement should therefore arise. Additionally,
chemotaxis-like eﬀects in the presence of spatial gradients inthe concentration or the activity of enzymes can take place.12
Remarkably, such phenomena persist even if mechanochemical
motions are reciprocal; they do not rely on the presence of selfpropulsion for proteins, which is predicted to be weak.5–7
Following the original publication,12 extensive further researchhas been performed.13–20 The effects of rotational diffusion and of
possible nematic ordering for enzymes were considered,14 thephenomena in biomembranes were extensively analyzed,15,16 andthe theory was extended to viscoelastic media as well.17,18 Recently,
it was shown that viscosity in dilute solutions of mechanochemically active enzymes should become also reduced.19 Multiparticle numerical simulations of active oscillatory colloids, explicitly
including hydrodynamic effects, were furthermore undertaken andprincipal theoretical predictions could thus be verified.20
At low Reynolds numbers, the flow distribution produced byan object, changing the shape due to internal forces within it,
can be characterized in the far field as that corresponding to ahydrodynamical force dipole. If the time-dependent stochastic
force dipole of an enzyme is known, the collective hydrodynamic

This journal is © The Royal Society of Chemistry 2020

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Paper

Soft Mattereﬀects in solution of such enzymes are predicted by the mean-field
theory.12 The diﬃculty, however, is that experimental measurements and precise theoretical estimates for intensities and statistical properties of the force dipoles corresponding to actual
enzymes are not available yet. Lacking this knowledge, only
rough quantitative estimates for the considered collective hydrodynamic eﬀects could be made so far.
Our present study has two aims and, respectively, it includes
two parts. Section 2 corresponds to the first part. Here, theactive dimer model is formulated. The active dimer represents aminimal model where ligand-induced mechanochemical motions
are reproduced.12,20–22 After presenting the model, we undertake
an approximate analytical investigation of statistical properties ofthe force dipoles corresponding to active dimers in Section 2.2,
followed by a numerical study in Section 2.3. Quantitative
estimates for the intensity of hydrodynamical force dipoles inreal enzymes are obtained in Section 2.4.
Section 3 corresponds to the second part. Based on theactive dimer results, we obtain in Section 3.1 more precise
analytical and numerical estimates for the maximal diﬀusion
enhancement for passive particles in solutions of active enzymes,
taking into account fast rotational diﬀusion of enzymes. Similar
estimates for diﬀusion enhancement of passive particles in lipid
bilayers are derived in Section 3.2.
The results are discussed in Section 4. There, we analyze theavailable experimental and computational data for diﬀusion
enhancement in, respectively, Sections 4.1 and 4.2. Conclusions
and an outline for the perspectives of further research are provided
in Section 5.

2 Statistical properties of force dipoles
2.1

The active dimer model

The simplest mechanical system that gives rise to a hydrodynamical force dipole is a dimer. It consists of two beads 1 and2 interacting via a potential u(r) that depends on the distancer = |r1  r2| between them. The forces acting on the particles are
f1 = qu/qr1 = f and f2 = f. If the dimer is immersed into aviscous fluid, the velocity V of the hydrodynamic flow far
enough from the dimer is approximately given by12
Va ¼

@Gab
eb eg m;@Rg

(1)where Gab(R) is the mobility tensor depending on the position Rof the dimer with respect to the observation point, e = (r1  r2)/ris the unit orientation vector of the dimer, and m = fr is themagnitude of the force dipole. Summation over repeatedindices is assumed. The force dipole is present only if there
are non-vanishing net interaction forces, i.e., if the distance
between the particles in a dimer continues to change. As in the
study,12 we assume that the Oseen approximation holds. For a
dimer, it is justified if the distance between the beads is muchlarger than their size.
The minimal active dimer model has been proposed12,21 (seealso review ref. 23) to imitate mechanochemical conformationalmotions accompanying a catalytic turnover cycle in an enzyme. Note

This journal is © The Royal Society of Chemistry 2020

Fig. 1 The turnover cycle and mechanochemical motions in the active
dimer model of an enzyme (see the text).

that the dimer model, with non-reactive dissociation of substrate
additionally included, was also considered in the study.22
The operation mechanism is illustrated in Fig. 1. Two
identical beads (green) of radius a are connected by an elastic
link with a certain natural spring length c0 and stiﬀness k0. A
substrate particle (red) arrives (A) and binds as a ligand to thedimer by forming an additional elastic link with stiﬀness k thatconnects the two beads (B). The natural length cc of thisadditional link is taken to be shorter than c0. Therefore, it tends
to contract the dimer until a new equilibrium conformation (C) witha certain distance c1 between the beads is reached. Once this hastaken place, a chemical reaction, that converts the ligand from thesubstrate to the product, occurs and the product (blue) is instantaneously released (D). Following the product release, the dimer is in
the state E with the spring length c1 that is shorter than the naturallength c0. Therefore, the spring expands and the domains moveapart until the equilibrium state (F) is approached again. After that,a new substrate can bind, repeating the turnover cycle.
It is assumed that products are immediately evacuated andtherefore we do not consider reverse product binding events.
Moreover, possible dissociation events for the substrate areneglected assuming that its aﬃnity is high. Note that, since theproduct is immediately released once it has been formed, theligand inside our model enzyme is always only in the substrateform. Therefore, the dimer can be either in the ligand-free(s = 0) or the ligand-bound (s = 1) states.
The elastic energies in these two states are
E0 ðxÞ ¼

k0
ðx  ‘0 Þ2 ;
2

(2)

and
E1 ðxÞ ¼

k0
kk1
ðx  ‘ 0 Þ2 þ ðx  ‘ c Þ2 ¼ A þ ðx  ‘ 1 Þ2 ;
2
2
2

(3)

where x is the distance between the beads andA¼

kk0ð‘0  ‘c Þ2 ;2ðk0 þ kÞ

k1 ¼ k0 þ k;(4)

k0 ‘0 þ k‘c‘1 ¼:
k0 þ kSoft Matter, 2020, 16, 10734--10749 | 10735

View Article Online

Soft Matter

Paper

The overdamped dynamics of the dimer in the ligand state s
is described by the Langevin equationdx@Esþ xðtÞ;
¼ g
@x
dt

(5)

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

where g is the mobility coeﬃcient. To account for thermal
fluctuations, this equation includes thermal noise,hxðt1Þxðt2Þi = 2gkBTdðt1  t2Þ,(6)

where kB is the Boltzmann constant and T is the temperature.
In eqn (5), we have omitted hydrodynamic interactions between
the beads. They were taken into account in the study22 of diﬀusion
enhancement for a single dimer itself. In the Oseen approximation, such interaction terms are proportional to the small
parameter a/c0, leading to corrections of the same order for theforce dipoles, neglected by us.
Stochastic transitions between the two ligand states takeplace at constant rates v0 and v1 within narrow windows of
width r near x = c0 and x = c1. If probability distributions ps(x,t)
are introduced, they obey a system of two coupled Fokker–
Planck equations@p0@@ 2 p0¼½gk0 ðx  ‘0 Þp0  þ gkB T 2
@x
@t
@x

(7)

þ u1 ðxÞp1 ðxÞ  u0 ðxÞp0 ðxÞ;
and
@p1
@
@ 2 p1½gk1 ðx  ‘1 Þp1  þ gkB T 2
¼
@x
@t
@x

(8)

þ u0 ðxÞp0 ðxÞ  u1 ðxÞp1 ðxÞ;where u0(x) = v0 for c0  r o x o c0 + r and vanishes outside ofthis interval; u1(x) = v1 for c1  r o x o c1 + r and zero outside
the interval. Note that the rate v0 of substrate binding is
proportional to the substrate concentration.
If the transition windows are very narrow, i.e., r { c0 andr { c1, one can use the approximationu0ðxÞ = n0dðx  c0Þ,

u1ðxÞ = n1dðx  c1Þ,

(9)where n0 = 2v0r and n1 = 2v1r.
Fig. 2 shows the energy diagram of the model. Within eachcycle, the dimer dissipates in mechanochemical motions theenergy DE = DE0 + DE1 which is furthermore equal to thediﬀerence Esub  Eprod of the energy Esub = E1(c0)  E0(c0)supplied with the substrate and the energy Eprod = E1(c1) E0(c1) removed with the product. We have
DE ¼

1ðk0 þ k1 Þð‘0  ‘1 Þ2 :2

(10)

The energy diﬀerence DE is always positive and, hence, theconsidered active dimer represents an exothermic enzyme.
The force dipole of the active dimer is m = k0(c0  x)x for s =0 and m = k1(c1  x)x for s = 1. Note that therefore m r k0c02/4
for s = 0 and m r k1c12/4 for s = 1.

10736 | Soft Matter, 2020, 16, 10734--10749

Fig. 2The energy diagram of the active dimer.

When the transition windows are narrow, the probability
rate w0 that substrate binding, i.e., a transition to state s = 1,
occurs per unit time in the state s = 0 is approximatelysﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
k0
w0 ¼ n 0:(11)
2pkB T
On the other hand, the probability rate w1 that product release,i.e., a transition to state s = 0, occurs per unit time in the states = 1 is then approximately given bysﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
k1
w1 ¼ n 1:(12)
2pkB T
These equations are derived in Appendix A. Moreover, thecharacteristic relaxation times of the dimer in the states s = 0
and s = 1 are, respectively, t0 = (gk0)1 and t1 = (gk1)1.
The parameter combinations w0t0 and w1t1 play an important role in determining the kinetic regimes. If the condition
w0t0 { 1 is satisfied, equilibration to thermal distribution in
the state s = 0 usually takes place before a transition to the states = 1, i.e., binding of a substrate, occurs. If the opposite
condition w0t0 c 1 holds, such transition takes place immediately
after the transition window at x = c0 is reached. If w1t1 { 1, theequilibration takes place in the state s = 1 before a transition to thestate s = 0, i.e., the reaction and the product release, occurs. In
the opposite limit with w1t1 c 1, the reaction takes place andproduct becomes released immediately once the respective
window at x = c1 is reached.
Note that, because the rate w0 is proportional to substrateconcentration, the condition w0t0 c 1 corresponds to thesubstrate saturation regime for the considered model enzyme.
The condition w1t1 { 1 implies that the enzyme waits a longtime before the product is released.
2.2

Approximate analytical results for force dipoles

At thermal equilibrium in the absence of substrate, p1(x) = 0 andsﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ


k0
k0
2
ðx  ‘ 0 Þ :
exp 
(13)
p0 ðxÞ ¼
2pkB T
2kB T

This journal is © The Royal Society of Chemistry 2020

View Article Online

Paper

Soft Matter

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Since m = k0(c0  x)x, one can easily find the equilibrium
statistical distribution for force dipoles by using the condition
Peq(m)dm = p0(x)dx. Using, for convenience, the dimensionlessforce dipole magnitude m̃ = m/(k0c02) and dimensionless temperature y = kBT/(k0c02), we get 
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ2
1
1
~ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ exp ~1 þ 1  4m
Peq ðmÞ
8y~2pð1  4mÞy
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ21
~
1  1  4m
þ exp 
:
8y

(14)

(15)

Using the distribution in eqn (14), one finds that the mean
force dipole ishmieq = kBT.

(16)

The correlation function C(t) = hDm(t)Dm(0)i for variations
Dm = m  hmi of force dipoles is20CeqðtÞ = k0c02kBTe|t|/t0 + 2ðkBTÞ2e2|t|/t0,

(17)

1

where t0 = (gk0) is the characteristic relaxation time for the
dimer in the state s = 0. As shown in Appendix B, the exact
relation hmi = kBT holds for the dimer in any steady state.
For an active dimer, approximate analytical estimates can beobtained in the four characteristic limits described below. The twoof them (A and C) correspond to low substrate concentrations,with rare turnover cycles controlled by the substrate supply. In
regime B, mechanochemical motions are limiting the overall
catalytic rate. In other words, product formation and its releaseoccur once an appropriate conformation (x = c1) has been reached.
In regime D, the overall kinetic rate is, on the other hand, limitedby the waiting time for product formation and release.
A. The limit of w0s0 { 1 and w1s1 { 1. If these conditions
are satisfied, binding of the substrate and product release havelarge waiting times. In this limit, there are two almost independent equilibrium subpopulations of dimers in the states s = 0 ands = 1. The relative weights of the subpopulations are w1/(w1 + w0)
and w0/(w1 + w0). Therefore, all statistical properties are given by
the sums of contributions from diﬀerent states taken with therespective weights. Particularly, the correlation function of force
dipoles isCðtÞ ¼i
w1 h
k0 ‘02 kB Tejtj=t0 þ 2ðkB T Þ2 e2jtj=t0
w0 þ w1i (18)w0 h
þ
k1 ‘12 kB Tejtj=t1 þ 2ðkB T Þ2 e2jtj=t1 :w0 þ w1

We can use the above equation to determine the nonequilibrium part of the fluctuation intensity of force dipoles
hDm2iA = hDm2i  hDm2ieq.

This journal is © The Royal Society of Chemistry 2020

Dm2 A ¼


w0
k1 ‘12  k0 ‘02 kB T:
w0 þ w1

(19)

(20)

As follows from eqn (17), the equilibrium fluctuation intensity is
hDm2ieq = k0c02kBT + 2ðkBTÞ2.

If y { 1, this distribution is approximately Gaussian andlocalized at m = 0, i.e.,~2m1~ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃ exp 
Peq ðmÞ
:2y2py

Because hDm2i = C(0), we have(21)

Since the eﬀective binding rate w0 of the substrate is proportional
to its concentration c, i.e., w0 = Zc, eqn (20) yields the Michaelis–Menten form of the dependence of hDm2iA on the substrateconcentration.
Remarkably, the catalytic activity of the model enzyme canthus lead not only to some enhancement, but also to reduction offluctuations of the force dipoles. According to eqn (20), reduction
should be observed if k1c12 o k0c02. Under this condition, theligand-bound dimer (s = 1) is characterized by a lower fluctuation
intensity of force dipoles than the free dimer (s = 0).
B. The limit of w0s0 c 1 and w1s1 c 1. In this limit,transitions take place once the respective transitions windows
are entered. If additionally the conditions k0c02 c kBT and k1c12 c
kBT are satisfied, thermal fluctuations can be neglected and thedimer essentially behaves as a deterministic oscillator. Then, thesolution can be obtained by integrating eqn (5) with appropriate
boundary conditions. This yields x(t) = c1 + (c0  c1  r)et/t1 for
0 o t o T1 and x(t) = c0 + (c1  c0 + r)e(tT1)/t0 for T1 o t o Tc.
Here, Tc is the oscillation period of the active dimer and T1 is theduration of the cycle time when the dimer is in the ligand-boundstate s = 1. If transition windows are narrow, i.e., the conditionr { (c0  c1) is satisfied, we approximately have
T1 ¼ t1 ln‘0  ‘1
;
r

(22)

and
Tc ¼ ðt0 þ t1 Þ ln

‘0  ‘1:r(23)

The respective time-dependent force dipole is m(t) = k1(c1  x)x
for 0 o t o T1 and m(t) = k0(c0  x)x for T1 o t o Tc. Hence, it isnegative for s = 1 and positive for s = 0.
The force dipole varies within the interval mmin o m ommax, where the minimum value mmin = k1c0(c0  c1) is taken
at t = 0, i.e., in the state s = 1 just after substrate binding, andthe maximum value mmax = k0c1(c0  c1) is reached at t = T1, inthe state s = 0 just after product release (here we again assume
that transition windows are narrow). Note that, if thermal
fluctuations were present, the force dipoles could however havealso taken the values outside of this interval.
It can be checked by direct integration that the periodaveraged force dipole for the deterministic active dimer ishm(t)idet = 0. The correlation function for the deterministic
oscillating dimer is defined as the period averageð
1 Tc
Cdet ðtÞ ¼
dhmðt þ hÞmðhÞ:
(24)
Tc 0
The explicit analytical form of this periodic correlation function
is too complicated and we do not give it (analytical results forSoft Matter, 2020, 16, 10734--10749 | 10737

View Article Online

Soft Matter

Paper

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

correlation functions are also omitted below in the limits
C and D).
The mean-square intensity of force dipoles is hm(t)2idet =Cdet(0). In the limit r - 0, we approximately have


k0 k1
‘0  ‘1 1ln
mðtÞ2 det ¼ð‘0  ‘1 Þ2
12ðk0 þ k1 Þr
 k0 ‘02 þ 2‘0 ‘1 þ 3‘12 þ k1 3‘02 þ 2‘0 ‘1 þ ‘12 :(25)
When k1 B k0, this equation yields the scaling hm(t)2idet B k02.
C. The limit of w0s0 { 1 and w1s1 c 1. If these conditions
are satisfied, the model enzyme waits a long time for binding ofa substrate (because the substrate concentration is low), butthen it performs a rapid reaction cycle. An approximate
solution in this regime can be obtained if, additionally, theconditions k0c02 c kBT and k1c12 c kBT are satisfied, i.e., that
thermal fluctuations are weak. Moreover, we shall assume thatthe transition window for substrate binding is narrow, i.e., theapproximation in eqn (9) holds for u0(x).
In this case, the dependence x(t) consists of a sum ofstatistically independent rare pulses, each corresponding to
one reaction cycle:
X(26)xðtÞ ¼z t  tj ;jwhere z(t) = c1 + (c0  c1)efor 0 o t o T1 and z(t) = c0 +(c1  c0)e(tT1)/t0 for t 4 T1, with T1 given by eqn (22). The
pulses appear at random time moments tj and the probabilityof their appearance per unit time is w0.
Moreover, we also haveX(27)z t  tj ;mðtÞ ¼
t/t1

jwhere z(t) = k1(c1  z(t))z(t) for 0 o t o T1 and z(t) =k0(c0  z(t))z(t) for t 4 T1.
Hence, this represents a random Poisson process. Its firsttwo statistical moments are approximately hm(t)i = 0 andð11
m2 ðtÞ ¼ w0 dtz2 ðtÞ ¼ w0 t0 ð‘0  ‘1 Þ2
12
0
 2 2


 k0 ‘0 þ 2‘0 ‘1 þ 3‘12 þ k0 k1 3‘02 þ 2‘0 ‘1 þ ‘12 :If we take into account eqn (12), it can be noticed that, when
k1 B k0, scaling hm2(t)i B k3/2
0 is again obtained.
2.3

Numerical simulations

Numerical simulations can yield statistical properties of force
dipoles for selected parameter values in the regions where there
are no approximate analytical results. Below in this section,
we focus on the situation under substrate saturation, but withthe waiting time for product release comparable to the conformational relaxation time (thus, lying between the limits B
and D). We shall consider a situation where the condition
kBT { k0c02 is satisfied, so that thermal fluctuations are
relatively weak.
Before proceeding to simulations, the model was nondimensionalized. The dimensionless variables were t̃ = t/t0,x̃ = x/c0, and m̃ = m/(k0c02). The dimensionless transition rates
were ṽ0 = v0t0 and ṽ1 = v1t0, while the dimensionless temperature
was y = kBT/(k0c02). Stochastic diﬀerential eqn (5) was numericallyintegrated, complemented by transitions between the ligand states.
In the simulations, we had c1 = 0.55c0, k1 = 2k0, and r =0.01c0. We have kept constant ṽ1 = 2, but varied the parameterṽ0. A relatively low dimensionless temperature y = 0.0018 waschosen to satisfy the condition kBT { k0c02. Under such choice,hm(t)2idet/hDm2ieq = 19.3 and w1t1 = 0.27.
Note that, because of the last condition, there was a significant
random variation in the waiting times for substrate conversion andproduct release. Moreover, waiting times for substrate binding,characterized by the rate w0, could also vary. These eﬀects kept
the model stochastic even when thermal noise was small.
Fig. 3 shows typical time dependences of the force dipoles. In
Fig. 3(a), the waiting time for substrate binding is long. Therefore, the dimer spends most of the time in the ligand-free state s= 0. Within the time shown, only one turnover cycle has takenplace. For the force dipole, the cycle consists of a negative
spike, just after binding of the substrate, and the following
positive spike, just after the product release. In Fig. 3(b), thesubstrate binding rate is increased. As a result, the dimer is(28)
Taking into account eqn (11), we notice that, when k1 B k0, thescaling hm2(t)i B k03/2 should hold.
D. The limit of w0s0 c 1 and w1s1 { 1. This situation
corresponds to substrate saturation and a long waiting time forthe reaction and product release in the ligand-bound state.
A derivation, similar to that given above, shows that, if k0c02 c
kBT and k1c12 c kBT, we approximately have hm(t)i = 0 andm2 ðtÞ ¼

1
w1 t1 ð‘0  ‘1 Þ2
12



 k12 3‘02 þ 2‘0 ‘1 þ ‘12 þ k0 k1 ‘02 þ 2‘0 ‘1 þ 3‘12 :(29)

10738 | Soft Matter, 2020, 16, 10734--10749

Fig. 3 Time dependence of dimensionless force dipoles m̃ = m/(k0c02) ontime for (a) ṽ0 = 0.03 and (b) ṽ0 = 3. Dashed lines show the lower boundm̃min = 0.9 for the deterministic oscillatory dimer and the absolute upperbound m̃max = 0.25 for force dipoles.

This journal is © The Royal Society of Chemistry 2020

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Paper

Fig. 4 Probability distributions of force dipoles m̃ for passive (black curve,v0 = 0) and active (red curve, ṽ0 = 0.03, and blue curve, ṽ0 = 3) dimers.

Soft Matter

Fig. 6 Normalized correlation functions of force dipoles at diﬀerent
substrate binding rates: v0 = 0 (absence of substrate, black), ṽ0 = 0.03
(red), and ṽ0 = 3 (blue). The correlation function for passive dimers (black)is given by eqn (17). Dashed curves are fits to the dependence in eqn (30).

frequently cycling, already resembling an oscillator. Nonetheless, therandom variation of the times between the cycles is relatively large.
Probability distributions of force dipoles are shown in Fig. 4.
The black curve is the distribution for passive dimers in theabsence of the substrate, given by eqn (14). It represents a
narrow Gaussian peak at m = 0. The distribution at ṽ0 = v0t0 =0.03 (red) is almost indistinguishable from it. The blue curve isthe distribution for active dimers corresponding to Fig. 3(b).
Now, the distribution is more broad and the central peak issmaller. The tail on the left side from the peak and the shoulder
on its right side are due to the non-equilibrium activity of force
dipoles.
The dependence of the non-equilibrium part of the fluctuation
intensity of force dipoles, eqn (19), on the substrate binding rate
v0, proportional to substrate concentration, is shown in Fig. 5. It
can be well fitted to the Michaelis–Menten function (the solid
curve). The saturation magnitude is close to the value of 0.033
predicted at such parameters for the deterministic dimer byeqn (25).
Normalized correlation functions of force dipoles at diﬀerent
substrate binding rates are shown in Fig. 6. In the absenceof the substrate (for v0 = 0) the dependence is monotonousFig. 7 shows how the dimensionless relaxation time 1/(Gt0), thedimensionless oscillation period 2p/(Ot0) and the phase shift adepend on the substrate binding rate. The oscillation period
under saturation conditions is still larger than Tc/t0 = 5.7 forthe deterministic dimer according to eqn (23). This is becauseof an additional waiting time for product release. The characteristic relaxation time is about 1/(Gt0) = 2.
It should be stressed that the form in eqn (30) of thecorrelation function would not hold in the deterministic limit.
Indeed, the oscillations stay harmonic in the limit of an infinite
correlation time. However, the deterministic oscillations are
actually non-harmonic, as seen in Fig. 3.

Fig. 5 Dependence of the non-equilibrium part hDm2iA of the fluctuation
intensity of force dipoles on the substrate binding rate v0 (dots). The solid
curve is a fit to the Michaelis–Menten function.

Fig. 7 The dependences of the relaxation time 1/G (circles), oscillation period
2p/O (triangles) and phase shift a (squares) on substrate binding rate v0.

This journal is © The Royal Society of Chemistry 2020

(it is given by eqn (17)). As the substrate concentration isincreased, damped oscillations in the correlation function
become observed, thus signaling the onset of the active oscillatory
behavior that prevails over the thermal noise.
The correlation functions could be fitted (dashed curves inFig. 6) to the dependenceCðtÞ=Cð0Þ ¼

1
expðGjtjÞ cosðOjtj  aÞ:
cos a

(30)

Soft Matter, 2020, 16, 10734--10749 | 10739

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Soft Matter
There are two eﬀects that make the dimer model stochastic,i.e., the thermal noise in the dynamical eqn (5) and randomtransitions between the ligand states s = 0 and s = 1. Wheny - 0, the thermal noise vanishes, but random transitions
between the states nonetheless remain. This second stochastic
eﬀect is responsible for the decay in the correlation function. As
shown in Appendix C, the dependence of the correlation
function in eqn (30) corresponds to an approximate solutionof the master eqn (7) and (8).
2.4

Estimates for hydrodynamic force dipoles of enzymes

Above, statistical properties of force dipoles were analyzed in
the framework of an idealized model of the active dimer. Now,
the obtained results can be applied to approximately estimate
the force dipoles for real enzymes and protein machines. To dothis, the relationship between such a simple model and theactual proteins needs to be first discussed.
Proteins fold into a definite conformation that however
incorporates many diﬀerent substates. Slow dynamics of proteins represents wandering over a Markov network of suchmetastable conformational substates.24 In all-atom molecular
dynamics (MD) simulations, transitions within tens of nanoseconds to the nearest metastable states can be clearly seen.
Long MD simulations show motions over a set of these states
extending to the millisecond timescales.25 Single-molecule
fluorescence correlation spectroscopy experiments with cholesteroloxidase revealed that thermal conformational fluctuations in thisenzyme, in the absence of the substrate, had correlations persistingeven over about 1.5 s time.26,27 In the coarse-grained structure-based
simulations of proteins, such as modeling based on elastic networks, the rugged atomic energy landscape becomes smoothed,28
thus yielding continuous slow conformational dynamics describedby a set of eﬀective collective coordinates.
In a detailed study of adenylate kinase,29 combining all-atom
MD simulations with single-molecule fluorescence resonanceenergy transfer and NMR, it was found that, in this characteristic
mechanochemical enzyme, conformational substates lie along a
trajectory that connects the initial open apo conformation to thefinal catalytically eﬃcient closed state. Thus, the energy landscape has a valley that guides towards the optimal protein state;
the motion along such a valley can be described by a single
coordinate. Similar organization of the energy landscape has
been noticed in structure-based coarse-grained modeling ofprotein machines and molecular motors, such as myosin V
and F1-ATPase30 and HCV helicase.31
Typically, mechanochemical enzymes and molecular machines
represent proteins with domain structure. Slow functional conformational dynamics in these proteins consists in relative
motions of the domains that can be often characterized by a single
coordinate, such as a hinge angle or a distance between the centersof mass of two protein domains. This leads to the reduced models
for proteins, with just one or a few mechanical coordinates.23 The
active dimer is a model belonging to this class. Note that previously a similar simple model with three beads was employed toestimate the magnitude of self-propulsion eﬀects in the enzymes.7
In the framework of the active dimer model, statistical properties

10740 | Soft Matter, 2020, 16, 10734--10749

Paperof force dipoles in diﬀerent kinetic regimes can be analyzed andcharacteristic order-of-magnitude estimates for the intensity ofsuch dipoles for typical enzymes and protein machines can bederived.
In Section 2.2, four kinetic regimes have been outlined. The twoof them (A and C) correspond to low substrate concentrations,with rare turnover cycles controlled by the substrate supply. Below,
we focus on the substrate saturation regimes B and D where highcatalytic activity and, thus, the strongest non-equilibrium force
dipoles can be expected.
In regime B, mechanochemical motions are limiting theoverall catalytic rate. In other words, product formation and its
release occur once an appropriate conformation (x = c1) has
been reached. Such regime is characteristic, for example, foradenylate kinase where the turnover time is limited by the time(about 1 ms) of the conformational transition from the open to
the closed state, with the reaction AMP + ATP - 2ADP rapidly
occurring once the latter state is reached.29
In regime D, the overall kinetic rate is, on the other hand,
limited by the waiting time for product formation and release.
This regime is typical for protein machines and motors such as
myosin V. In each operation cycle of this molecular motor, catalytichydrolysis of substrate ATP into product ADP takes place. The cycle
duration of 66 ms under ATP saturation is limited by waiting for
ADP release.32 The conformational transition from the open to theclosed state, i.e., the lever-arm swing after ATP binding, takes place
within a much shorter millisecond time.
The principal parameters of the active dimer model are
stiﬀness constants k0 and k1 and inter-domain distances c0
and c1 in the open (s = 0) and closed (s = 1) conformations,
respectively. The typical size of a protein is of the order of tens ofnanometers and this would be also the characteristic distances
c0 and c1 between the domains. Moreover, if the open andclosed states are distinctly diﬀerent, as, for example, in adenylate kinase or myosin, the change Dc = c0  c1 is comparable inmagnitude to c0 and c1. As characteristic values for order-ofmagnitude estimates, one can, for example, choose c0 = 10 nm
and c1 = 5 nm in the open and the closed states, respectively.
In the active dimer, two domains (beads) are connected by aspring. In real proteins, they can be, instead, connected by ahinge with the elastic energy
1
E ¼ K ðY  Y0 Þ2 ;2(31)which depends on the deviation of the hinge angle Y from theequilibrium angle Y0. This can also be approximately written asE¼

1
kðx  ‘0 Þ2 ;2

(32)

so that the hinge is described as an elastic spring with x = LY,
c0 = LY0 and the eﬀective stiﬀness k = K/L2, where L is thecharacteristic linear size of the domains connected by the hinge.
The stiﬀness of the converter hinge in myosin V was estimatedin single-molecule experiments by Kinoshita with coworkers33 to
be about K = 5 kBT rad2 both in the open and the closed states.
On the other hand, the data of high-speed AFM observations

This journal is © The Royal Society of Chemistry 2020

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Paper

Soft Matterby Ando with coworkers34 corresponds to a higher value of
K = 23 kBT rad2. The diﬀerence may be due to the fact thatthe hinge becomes softer for larger angles. In our estimates
below, we take K = 10 kBT rad2. Choosing L = 10 nm, this leads
to k = 0.1 kBT nm2.
As noticed above, in adenylate kinase, the overall turnoverrate under substrate saturation in an enzyme is limited by
conformational transitions between the open and closed states(and hence the turnover rate is about 103 s1). The maximum
intensity of force dipoles can be estimated by using eqn (25).
If the parameter values k0 = k1 = 0.1 kBT nm2, c0 = 10 nm,
c1 = 5 nm and r = 1 nm are chosen, the non-equilibrium meansquare fluctuation intensity hDm2iA of force dipoles in such
enzymes is estimated as approximately 80 [pN nm]2. This issimilar to the previous estimate of 100 [pN nm]2 in ref. 12 basedon typical stall forces in molecular motors.
In more slow enzymes and protein machines with the turnover
numbers of tens per second, the turnover is limited by productformation and its release. In this case, the intensity of force dipoles
can be estimated using eqn (29). There, the rate w1 of productformation and release is approximately the same as the overallturnover rate, whereas t1 corresponds to the conformationaltransition time. Choosing the turnover rate of 15 s1, as in myosin
V, and the conformational transition time of 1 ms and keeping thesame other parameters as above, the non-equilibrium meansquare fluctuation intensity of force dipoles can then be estimatedas about 4 [pN nm]2. This is much smaller than the above estimate
for fast enzymes because non-thermal mechanical forces are only
generated in conformational transitions of about 1 ms in duration.
It represents only a small fraction of the entire cycle time of tens of
milliseconds in such enzymes or protein machines.
While typical enzymes have turnover times between milliseconds
and tens of milliseconds, there are also very slow enzymes, such as
tryptophan synthase with the turnover time of 0.5 s,35 and enzymes
that are very fast, such as catalase (17 ms) or urease (59 ms).36
Moreover, transition times from open to close confirmations can bealso very short in some enzymes. For example, for phosphoglycerate
kinase (PGK), neutron spin-echo spectroscopy yields conformationaltransition times of the order of tens of nanoseconds.37 This was also
observed in coarse-grained MD simulations for PGK.38 Therefore, it
is interesting to discuss under what general conditions stronger
force dipoles can be expected in enzymes.
Eqn (10) relates the energy (generally, enthalpy) dissipated inmechanochemical motions within the turnover cycle of an active
dimer to the stiﬀness of the dimer and the magnitude of conformational changes in it. While it has been derived for an idealized
model, it can also be used for order-of-magnitude estimates in real
enzymes. Taking, for example, k0 = k1 = 0.1 kBT nm2 and c0  c1 =10 nm for myosin, we obtain DE = 10 kBT, which is in reasonableagreement with the energy of about 20 kBT supplied to this molecularmotor with ATP (only half of this energy is used in the power stroke).
Note that, assuming for simplicity that k0 = k1 = k, eqn (10)can be also written asDE
k ¼ 2;
D‘

This journal is © The Royal Society of Chemistry 2020(33)

thus expressing the stiﬀness in terms of the energy DE dissipated in mechanochemical motions and the conformationalchange Dc = c0  c1. An enzyme is stiﬀer if the same energy is
dissipated within a conformational transition of a smaller
magnitude.
Suppose that conformational changes are indeed small in an
enzyme and, moreover, its turnover rate is limited by conformationaltransitions within the cycle. Then, eqn (25) can be used to estimate
the intensity of force dipoles. For approximate numerical estimates,
it can be written in the formhDm2i = z0k2c02ðc0  c1Þ2,(34)where z0 is a dimensionless factor of order unity that also
includes the logarithmic term and we have taken k0 = k1 = k.
Note that this estimate holds assuming that the force dipoles in
the catalytically active enzyme are much stronger than thosedue to thermal fluctuations in the absence of substrate.
Substituting k from eqn (33), a simple order-of-magnitude
estimate is obtained
Dm2 ¼ z1

‘0D‘

2

DE 2 ;(35)

where z1 is another dimensionless factor of order unity. Moreover, by using eqn (21) and (33), we furthermore get
Dm2
DE
¼ z1
;
kB T
hDm2 ieq(36)

if the condition kc02 c kBT holds.
These results show that the intensity of force dipoles isstrongly sensitive to the magnitude of mechanochemical motions
within the turnover cycle. Moreover, they show that, in strongly
exothermic enzymes, force dipoles are greatly enhanced whencatalytic activity takes place.
The above-mentioned catalase and urease enzymes are not
only exceptionally fast, but also highly exothermic, with DH =100 kJ mol1 for catalase and DH = 59.6 kJ mol1 for urease.36
Hence, large energies of 42 kBT or 25 kBT are released in themand dissipated into heat within very short microsecond cycle times.
Furthermore, at least for catalase, it is known that functional
conformational changes are involved within the turnover cycle,
but their magnitude is small.39 It has been previously proposed36
that chemoacoustic intramolecular eﬀects caused by strong heatrelease may even lead to hydrodynamic self-propulsion of these
enzymes, although subsequent examination could not confirm
this.7 These enzymes do not have a domain structure and therefore
the results of our analysis based on the dimer model are not directly
applicable to them. Nonetheless, they suggest that hydrodynamic
force dipoles in them may be very strong.

3 Diﬀusion enhancement for passive
particles in solutions in active enzymes
The most important application of the obtained results for force
dipoles is that they allow to obtain more accurate analytical andnumerical estimates for diﬀusion enhancement of passive

Soft Matter, 2020, 16, 10734--10749 | 10741

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Soft Matter

Paper

particle in solutions of active enzymes. In the previous studies,12,14
the magnitude of diﬀusion enhancement has been expressed interms of the statistical properties of hydrodynamic force dipoles.
However, because these properties were only poorly known, precise
estimates for such magnitude could not be obtained. This has led
to diﬃculties in the interpretation of experimental results and in
the analysis of the computational data. In this section, we use thestatistical properties of force dipoles for active dimers, determined
above in Section 2, to estimate the diﬀusion enhancement eﬀects
for enzymes in water solution and for active protein inclusions inbiomembranes.
3.1

12,14

As previously shown,the change DA in the diﬀusion coeﬃcientof passive tracer particles in a three-dimensional (3D) solution is
given byn
DA ¼
w  weq ;(37)60pm2 ‘cut
where n is the concentration of active enzymes, m is viscosity,and ccut is a microscopic cut-oﬀ length of the order of a protein
size. Moreover, we haveð1dtCðtÞsðtÞ;
(38)
w¼
0

where C(t) is the correlation function of force dipoles corresponding to the enzymes and s(t) is the orientational correlation
function for them; weq is given by the same equation, but with
C(t) replaced by Ceq(t).
The orientational correlation function has the form
(39)

where trot is the orientational correlation time. As seen fromeqn (37)–(39), the magnitude of diffusion enhancement issensitive to the relationship between the correlation time offorce dipole and the orientational correlation time.
According to the Stokes equation, rotational diﬀusion coeﬃcient
for a spherical particle of radius R is
Drot ¼

kB T
:
8pmR3

(40)

The orientational correlation time is defined as trot = 1/Drot.
Since proteins are not spheres, their orientational correlation
times are shorter than given by the Stokes estimate. Even incrowded solutions, they do not exceed a microsecond.40–42 For
active dimers, the orientational correlation times can be estimatedby using eqn (40) with R = c0.
For the active dimer model, the correlation functions of force
dipoles and, therefore, their correlation times were analyticallydetermined at equilibrium [see eqn (17)] and in the limit A [seeeqn (18)]. Moreover, they were also numerically determined, as
shown in Fig. 6. As follows from these results, the correlation times
for force dipoles are determined by conformational relaxation
times t0 and t1. In the discussion of conformational relaxation
phenomena in proteins in Section 2.4, we have noticed that slow
conformational relaxation processes, involving relative domain

10742 | Soft Matter, 2020, 16, 10734--10749

w  weq = trothDm2iA.

(41)

Hence, the change in the diﬀusion coeﬃcient can be estimated as
DA ¼

Diﬀusion eﬀects of enzymes in water solutionssðtÞ = expðt/trotÞ,motions in real enzymes, would usually lie in the microsecond
to millisecond range. Hence, correlation times for force dipolesof the enzymes would be typically longer than their orientational correlation times.
Below, we assume that the orientational correlation time ismuch shorter than the correlation time for force dipoles. By
using eqn (39) and putting C(t) E C(0) = hDm2i in eqn (38), weapproximately findtrot nDm2 A :
60pm2 ‘cut

(42)

Substituting approximate analytical expressions for hDm2iA in
diﬀerent limiting regimes, obtained in Section 2.2, or using thenumerical data from Section 2.3, we can now estimate and analyze
diﬀusion enhancement.
Particularly, it was found in Section 2.3 that hDm2iA for
active dimers has a Michaelis–Menten dependence on thesubstrate concentration. As follows from the above equation,
the same dependence should hold for DA. As could have been
expected, the highest diffusion enhancement is reached under
the condition of substrate saturation for the enzymes. Under
substrate saturation, the intensity of force dipoles depends, as
demonstrated in Sections 2.2 and 2.3, on the properties of theturnover cycle of an enzyme. The highest intensity is found in
regime B, i.e., when there is no long waiting time for productrelease, and the turnover rates are limited by conformationaltransitions within the catalytic cycle (adenylate kinase is an
example of such an enzyme).
In such asymptotic regime, hDm2iA is given by eqn (35)provided that the condition kc02 c kBT holds. Substituting thisexpression into eqn (42), we obtain
DA
nR0
‘0
¼
DT R0 þ ‘0 D‘

2

‘0
‘

3

DE
kB T

2

:

(43)

Here, R0 is the radius of a tracer particle, c0 is the characteristic
size of an enzyme (i.e., the dimer length), Dc specifies themagnitude of the conformational change, DE is the free energy
supplied with the substrate and dissipated within each cycle,
and c = n1/3 is the mean distance between active enzymes. For
the equilibrium diﬀusion constant DT, the Stokes equation
DT ¼

kB T
;
6pmR0

(44)

has been used. Moreover, the microscopic cut-oﬀ length ischosen as ccut = c0 + R0 and n = 4pz1/5 is a numerical factorof order unity.
Eqn (43) can be employed to estimate the maximum relative
diﬀusion enhancement for passive particles that can be obtained,
under substrate saturation, in water solutions of enzymes. For
numerical order-of-magnitude estimates, we consider exothermic
enzymes with DE = 10 kBT and Dc = 0.1c0. As the enzyme
concentration, we take n = 1 mM. This corresponds to a noncrowded solution where the mean distance between the enzymes

This journal is © The Royal Society of Chemistry 2020

View Article Online

Paper

Soft Matter

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

is about ten times larger than their size (c B 10c0). Moreover, weconsider passive particles with the sizes comparable to that of an
enzyme (R0 B c0). Under these conditions, we have DA B 10DT,
i.e., diﬀusion of tracer particles is ten times faster in the solutionof catalytically active enzymes.
Dependence of diﬀusion enhancement on the orientationalcorrelation time is additionally discussed in Appendix D.
3.2 Diﬀusion eﬀects of active protein inclusions inbiomembranes
It is known that, on the length scales shorter than the Saﬀman–
Delbrück length of about a micrometer, lipid bilayers behave as
2D fluids.8 Similar to enzymes in water solutions, active protein
inclusions (such as ion pumps or transporters) can cyclicallychange their shapes inside a lipid bilayer within each ligand
turnover cycle. Hence, they behave as hydrodynamical force
dipoles within a fluid lipid bilayer. Therefore, diffusion enhancement is expected for biomembranes when non-equilibrium conformational activity of proteins takes place.12
A significant diﬀerence to water solution is that, for thebiomembranes as 2D fluids, hydrodynamic diﬀusion enhancementeﬀects are non-local. For such systems, eqn (37) is replaced by12,14
ð ra ra01
DA;aa0 ðRÞ ¼
w

w
(45)
eq dr 4 n2D ðR þ rÞ:
32p2 m2D 2rHere, w is again given by eqn (38) with s(t) being the planar
orientational correlation function for protein inclusions. Moreover,m2D is the 2D viscosity of the lipid bilayer, related as m2D = hm3D to its
3D viscosity m3D (where h is the bilayer thickness); n2D is the 2D
concentration of active inclusions within the membrane.
For numerical estimates, we assume that active proteins
occupy a small circular region (a raft) of radius Rm (shorter thanthe Saﬀman–Delbrück length) within a membrane. Then, diffusion
enhancement for a passive particle of radius R0 located in the centerof the disc is12,14
DA ¼ zmn2Dw  weq ;
2
m2D

(46)

where zm = (1/32p)ln(Rm/ccut), ccut = R0 + c0, and w is given by the
integral in eqn (38) where, however, s(t) is the planar orientationalcorrelation function for proteins inside a membrane.
The viscosity m3D of lipid bilayers is about 103 times higher
than that of water and, therefore, both translational and rotational
diﬀusion is much slower in them. From experiments, it is known
that diﬀusion constants for proteins in lipid bilayers are about
DT = 1010 cm2 s1, i.e., about 103 times smaller than in water
for similar proteins. One can therefore expect that rotational
diﬀusion of proteins in lipid bilayers would be slowed by abouta factor of 103 too, yielding orientational correlation times trotthat might approach a millisecond, still being shorter than theturnover time of an enzyme.
The magnitude of diﬀusion enhancement in eqn (46) can be
determined by modeling protein inclusions as active dimersthat lie flat in the membrane. Then, the same estimate (35)for hDm2iA can be used. Combining all terms, diﬀusion

This journal is © The Royal Society of Chemistry 2020

enhancement in eqn (46) for a passive particle in the centerof a protein raft approximately is
DA ¼ n m trot

‘0
D‘

2

DE
h‘2D m3D

2

;

(47)

where the dimensionless prefactor is nm = z1zm and c2D = n2D1/2
is the mean distance between inclusions in the membrane.
To obtain a characteristic order-of-magnitude estimate, the3D viscosity of the lipid bilayer is chosen as m3D = 1 Pa s and thethickness of the bilayer as h = 1 nm. For protein inclusions, we
assume that DE = 10 kBT and Dc B c0. The orientationalcorrelation time is taken to be trot = 100 ms and the mean
lateral distance between the proteins is c2D = 10 nm. For suchparameter values, the maximal possible diﬀusion enhancement
under substrate saturation conditions is about DA = 109 cm2 s1.
For comparison, Brownian diﬀusion constants for proteins in lipid
bilayers are of the order of 1010 cm2 s1 and diﬀusion constants
for lipids are about 108 cm2 s1.

4 Discussion
Using the results of our study, available experimental andcomputational data on diﬀusion enhancement in solutions ofcatalytically active enzymes can be discussed.
4.1

Experimental data

Diﬀusion enhancement for the enzymes has been reported insolutions of several catalytically active enzymes, at the concentrations varying between 1 nM and 10 nM.36,43–46 With the
exception of aldolase43 (for which, however, the enhancement
could not be independently confirmed47), all these enzymes
were exothermic and had high turnover rates of about 104 s1.
The enhancement was reported not only for the enzymes
themselves, but also for inert molecules (tracers) surrounding
them.44,45 The enzyme concentration dependence of the diffusion
enhancement effects could not however be detected.46
It does not seem plausible that such experimental data canbe understood in the framework of the original theory12 and its
subsequent extensions, including the present work. The factthat a significant diﬀusion enhancement (by tens of percent)
was observed already at low nanomole concentrations can stillbe perhaps explained by assuming that, for some reasons, the force
dipoles of specific enzymes with high catalytic turnover rates were
exceptionally strong. However, the absence of a dependence of theexperimentally observed diﬀusion enhancement on the enzyme
concentration clearly contradicts the theory12 where diﬀusion
enhancement arises as a collective hydrodynamic eﬀect. Eﬀectively,
diﬀusion enhancement was observed in the experiments36,43–46
already for single molecules of enzymes.
When our study was completed, an interesting publication48
has appeared where diﬀusion enhancement was demonstrated
by a diﬀerent method for several other reactions. Since catalyst’s
diﬀusion was not aﬀected by its concentration, this was again asingle-particle property not covered by the theory.12

Soft Matter, 2020, 16, 10734--10749 | 10743

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Soft Matter
Experiments on optical tracking of particles in animal cells49
and in bacteria or yeast50 have been furthermore performed. They
have shown that, when metabolism was suppressed (by depletion ofATP), diﬀusion dropped to undetectable levels49,51 or it was muchslowed down and replaced by subdiﬀusion characteristic for a
colloidal glass.50 Strong reduction of diﬀusion under metabolism
suppression was moreover found in various cytoplasm extracts.52
It should be also noted that diﬀusion enhancement has been
experimentally observed within chromatin in a living biologicalcell.53 This was explained by active operation of molecularmachines involved in transcription and translation of DNA.54
The cytoplasm of a living cell represents a crowded solutionof proteins. In bacteria, the volume fraction of proteins incytosol is about 30 percent,55 with the highest concentrationsof the order of 100 mM reached for glycolysis enzymes. Most of
the enzymes in the cell are mechanochemical, i.e., they exhibit
conformational changes in their catalytic cycles. Typical turnovertimes of enzymes in a biological cell are of the order of 10 ms.
According to the previous12 and current estimates, substantial
diﬀusion enhancement due to hydrodynamic collective eﬀects
should thus be expected under metabolism in the cytoplasm.
There are, however, also other mechanisms that can contributeto diﬀusion enhancement in the cells.
The cytoskeleton of animal cells represents an active gel, withnumerous myosin molecular motors operating within it. It is known
that the activity of the motors can lead to development of nonequilibrium fluctuations in the cytoskeleton which induce in turnfluctuations and diﬀusion enhancement in the cytosol.17,51,56 The
skeleton of bacteria and yeast is however passive; moreover, metabolic diﬀusion enhancement in such cells could also be observed
when their skeleton was chemically resolved.50 Therefore, the active
gel mechanism56 cannot account for the eﬀects observed in them.
On the other hand, under high crowding characteristic for
cytoplasm, proteins are frequently colliding and direct interactions
between them often take place.40,42 It is known that, for dense
colloids, glass behavior can be expected, with the transport andrelaxation phenomena strongly slowed down in them.57 Indeed,such behavior could be observed both in the cells50 and in theextracts52 in the absence of metabolism.
It has been recently shown that, when the particles forming
a glass-like colloid cyclically change their shapes, the colloid gets
fluidized and classical transport properties become restored.58,59
Even in the absence of hydrodynamic interactions, conformational
activity of proteins, at the rates of energy supply of about 10 kBT pera protein molecule per a cycle, can lead to diﬀusion enhancement
by one order of magnitude.58 This provides an additional, nonhydrodynamic, mechanism that can contribute to the experimentally observed diﬀusion enhancement in living biological cells.
In summary, the analysis of the available experimental data
reveals that the predicted diﬀusion enhancement12 for passive
particles caused by collective catalytic activity of enzymes could
not have been reliably confirmed.
4.2

Paper
simulations, the solvent was explicitly included and the multiparticle collision dynamics (MPCD) approximation60 was employed,thus allowing to fully account for hydrodynamic eﬀects.
To facilitate the comparison, we first give a summary of theessential parameter values in the study,20 using the current
notations employed by us. The natural lengths of the dimer in
two ligand states were c0 and c1 = c0/2, and the spring constants
were k0 and k1 = 2k0. The dimensionless spring constant kc02/(kBT), characterizing stiﬀness of the dimer, was varied between
144 and 1440. The energy DE = (1/2)(k0 + k1)(c0  c1)2, supplied
to a dimer and dissipated by it as heat within a single cycle, waschanging therefore between 121.5 kBT and 1215 kBT. The simulations were performed under substrate saturation conditions. Product formation and release were possible within a window of halfwidth r = 0.025c0 near x = c1. The rate v1 of this transition could bevaried in the simulations by a factor of 5.
The Langevin eqn (5) with viscous friction and thermal noise
was not used. Instead, collisions between the two beads of the
dimer and the solvent particles were explicitly taken into
account in the framework of MPCD. For a single passive dimer,
the equilibrium correlation function of force dipoles Ceq(t) was
computed yielding the correlation time for fluctuations of its
force dipole; this function could be well fitted to the theoreticaldependence in eqn (17). Note that, when k0c02/(kBT) c 1, therelaxation time t0 = (gk0)1 of the dimer should be close to thiscorrelation time. Moreover, we have t1 = (gk1)1 = t0/2. Usingsuch estimates, it can be shown that w1t1 varied between 0.001
and 0.1 in the simulations.20 Because substrate saturation was
assumed, conditions w0t0 c 1 and w1t1 { 1 corresponding to
the limit D in Section 2.2 were therefore approximately satisfied.
For single active dimers, correlation functions C(t) of force
dipoles were determined.20 They showed damped oscillations
and were similar to the correlation function for v0t0 = 3 inFig. 6. The correlation times varied, but remained of the sameorder of magnitude as the correlation time of the passive dimer.
The force-dipole intensity hDm2i of active dimers was by about
an order of magnitude larger than hDm2ieq for the passive ones.
Depending on the parameters, it scaled as ka0 with the exponent
a in the range between 1.2 and 1.6, comparable with theexponent of 1.5 in eqn (29).
Orientational correlation functions s(t) were furthermore
computed for single dimers.20 Remarkably, it was found that theorientational correlation time trot was sensitive to the conformational activity of the dimer, getting shorter by about an orderof magnitude when such activity was switched on. Nonetheless,in all simulations trot was larger than the force dipole
correlation time.
Multiparticle 3D computer simulations of colloids formedby active dimers were further performed.20 In the simulations,the truncated potential"#
2r0 48
2r0 24 1(48)
þ ;uðrÞ ¼ 4e
4rr

Computational data

Large-scale computer simulations for colloids of active dimers
have been performed by Dennison, Kapral and Stark.20 In these10744 | Soft Matter, 2020, 16, 10734--10749

for r o 21/24(2r0) and zero otherwise, with e = 2.5 kBT and r0 =1.075c0, was used to describe steep repulsive interactions

This journal is © The Royal Society of Chemistry 2020

View Article Online

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Paper
between the beads belonging to diﬀerent dimers. The interaction radius r0 was chosen as defining the radius of a bead.
Since distances c0 and c1 = 0.5c0 in the open and closed
dimer conformations were both smaller than 2r0 = 2.15c0, large
overlaps between the beads in a dimer were present in thesimulations. However, this did not aﬀect the internal dimer
dynamics because there were no repulsive interactions between
the beads in the same dimer. Additionally, the simulated
system included one passive tracer particle of radius 0.5c0.
The volume fraction f occupied by dimers was determinedby taking into account the overlaps, but assuming that all
dimers were in the equilibrium open state with the length ofc0. Because, under substrate saturation conditions, they were
however mainly found in the closed state with an even stronger
overlap, such definition overestimated the actual volume fraction
by a factor of up to two.
Due to the crowding eﬀects, diﬀusion of a passive particle inthe system of inactive dimers decreased with the volume
fraction of them. The diﬀusion reduction at the highest taken
volume fraction f = 0.266 was less than ten percent, indicatingthat this colloidal system was still far from the glass transition
threshold.57
When the dimers were active, diﬀusion of tracers wasincreasing instead with the dimer volume fraction f. For themost stiﬀ active dimers with k0c02/(kBT) = 1440 and the kinetic
regime with w1t1 about 0.1, relative diﬀusion enhancement of
DA/DT = 0.3 could be observed20 at the dimer volume fraction of
f = 0.266. For the least stiﬀ dimers with k0c02/(kBT) = 144,
diﬀusion enhancement by 5 percent was seen at f = 0.133.
Thus, collective hydrodynamic eﬀects of active enzymes on
diﬀusion of passive particles could be computationally confirmed. To speed up the calculations, model enzymes in the
study20 were chosen however to be unusually rapid (with theturnover times shorter than the rotational diﬀusion time) andunusually exothermic (with the heat release of hundreds of kBT
per a turnover cycle). It would be therefore important to undertake such simulations also for the parameters closer approaching
those of the real enzymes.

Soft Matter
Using the derived statistical properties of force dipoles inthe dimer model, more accurate estimates for diﬀusion enhancement for surrounding passive particles in solutions of active
enzymes were obtained.
Based on these results, currently available experimental andcomputational data has been examined. We have concluded
that, while the collective hydrodynamic eﬀects of diﬀusion
enhancement have been principally confirmed in the computational study,20 further work is needed to bring simulations
closer to the parameter region corresponding to real enzymes.
On the experimental side, we have concluded that the dataon diﬀusion enhancement in weak nanomole solutions of several
fast exothermic enzymes cannot be explained in the framework of
the theory12 and alternative explanations for them should besought. In experimental studies of diﬀusion phenomena in living
cells and in cellular extracts, additional work is needed to distinguish possible hydrodynamic contributions from the eﬀects ofdirect collisions between active proteins and the resulting kinetic
crowding eﬀects. Large-scale numerical simulations of crowded
active colloids including hydrodynamic interactions between theparticles are to be performed. It should be also pointed out that,
although the eﬀects of diﬀusion enhancement are also predicted
for biomembranes crowded with active protein inclusions, experiments and numerical multiparticle simulations of such phenomena
are still missing today; it would be interesting to carry them out.
It was not the aim of the present work to provide a review of
all proposed mechanisms for diﬀusion enhancement eﬀects.
Especially, we have not considered possible origins of diﬀusion
enhancement for single catalytically active enzymes, eventhough this question attracts much attention in view of therecent research.45,46,48 Our focus was on diﬀusion enhancement for passive particles caused by hydrodynamic collective
stirring of the solution by a population of active particles thatcyclically change their shapes, but do not propel themselves.
In the future, the active dimer model can be used to develop
stochastic thermodynamics of mechanochemical enzymes. It
would be important to investigate in detail hydrodynamic
eﬀects, accompanying functional conformational transitions,in all-atom or coarse-grained molecular dynamics simulations
for specific enzymes and protein machines.

5 Conclusions and outlook
To our knowledge, the present work is the first study wherehydrodynamic force dipoles of mechanochemical enzymes havebeen systematically analyzed. Although the analysis has been
performed for an idealized model, order-of-magnitude estimates
for the intensity of such dipoles for characteristic enzymes, suchas adenylate kinase, and for protein machines, such as myosin,
have been obtained.
We have also examined for what kinds of enzymes stronghydrodynamic eﬀects may be expected. Our analysis reveals
that, in the framework of the investigated model, these shouldbe very rapid, highly exothermic and stiﬀ enzymes, where theenergy is dissipated in mechanical motions of a small amplitude. It is interesting to note that these general conditions areindeed satisfied, for example, for catalase or urease.

This journal is © The Royal Society of Chemistry 2020

Conflicts of interest
There are no conflicts to declare.

A Transition probabilities
When transitions between the states s = 0 and s = 1 are rare, thesolution of the master eqn (7) and (8) can be approximatelysought in the form
psðx,tÞ = psðtÞpðx,t|sÞ,(49)

where ps(t) is the probability to find the dimer in the ligand
state s and p(x,t|s) is the probability distribution for distance x
provided that the dimer is (permanently) in the state s.

Soft Matter, 2020, 16, 10734--10749 | 10745

View Article Online

Soft Matter

Paper

Substituting these expressions into eqn (7) and (8) andintegrating over the variable x, one finds that the probabilities
ps obey classical master equations for a two-level system,
dp0¼ w1 p1  w0 p0 ;
dt

(50)

by using eqn (7) and (8) and integrating by parts, we
find
ð1
d x2
dx½k0 xð‘0  xÞp0 ðxÞ þ k1 xð‘1  xÞp1 ðxÞ
¼ 2g
dt
1
ð1
(60)
þ 2gkB T
dx½p0 ðxÞ þ p1 ðxÞ
1

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

anddp1¼ w0 p0  w1 p1 :dt

¼ 2ghmi þ 2gkB T ¼ 0:
(51)Here w0 and w1 are eﬀective rates of transitions between thestates given by
ð1
w0 ¼
dxu0 ðxÞpðxjs ¼ 0Þ;(52)Thus, we straightforwardly obtain that, for an active dimer inany statistically steady state, hmi = kBT.
Note that here and also in the equations below, the integration
limits over x are taken as +N and N. The actual limits areautomatically selected by probability distributions p0(x) and p1(x).

1

C Force-dipole correlation function

and
w1 ¼

ð1

Introducing
dxu1 ðxÞpðxjs ¼ 1Þ:(53)

p0 ðx; tÞ1!

pðx; tÞ ¼

The involved probability distributions in the statisticallystationary case aresﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ


k0
k0
ðx  ‘0 Þ2 ;exp pðxjs ¼ 0Þ ¼
2pkB T
2kB T

;

(61)

p1 ðx; tÞ
we can write the system of two master eqn (7) and (8)concisely as
dp^¼ Lp;
dt

(54)

(62)

where

and
sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ


k1
k1
exp 
ðx  ‘1 Þ2 :pðxjs ¼ 1Þ ¼
2pkB T
2kB T

0
^¼@
L

(55)

L^00

L^01

L^10

L^11

1
A;

(63)

and
If the transition windows are narrow, approximations in
eqn (9) can furthermore be used, so that we obtainw0 = n0 pðx = c0|s = 0Þ,

w1 = n1 pðx = c1|s = 1Þ.

(56)

2

@
@
L^00 ¼ gk0 ð‘0  xÞ  gkB T 2 þ u0 ðxÞ;
@x
@x
and
2

Thus, using the above expressions for distance distributions, we finally getsﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
k0;
w0 ¼ 2rv0
2pkB T(64)@@
L^11 ¼ gk1 ð‘1  xÞ  gkB T 2 þ u1 ðxÞ;
@x
@x

(65)

and
(57)

L̂01 = u1ðxÞ,L̂10 = u0ðxÞ.

(66)

The general solution of eqn (62) isandsﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
k1:
w1 ¼ 2rv1
2pkB Tps ðx; tÞ ¼
(58)

w1
;
w0 þ w1

p1 ¼

w0:
w0 þ w1

L̂qðnÞ = lnqðnÞ,(59)

B Average force dipole
Let us consider the second statistical moment hx2i. In a
steady state, its time derivative is zero. On the other hand,

10746 | Soft Matter, 2020, 16, 10734--10749

ln t
An qðnÞ
þ c:c:s ðxÞe(67)n¼0

where ln and q(n) are eigenvalues and eigenvectors of the linear
operator L̂,In the steady state, the probabilities are
p0 ¼

1
X

(68)

and decomposition coeﬃcients An are determined by initialconditions.
Because the master equation must have a stable stationary
solution, the operator L̂ should always possess a zero eigenvalue
l0 = 0 and, furthermore, condition Reln 4 0 should hold forall other eigenvalues n.61 Generally, the eigenvectors can be

This journal is © The Royal Society of Chemistry 2020

View Article Online

Paper

Soft Matter

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

ordered according to the increase of Reln (and therefore we canenumerate the eigenvalues in such a way that 0 o Re l1 r
Re l2 r Re l3 r . . .). The stationary probability distribution p(x)%coincides with the eigenvector q(0)(x).
The conditional probability G(x,s,t|x0,s0) gives the probability
to find the dimer in various states (x,s) at time t provided that it
was in the state (x0,s0) at time t = 0. It represents a special
solution of the master eqn (62) given by
Gðx; s; tjx0 ; s0 Þ ¼

1
X

ln t
an ðx0 ; s0 ÞqðnÞ
þ c:c:s ðxÞewhere an(x0,s0) are the coeﬃcients of decomposition of thisinitial condition over eigenvectors q(n).
The force dipole m depends on the distance x between thedomains and on the dimer state s, i.e., m(t) = m(x(t),s(t)).
Therefore, in the statistically stationary state we have
X ð1

ð1
dx0

s;s0 ¼0;1 1dxmðx0 ; s0 Þmðx; sÞ

1

(70)

 ps0 ðx0 ÞGðx; s; tjx0 ; s0 Þ:

1X

Bn eln jtj þ c:c:

(71)

n¼1

where the complex coeﬃcients Bn are
Bn ¼

X ð1

ð1
dx0

s;s0 ¼ 0;1 1dxmðx0 ; s0 Þmðx; sÞ

1

(72)

 ps0 ðx0 Þan ðx0 ; s0 ÞqðnÞs ðxÞ:
If we retain in this decomposition only the first, most slowly
decaying term, this yieldsCðtÞ  B1 el1 jtj þ c:c: ¼Cð0Þ GjtjecosðOjtj  aÞ:cos a

ð1=trot þ GÞ2 þO2

Dm2 :

(75)

This yields a non-monotonous dependence of w on the orientationalcorrelation time. If the phase shift a is small and can be neglected(cf. Fig. 7), the maximum value wmax is reached at trot = (O  G)1
and we havewmax G2 þ O2
;
¼
w1
2GO

(76)

where
G
Dm2 ;
G2 þ O2(77)

is the limit of w when trot c G1 and trot c O1.
These results allow us to discuss how the diﬀusion enhancement
would depend on the orientational correlational time trot, notassuming that it is much shorter than the correlation time for force
dipoles. If the approximation in eqn (30) holds, diﬀusion enhancement is determined by eqn (37) where w is given by eqn (75). The
diﬀusion enhancement depends non-monotonously on the orientational correlation time. It increases linearly with trot at short times,
then reaches a maximum at trot = (O  G)1 and finally saturates atlarge orientational correlation times. For example, if we take thevalues G E 1/(2t0) and O E p/(3t0) corresponding to substratesaturation in Fig. 7, the maximum diﬀusion enhancement
would be reached at trot = 1.8t0 and, at the maximum, it will
be larger by about 30 percent than in the limit trot c t0.

Acknowledgements(74)

where G = Re l1, O = Im l1, and B1 = C(0)eia/cos a.
Our numerical simulations, described in Section 2.3, haveshown that, in the regimes approaching a deterministic oscillatory
dimer, the correlation functions of force dipoles could be wellfitted to the above dependence. This suggests that contributions
from the higher, more rapidly decaying relaxation modes n 4 1
have been indeed relatively small. As generally known,62 noisy
oscillators possess a slowly relaxing mode that corresponds to
diﬀusion of the oscillation phase. It can be expected that, underchosen conditions, such a mode has been dominating the correlation functions for oscillatory dimers.

This journal is © The Royal Society of Chemistry 2020

1=trot þ G þ O tan a(73)

Therefore, the normalized correlation function isCðtÞ
1 Gjtj¼ecosðOjtj  aÞ;
Cð0Þ cos a

w¼

w1 ¼

By using eqn (69) and (70), we find that, in the statisticallystationary state, the correlation function of force dipoles isCðtÞ ¼ hmðtÞmð0Þi  m2 ¼

Suppose that the force-dipole correlation function C(t) and theorientational correlation function s(t) are given by eqn (74) and(39). By taking the integral in eqn (38), we find

(69)

n¼0

hmðtÞmð0Þi ¼

D Dependence on orientationalcorrelation time

Stimulating discussions with K. K. Dey, R. Kapral, H. Kitahata
and Y. Koyano are gratefully acknowledged. Y. H. acknowledges
support by a Grant-in-Aid for JSPS Fellows (Grant No. 19J20271)
from the Japan Society for the Promotion of Science (JSPS).
Y. H. also thanks for the hospitality of the Fritz Haber Instituteof the Max Planck Society, where part of this research was
conducted in the support of the scholarship from TMU. S. K.
acknowledges support by a Grant-in-Aid for Scientific Research(C) (Grant No. 18K03567) from the JSPS, and support by
a Grant-in-Aid for Scientific Research on Innovative Areas‘‘Information Physics of Living Matters’’ (Grant No. 20H05538)
from the Ministry of Education, Culture, Sports, Science andTechnology of Japan. S. K. and A. S. M. acknowledge the support
by Grant-in-Aid for Scientific Research (C) (Grant No. 19K03765)
from the JSPS.

Soft Matter, 2020, 16, 10734--10749 | 10747

View Article Online

Soft Matter

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

References
1 A. Cressman, Y. Togashi, A. S. Mikhailov and R. Kapral, Phys.
Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 050901.
2 C. Echeverria, Y. Togashi, A. S. Mikhailov and R. Kapral,
Phys. Chem. Chem. Phys., 2011, 13, 10527.
3 R. Golestanian and A. Ajdari, Phys. Rev. Lett., 2008, 100, 038101.
4 R. Golestanian and A. Ajdari, J. Phys.: Condens. Matter, 2009,21, 204104.
5 M. Iima and A. S. Mikhailov, EPL, 2009, 85, 44001.
6 T. Sakaue, R. Kapral and A. S. Mikhailov, Eur. Phys. J. B,
2010, 75, 381.
7 X. Bai and P. Wolynes, J. Chem. Phys., 2015, 143, 165101.
8 H. Diamant, J. Phys. Soc. Jpn., 2009, 78, 041002.
9 M.-J. Huang, A. S. Mikhailov, R. Kapral and H.-Y. Chen,
J. Chem. Phys., 2012, 137, 055101.
10 M.-J. Huang, A. S. Mikhailov, R. Kapral and H.-Y. Chen,
J. Chem. Phys., 2013, 138, 195101.
11 M.-J. Huang, H.-Y. Chen and A. S. Mikhailov, Eur. Phys. J. E:
Soft Matter Biol. Phys., 2012, 35, 119.
12 A. S. Mikhailov and R. Kapral, Proc. Natl. Acad. Sci. U. S. A.,
2015, 112, E3639.
13 R. Kapral and A. S. Mikhailov, Phys. D, 2016, 318–319, 104.
14 A. S. Mikhailov, Y. Koyano and H. Kitahata, J. Phys. Soc. Jpn.,
2017, 86, 101013.
15 Y. Koyano, H. Kitahata and A. S. Mikhailov, Phys. Rev. E,
2016, 94, 022416.
16 Y. Hosaka, K. Yasuda, R. Okamoto and S. Komura, Phys. Rev.
E, 2017, 95, 052407.
17 K. Yasuda, R. Okamoto, S. Komura and A. S. Mikhailov, EPL,
2017, 117, 38001.
18 K. Yasuda, R. Okamoto and S. Komura, Phys. Rev. E, 2017,95, 032417.
19 Y. Hosaka, S. Komura and D. Andelman, Phys. Rev. E, 2020,
101, 012610.
20 M. Dennison, R. Kapral and H. Stark, Soft Matter, 2017,13, 3741.
21 F. Kogler, Interactions of artificial molecular machines,
Diploma thesis, Technical Univ. of Berlin, Berlin, 2012.
22 P. Illien, T. Adeleke-Larodo and R. Golestanian, EPL, 2017,
119, 40002.
23 H. Flechsig and A. S. Mikhailov, J. R. Soc., Interface, 2019,16, 20190244.
24 J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held,J. D. Chodera, C. Schütte and F. Noé, J. Chem. Phys., 2010,
134, 174105.
25 D. E. Shaw, P. Maragakis, K. Lindorﬀ-Larsen, S. Piana,
R. O. Dror, M. P. Eastwood, J. A. Bank, J. M. Jumper,
J. K. Salmon, Y. Shan and W. Wriggers, Science, 2010,
330, 341.
26 H. P. Lu, L. Xun and X. S. Xie, Science, 1998, 282, 1887.
27 H.-P. Lerch, R. Rigler and A. S. Mikhailov, Proc. Natl. Acad.
Sci. U. S. A., 2005, 102, 10807.
28 M. Gur, E. Zomot and I. Bahar, J. Chem. Phys., 2013, 139, 121912.
29 K. A. Henzler-Wildman, V. Thai, M. Lei, M. Ott, M. WolfWatz, T. Fenn, E. Pozharski, M. A. Wilson, G. A. Petsko,

10748 | Soft Matter, 2020, 16, 10734--10749

Paper
M. Karplus, C. G. Hübner and D. Kern, Nature, 2007,
450, 838.
30 Y. Togashi and A. S. Mikhailov, Proc. Natl. Acad. Sci. U. S. A.,
2007, 104, 8697.
31 H. Flechsig and A. S. Mikhailov, Proc. Natl. Acad. Sci. U. S. A.,
2013, 107, 20875.
32 M. Hinczewski, R. Tehver and D. Thirumalai, Proc. Natl.
Acad. Sci. U. S. A., 2013, 110, E4059.
33 K. Shiroguchi, H. F. Chin, D. E. Hannemann, E. Muneyuki,
E. M. D. L. Cruz and K. Kinosita Jr., Proc. Natl. Acad. Sci. U. S. A.,
2011, 9, e1001031.
34 N. Kodera, D. Yamamoto, R. Ishikawa and T. Ando, Nature,
2010, 468, 72.
35 J. M. Berg, L. Stryer, J. Tymoczko and G. Gatto, Biochemistry,
WH Freeman, 9th edn, 2019.
36 C. Riedel, R. Gabizon, C. A. M. Wilson, K. Hamadani,
K. Tsekouras, S. Marqusee, S. Pressé and C. Bustamante,
Nature, 2015, 517, 227.
37 R. Inoue, R. Biehl, T. Rosenkranz, J. Fitter, M. Monkenbusch,
A. Radulescu, B. Farago and D. Richter, Biophys. J., 2010,
99, 2309.
38 J. Schonfield, P. Inder and R. Kapral, J. Chem. Phys., 2012,
136, 205101.
39 P. Gouet, H.-M. Jouve, P. A. Williams, I. Andersson,
P. Andreoletti, L. Nussaume and J. Hajdu, Nat. Struct. Biol.,
1996, 3, 951.
40 S. von Bülow, M. Siggel, M. Linke and G. Hummer, Proc.
Natl. Acad. Sci. U. S. A., 2019, 116, 9843.
41 Z. Bashardanesh, J. Elf, H. Zhang and D. van der Spoel, ACS
Omega, 2019, 4, 20654.
42 G. Nawrocki, A. Karaboga, Y. Sugita and M. Feig, Phys.
Chem. Chem. Phys., 2019, 21, 876.
43 P. Illien, X. Zhao, K. Dey, P. J. Butler, A. Sen and
R. Golestanian, Nano Lett., 2017, 17, 4415.
44 K. K. Dey, Angew. Chem., Int. Ed., 2019, 58, 2208.
45 A.-Y. Jee, S. Dutta, Y.-K. Cho, T. Tlusty and S. Granick, Proc.
Natl. Acad. Sci. U. S. A., 2018, 115, 14.
46 M. Xu, J. L. Ross, L. Valdez and A. Sen, Phys. Rev. Lett., 2019,
123, 128101.
47 Y. Zhang, M. J. Armstrong, N. M. B. Kazeruni and H. Hess,
Nano Lett., 2018, 18, 8025.
48 H. Wang, M. Park, R. Dong, J. Kim, Y.-K. Cho, T. Tlusty and
S. Granick, Science, 2020, 369, 537.
49 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, 2014, 158, 822.
50 B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O’Hem,
E. R. Dufresne and C. Jacobs-Wagner, Cell, 2014, 156, 183.
51 E. Fodor, M. Guo, N. S. Gov, P. Visco, D. A. Weitz and F. van
Wijland, EPL, 2015, 110, 48005.
52 K. Nishizawa, K. Fujiwara, N. Ikenaga, N. Nakajo and
D. Mizuno, Sci. Rep., 2017, 7, 15143.
53 S. C. Weber, A. J. Spakowitz and J. A. Theriot, Proc. Natl.
Acad. Sci. U. S. A., 2012, 109, 7338.
54 R. Bruinsma, A. Y. Grosberg, Y. Rabin and A. Zidovska,
Biophys. J., 2014, 106, 1871.

This journal is © The Royal Society of Chemistry 2020

View Article Online

Paper

59 N. Oyama, T. Kawasaki, H. Mizuno and A. Ikeda, Phys. Rev.
Res., 2019, 1, 032038.
60 R. Kapral, Adv. Chem. Phys., 2008, 140, 89.
61 H. Risken, The Fokker-Planck Equation: Methods of Solution
and Applications, Springer, Berlin, 1989.
62 Y. Kuramoto, Chemical Oscillations, Waves and Turbulence,
Springer, Berlin, 1984.

Open Access Article. Published on 27 October 2020. Downloaded on 12/16/2020 1:31:29 PM.
This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

55 A. Vendeville, D. Lariviere and E. Fourmentin, FEMS Microbiol. Rev., 2010, 35, 395.
56 F. C. MacKintosh and A. J. Levine, Phys. Rev. Lett., 2008,
100, 018104.
57 G. L. Hunter and E. R. Weeks, Rep. Prog. Phys., 2012, 75, 066501.
58 Y. Koyano, H. Kitahata and A. S. Mikhailov, EPL, 2019,
128, 40003.

Soft Matter

This journal is © The Royal Society of Chemistry 2020

Soft Matter, 2020, 16, 10734--10749 | 10749