PHYSICAL REVIEW E 103, 062113 (2021)

Nonreciprocality of a micromachine driven by a catalytic chemical reaction
Kento Yasuda* and Shigeyuki Komura†

Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 28 November 2020; revised 4 March 2021; accepted 20 April 2021; published 9 June 2021)
We propose a model that describes cyclic state transitions of a micromachine driven by a catalytic chemical
reaction. We consider a mechanochemical coupling of variables representing the degree of a chemical reaction
and the internal state of a micromachine. The total free energy consists of a tilted periodic potential anda mechanochemical coupling energy. We assume that the reaction variable obeys a deterministic stepwise
dynamics characterized by two typical timescales, i.e., the mean first passage time and the mean first transition
path time. To estimate the functionality of a micromachine, we focus on the quantity called “nonreciprocality”
and further discuss its dependence on the properties of catalytic reaction. For example, we show that thenonreciprocality is proportional to the square of the mean first transition path time. The explicit calculationof the two timescales within the decoupling approximation model reveals that the nonreciprocality is inversely
proportional to the square of the energy barrier of catalytic reaction.
DOI: 10.1103/PhysRevE.103.062113

I. INTRODUCTION

In recent years, physics of micromachines such as bacteria,
motor proteins, and artificial molecular machines has been
intensively studied [1,2]. Generally, a micromachine can be
defined as a small object that extracts energy from chemical substances in the system and further exhibits mechanicalfunctions. The interplay between the structural dynamics ofsuch a small object and the associated chemical reaction is
crucial for the operation of a micromachine [3,4]. Owing to
the developments in nonequilibrium statistical mechanics andexperimental techniques, various studies have been conducted
to reveal the energetics of a single micromachine. For example, energy efficiencies of F1 -ATPase and kinesin motors
have been experimentally measured by using the Harada-Sasa
relation [5–7].
Furthermore, attention has been paid to the dynamics of
micromachines. For instance, several works reported that
diffusion coefficients of metabolic enzymes increase due
to enzymatic reactions [8–10]. Although various possible
scenarios have been proposed such as self-thermophoresis,
stochastic swimming, or collective heating, the main physicalmechanism for the enhanced diffusion is not yet specified [9].
Moreover, the experiment by Jee et al. showed that metabolic
enzymes can move in a directional manner in the presence of
catalytic reactions [11]. Although such a swimming behavior
can be explained by a specific theoretical model [12], a more
fundamental understanding concerning the interplay between
the dynamics and function of a micromachine is necessary.
Biological functions of a micromachine is intimately related to the transitions between different internal states. As
depicted in Fig. 1(a), one can use time-dependent state*
Present address: Research Institute for Mathematical Sciences,
Kyoto University, Kyoto 606-8502, Japan.
†
komura@tmu.ac.jp

2470-0045/2021/103(6)/062113(11)variables si (t ) (i = 1, 2, 3, . . . ) to characterize such as conformational structure or adhesion state when a micromachineis interacting with a substrate. The state variables si (t ) change
dynamically when a micromachine catalyzes a chemical reaction of substrate molecules. As long as it acts as a catalyst,
however, the internal state should return to the initial state
after one cycle of reaction. Hence, si (t ) should change periodically in time as the chemical reaction proceeds repeatedly.
In overdamped systems, such a cyclic change of internal
state is related to biological functions of a micromachine. For
example, microswimmers in a viscous fluid have been investigated by using specific models such as connected spheres[13,14] or a spherical squirmer [15]. It was shown that the
average swimming velocity is proportional to the closed loop
area in the corresponding deformation space. In order toextend this concept and to generally characterize the functionality of a generic micromachine, we focus on the following
quantity τc
Ri j =
dt ṡi s j ,(1)0

where ṡi = dsi /dt and τc is the period of the cycle. In thispaper, we shall call Ri j the “nonreciprocality” representing thearea enclosed by a trajectory in the state space, as shown in
Fig. 1(b).
For a three-sphere microswimmer [14], s1 and s2 correspond to the lengths of the two arms, and its averageswimming velocity is directly proportional to the nonreciprocality, i.e., V ∼ R12 . This relation was also discussed in
the experimental realization of a three-sphere microswimmer[16]. Such a relation holds not only for deterministic microswimmers, but also for stochastic microswimmers [17,18].
The nonreciprocality Ri j is also relevant to crawling motions of a cell on a substrate [19,20]. Within a gauge theory,
the average velocity of a deformable body is given by V =i, j Wi j Ri j , where Wi j is a weighting tensor that connects
the velocity and the nonreciprocality [21]. It should be noted

062113-1©2021 American Physical Society

KENTO YASUDA AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 103, 062113 (2021)

time τt . Solving the equations for the state variables si , weanalytically obtain the nonreciprocality R12 as a function of
the above two timescales. Furthermore, we obtain analyticalexpressions of τp and τt within the decoupling approximation and relate them to the properties of catalytic reaction.
Combining these results, we show that the nonreciprocality
is inversely proportional to the square of the energy barrier ofcatalytic reaction.
In the next section, we explain the model of a micromachine driven by a catalytic chemical reaction. In Sec. III,
we argue the dynamics of the state variables si . In Sec. IV,
we calculate the nonreciprocality R12 analytically. In Sec. V,
we obtain the mean first passage time τp and the mean firsttransition path time τt . Finally, a summary of our work andsome further comments are given in Sec. VI.
II. MODEL
A. Catalytic chemical reaction

Consider a system which contains one enzyme molecule(E) that acts as a micromachine, nS substrate molecules (S),and nP product molecules (P). The enzyme molecule plays therole of a catalyst and the corresponding chemical reaction is
written as [22]S + E  ES → P + E

FIG. 1. (a) Schematic picture of a micromachine characterized
by the conformational state variables s1 , s2 , and s3 . Moreover, theadhesion between the domains and the substrate is described by thevariables s4 and s5 . (b) The state transition of a micromachine isrepresented by a trajectory in the state space spanned by the variables
si and s j . The nonreciprocality Ri j [see Eq. (1)] represents the area
enclosed by the trajectory.

that the nonreciprocality Ri j is a universal quantity that does
not depend on specific self-propulsion models. Moreover, Ri jquantifies how much a micromachine breaks the detailed balance that needs to be satisfied in thermal equilibrium.
Although the functionality of a micromachine can be quantified by Ri j , it is important to clarify how the nonreciprocality
is regulated within a micromachine. Considering a micromachine that utilizes a catalytic chemical reaction, we investigate
the relationship between the properties of the reaction andthe nonreciprocality Ri j . Our main purpose is to propose aminimum model of a micromachine undergoing cyclic statetransitions which are driven by repeated catalytic reactions.
Hence, in addition to the state variables si (t ), we introduce
another variable θ (t ) to describe the degree of a catalytic
chemical reaction. These variables are related to each otherthrough the mechanochemical coupling mechanism.
We employ Onsager’s phenomenological equations for thetime evolutions of si (t ) and θ (t ). In order to solve the equations analytically, we consider the weak coupling limit andassume that the dynamics of θ (t ) is described by a step function characterized by two characteristic time scales, i.e., themean first passage time τp and the mean first transition path(2)

where ES indicates a complex molecule. The reaction rateṅP is often analyzed by the Michaelis-Menten equation [22].
Although the above catalytic chemical reaction is relevant to
the present study, our purpose is to adopt the simplest model
for such a chemical reaction and not to reproduce it.
The extent of a catalytic reaction is commonly describedby the number of product molecules nP . However, since our
purpose is to investigate a single molecular reaction process,we introduce a reaction variable θ (t ) to quantify the extent of
catalytic reaction. Unlike the quantity nP , the reaction variableθ is a continuous number and increases 2π for each reaction.
Under this assumption, θ represents the reaction phase of aperiodic catalytic reaction.
According to the Kramers theory, the free energy Gr describing a chemical reaction is given by a tilted periodic
potential [23]Gr (θ ) = Gp (θ ) − F θ ,(3)where Gp is a periodic potential with a period of 2π , i.e.,
Gp (θ + 2π ) = Gp (θ ), as shown schematically in Fig. 2(a).
This is because θ increases by 2π for one cycle of chemical
reaction and should experience the same potential. We also
require that Gp takes minimum values at θ = 2nπ (n being an
integer) because the chemical states should be stable before
and after the catalytic reaction. As shown in Fig. 2(a), theamplitude of Gp , denoted by A, represents the energy barrier inthe chemical reaction and is regarded as the activation energy.
The explicit form of Gp will be presented later in Eq. (23).
On the other hand, F in Eq. (3) represents the chemical potential difference (such as between ATP and ADP molecules)
that drives catalytic reaction. Physically, it corresponds to anonequilibrium force even though F has the dimension ofenergy. The system is in chemical equilibrium when F = 0,
whereas it is in an out-of-equilibrium situation when F = 0.

062113-2

NONRECIPROCALITY OF A MICROMACHINE DRIVEN …

PHYSICAL REVIEW E 103, 062113 (2021)

that is a function of the reaction variable θ . For a catalytic
reaction, the internal state of a micromachine should return to the initial state after one cycle of reaction and thesame process takes place repeatedly. Hence, we consider thatthe natural state i (θ ) changes also periodically and assume the simplest periodic form i (θ ) = di sin(θ + φi ), where
di is the amplitude and φi is the constant phase difference
relative to the reaction phase θ . Under these assumptions, weconsider the following mechanochemical coupling energy Gc
between θ and si :
 Ki
Gc (θ , {si }) =(4)[si − di sin(θ + φi )]2 .
2iThen the total free energy Gt in our model is simply given by
Gt (θ , {si }) = Gr (θ ) + Gc (θ , {si }).

(5)

C. Dynamic equations

For the time evolutions of θ and si , we employ the Onsager’s phenomenological equations [24],∂Gt+ ξ (t ),∂θ∂Gts˙i = −μi j
+ ξi (t ),∂s jjθ̇ = −M

FIG. 2. (a) The periodic potential Gp as a function of the catalytic
reaction variable θ with a period of 2π . Here we show a linear
function given by the Eq. (23) with an energy barrier A. (b) The
tilted periodic potential Gr as a function of the reaction variable θ.
As shown in Eq. (3), Gr consists of the periodic part Gp in panel(a) and the linear part −F θ , where F is the nonequilibrium force. A
possible trajectory of θ is shown by the blue (gray) arrow. The valueof θ fluctuates around the minimum of the potential and the transition
to the next minimum takes place occasionally.

In this paper, we shall consider the case of F > 0. With theadded nonequilibrium force F , the free energy Gr for catalytic
reaction becomes a tilted periodic potential, as schematicallyshown in Fig. 2(b).
B. Mechanochemical coupling

Next, we introduce the state variables si (t ) (i =1, 2, 3, . . . ) characterizing the conformation of a micromachine. As shown in Fig. 1(a), examples of the state variables
are distances between the domains in a micromachine or distances between the domains and the substrate (if it exists). In
the molecular dynamics simulation of myosin V, for example,
the protein structure is characterized by the relative distances
between the three amino acids [3]. In principle, there are alarge number of degrees of freedom of a micromachine, andhence the number of the state variables si can also be large.
Next, we explain the mechanochemical coupling mechanism in our model. We assume that each state variable si
experiences a harmonic potential, Ki [si − i (θ )]2 /2, where
Ki is the coupling parameter and i (θ ) is the natural state

(6)
(7)

where M and μi j are the Onsager coefficients for θ and si ,respectively. These coefficients represent energy dissipation,and μi j is given, for example, by the inverse of the friction
coefficient of a domain due to the surrounding viscous fluid.
Moreover, ξ and ξi represent thermal fluctuations which satisfy the fluctuation-dissipation theorem:ξ (t ) = 0,


(8)


ξ (t )ξ (t ) = 2MkB T δ(t − t ),(9)ξi (t ) = 0,

(10)





ξi (t )ξ j (t ) = 2μi j kB T δ(t − t ),(11)

where kB is the Boltzmann constant and T is the temperature.
In the above equations, the Onsager coefficients are assumed to be constant and thermal fluctuations are given by
Gaussian white noise. In the presence of a memory effect suchas viscoelasticity, the Onsager coefficients depend on time andthermal fluctuations are given by colored noise in the form ofgeneralized fluctuation-dissipation relations. In the absence ofthermal fluctuations, the reaction variable θ does not evolvein time because of the energy barrier A in the potential Gr[see Fig. 2(b)]. If thermal fluctuations are present, the valueof θ fluctuates around the minimum of the potential andthe transition to the next minimum takes place occasionally[blue (gray) trajectory in Fig. 2(b)]. Hence, thermal fluctuations are necessary to drive time evolutions of θ and si in ourstochastic model.
Although our model is general, we make several simplifications in order to solve the coupled equations analytically.
First, we only take into account two degrees of freedom, i.e.,s1 and s2 . Second, the mobility coefficients μi j is assumed tohave the form μ11 = μ22 = μ and μ12 = μ21 = 0. Third, thecoupling free energy is symmetric between the two degrees of

062113-3

KENTO YASUDA AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 103, 062113 (2021)

by thermal fluctuations, it is difficult to isolate the peaks ofsi because they are almost comparable to the backgroundfluctuations. Moreover, the state transitions become very rare
when A  kB T . Therefore, instead of performing further simulations, we investigate our model analytically by using thedecoupling approximation.
D. Decoupling approximation

The decoupling approximation relies on two assumptions:A, and (ii) assum(i) taking the weak coupling limit, Kd 2
ing a deterministic dynamics for θ (t ). Under the assumption(i), Eq. (12) can be simplified to

s2 /d

s1 /d

θ̇ = −M∂θ Gr (θ ) + ξ ,Dt
FIG. 3. Time evolutions of θ (top), s1 (middle), and s2 (bottom)
obtained by numerically solving Eqs. (12)–(14). We use the dimensionless time by Dt, where D = MkB T has the dimension of inverse
time. The parameters are a = A/(kB T ) = 10, f = F/(kB T ) = 8,
Kd 2 /(kB T ) = 10, φ1 = 0, φ2 = π /2, and μ/(d 2 M ) = 1. In this simulation, the periodic potential is taken to be Gp (θ ) = −A cos θ and isdifferent from the linear function in Eq. (23). The red (gray) verticallines indicate the moments when the catalytic chemical reactions takeplace.

freedom, i.e., K1 = K2 = K and d1 = d2 = d. Then Eqs. (6)and (7) reduce toθ̇ = −M[∂θ Gr (θ ) − Kd cos(θ + φ1 )δ1− Kd cos(θ + φ2 )δ2 ] + ξ ,(12)δ̇1 = −γ δ1 − d cos(θ + φ1 )θ̇ + ξ1 ,(13)

δ̇2 = −γ δ2 − d cos(θ + φ2 )θ̇ + ξ2 ,(14)

where we have introduced δi = si − d sin(θ + φi ) and defined
the relaxation rate γ = μK.
As a demonstration of our model, we have numericallysolved Eqs. (12)–(14) by using the Euler method and plottedin Fig. 3 the time evolutions of θ (top), s1 (middle), and s2(bottom) for certain parameters. We see that the reaction variable θ increases stochastically in a stepwise manner, whereas
the state variables s1 and s2 undergo almost random fluctuations. Corresponding to the stepwise increase of θ , both s1 ands2 tend to show peaks as indicated by the vertical red (gray)lines. Although the simulation result in Fig. 3 demonstrates
that the present minimum micromachine is indeed driven(15)

where θ is now decoupled from si . If we were able to solve
Eq. (15) for θ , we can further solve Eqs. (13) and (14) for δi(si ). However, it is still difficult to solve Eq. (15) analyticallybecause θ is mostly trapped in the local minimum of Gr andthermal fluctuations are necessary to overcome the energy
barrier A.
To tackle this problem, we further employ the assumption(ii) for θ (t ). As we shall explain in the next section, we assume
that θ (t ) is described by a deterministic function characterized by two timescales, i.e., the mean first passage time τpand the mean first transition path time τt [see later Eq. (16)and Fig. 4(b)]. Then, we can first solve Eqs. (13) and (14)analytically and obtain the explicit expression for δi , as weshow in Sec. III. This result will be used to compute thenonreciprocality R12 analytically in Sec. IV. In Sec. V, on
the other hand, the two characteristic timescales τp and τtwill be separately calculated by using Eq. (15) within thedecoupling approximation.
III. DYNAMICS OF STATE VARIABLES

In this section, we discuss the dynamics of the statevariables δ1 (s1 ) and δ2 (s2 ) that obey Eqs. (13) and (14),respectively. To solve these equations, we make an assumption
for the time dependence of θ (t), as we mentioned in theprevious section. With a tilted periodic potential given by
Eq. (3) and shown in Fig. 2(b), the reaction variable θ changes
stochastically and increases in a stepwise manner as we saw inthe numerical simulation [see Fig. 3 (top)] and also schematically depicted in Fig. 4(a). Such a time evolution of θ (t ) can
be characterized by two characteristic timescales. The firstone is the “first passage time” tp , which is the time required
to change from one local minimum to the neighboring lower
local minimum [23]. The second one is the “first transition
path time” tt which is the time needed for the actual transition[25–27]. It should be noticed that both tp and tt are stochastic
quantities.
To discuss the dynamics of s1 and s2 , let us assume that
θ can be represented by a deterministic stepwise function
characterized by the “mean first passage time” τp and the“mean first transition path time” τt , which are the averages
of tp and tt , respectively. As depicted in Fig. 4(b), the assumed
functional form of θ is
2πt/τt for 0  t < τt
θ (t ) =.
(16)
2π
for τt  t < τp062113-4

NONRECIPROCALITY OF A MICROMACHINE DRIVEN …

PHYSICAL REVIEW E 103, 062113 (2021)

for 0  t < τt andeγ τp (eγ τt − 1)d(γ τt /2π )2 + 1 eγ τp − 1γ τt
cos φi + sin φi e−γ t ,×2πδi (t ) = −FIG. 4. (a) Schematic example of a stochastic time evolution of
the reaction variable θ . Such a time evolution is characterized by the
first transition path time tt and the first passage time tp , which are bothstochastic quantities. (b) After averaging over these quantities, we
obtain the average time evolution of θ as assumed in Eq. (16). Hereτt and τp are the mean first transition path time and the mean firstpassage time, respectively. The reaction variable increases linearlyfor 0  t < τt and remains constant for τt  t < τp . The reaction
phase θ (t ) satisfies θ (t + τp ) = θ (t ) + 2π .

Furthermore, we require that θ increases by 2π after one
cycle of catalytic reaction τp , i.e., θ (t + τp ) = θ (t ) + 2π . The
explicit expressions of τp and τt under the decoupling approximation will be given in Sec. V, where we focus on theirdependencies on the energy barrier A and the nonequilibrium
force F .
Substituting Eq. (16) into Eqs. (13) and (14), we solve
them in the absence of thermal noise, i.e., ξ1 = ξ2 = 0 (seeAppendix A for the details). Then the stationary solution forδi (i = 1, 2) can be obtained as

δi (t ) =d(γ τt /2π )2 + 1





2πt
γ τt
2πt
× −
cos
+ φi − sin
+ φi
2π
τt
τt
+

eγ τp − eγ τt γ τtcos φi + sin φi e−γ t ,eγ τp − 1 2π

(17)

(18)

for τt  t < τp .
In Figs. 5(a) and 5(b), we plot the time evolutions of s1 ands2 , respectively, when φ1 = 0 and φ2 = π /2. We set τt /τp =0.5 and change γ τp = 1, 10, and 102 . Notice that γ τp = μKτp
is the dimensionless relaxation rate. When the relaxation isfast, γ τp > 1, the internal state of a micromachine can be
sufficiently relaxed to the initial state within the reaction cycle
τp . When the relaxation is slow, γ τp < 1, on the other hand,the next reaction starts before the internal state is fully relaxed.
The fast relaxation cases, γ τp = 10, 102 adopted in Fig. 5
do not contradict with the weak coupling limit, Kd 2
A,
as we have discussed before. These two conditions can be
simultaneously satisfied when μ  d 2 /(τp A), namely, when
the Onsager coefficient μ is sufficiently large.
When γ τp = 102 (green or light gray), the behaviors for0  t < τt are well described by sinusoidal functions s1 /d =sin(2πt/τt ) and s2 /d = cos(2πt/τt ). In this fast relaxationcase, we see a sufficiently large state change within a micromachine. When γ τp = 1 (black) or 10 (red or gray), si cannot
follow the change in θ and the functionality of a micromachine is diminished. In Fig. 5(c), we plot the trajectories ofs1 and s2 over one cycle of reaction for different values ofγ τp . For each cycle, 20 equal time intervals are marked by thefilled circles. As mentioned before, the enclosed area of eachtrajectory gives the nonreciprocality R12 .
In Figs. 5(d)–5(f), we show the corresponding plots whenτt /τp = 0.1 (smaller τt ). When γ τp = 102 (green or lightgray) and 0  t < τt , both s1 and s2 are well described by thesame sinusoidal functions as in Figs. 5(a) and 5(b), respectively. On the other hand, the black and red (gray) curves forγ τp = 1 and 10, respectively, deviate significantly from thegreen (light gray) curve for γ τp = 102 , and the magnitudes are
significantly suppressed. The reduced magnitudes can also beseen in Fig. 5(f), where the areas enclosed by the black and red(gray) lines are much smaller than that of the green (light gray)line. This means that, for γ τp = 1 and 10, the nonreciprocality
is further decreased as τt /τp is made smaller.
IV. NONRECIPROCALITY OF A MICROMACHINE

We have mentioned in introduction that the nonreciprocality defined in Eq. (1) provides us with a useful quantity
to evaluate the functionality of a micromachine [14,17,19–
21]. Previously, the nonreciprocality was obtained along a
deterministic state change when the period of deformationis constant. However, this is not always possible when fluctuations are present [28]. For a stochastic micromachine, itis necessary either to take a long time limit or to estimate
the statistical average to estimate the nonreciprocality. In thepresent model and analysis, on the other hand, one can calculate the nonreciprocality directly from Eq. (1) because wehave assumed a deterministic dynamics for θ as in Eq. (16).

062113-5

KENTO YASUDA AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 103, 062113 (2021)γτp = 1γτp = 10γτp = 102

s2 /d

s1 /d

0.5

0.0

1.0

1.0

0.5

0.5

s2 /d

1.0

0.0

γτp = 1γτp = 10γτp = 102

-0.5

-0.5

(b)

(a)

-1.0
0.0

0.2

0.4

0.6

0.8

-1.0
0.0

1.0

(c)

0.0

-0.5

-1.0
0.2

0.4

0.6

0.8

-1.0

1.0

-0.5

1.0

1.0

(d)

0.0

0.5

1.0

0.5

1.0

s1 /d

t/τp

t/τp

1.0

(f)

(e)
0.5

0.0

γτp = 1γτp = 10γτp = 102

-0.5

-1.0
0.0

0.2

0.4

0.6

0.8

0.5

s2 /d

s2 /d

s1 /d

0.5

0.0

γτp = 1γτp = 10γτp = 102

-0.5

1.0

-1.0
0.0

0.0

-0.5

-1.0

0.2

0.4

t/τp

0.6

0.8

1.0

-1.0

-0.5

t/τp

0.0

s1 /d

FIG. 5. Time evolutions of (a) s1 and (b) s2 when φ1 = 0 and φ2 = π /2. We set τt /τp = 0.5 (shown by the dashed line) and change γ τp = 1(black), 10 (red or gray), and 102 (green or light gray). (c) The trajectory of the state variables s1 and s2 shown in panels (a) and (b), respectively,
over one cycle. For each cycle, 20 equal time intervals are marked by the filled circles. The black arrow indicates the direction of the statetransition. The enclosed area of the trajectory corresponds to the nonreciprocality R12 . Plots (d)–(f) are similar plots to (a)–(c), respectively,
when τt /τp = 0.1 (shown by the dashed line).

With the use of Eqs. (17) and (18), the nonreciprocality R12
can be analytically obtained in terms of τp and τt as τp
dt ṡ1 s2
R12 =In Fig. 6, we plot the dimensionless nonreciprocality r12 , i.e.,
R12 scaled by the right-hand side of Eq. (20), as a function of
the ratio τt /τp for different values of γ τp . The dotted lines

0=d 2 γ 2 τt2sin(φ2 − φ1 )4π [(γ τt /2π )2 + 1]

2[1 + eγ τp − eγ τt − eγ (τp −τt ) ].
× 1+
γ τt [(γ τt /2π )2 + 1](eγ τp − 1)

1.0
0.8

(19)
0.6

R12 ≈
R12 ≈

d 2 γ 2 τp2 sin(φ2 − φ1 )4π [(γ τp /2π )2 + 1]

(τt /τp ≈ 1),

3d 2 γ 2 τt2sin(φ2 − φ1 ) (τt /τp
4π

1).

(20)
(21)

r12

This is the main result of this paper. Since R12 is proportional
to sin(φ2 − φ1 ), it vanishes when φ1 = φ2 . In other words,
the state variables s1 and s2 should be out of phase (φ1 = φ2 )with respect to each other in order to exhibit a functionality.
This result is in accordance with the scallop theorem for a microswimmer [14,21]. Moreover, the nonreciprocality satisfies
the symmetry property such that R12 = −R21 .
From Eq. (19), the asymptotic expressions of R12 can beobtained as

0.4γτp = 1γτp = 10γτp = 102

0.2
0.0
0.0

0.2

0.4

0.6

0.8

1.0

τt /τpFIG. 6. The dimensionless nonreciprocality r12 (defined in thetext) as a function of τt /τp for γ τp = 1 (black), 10 (red or gray),and 102 (green or light gray). The dotted lines are the asymptotic
1 given by Eq. (21).
expression for τt /τp062113-6

NONRECIPROCALITY OF A MICROMACHINE DRIVEN …

PHYSICAL REVIEW E 103, 062113 (2021)

represent the asymptotic expression in Eq. (21). From thisplot, one can confirm the scaling behavior R12 ∼ (γ τt )2 whenτt /τp1. When τt /τp ≈ 1, on the other hand, r12 approaches
unity as we see in Eq. (20).
V. TWO CHARACTERISTIC TIMESCALESAs discussed in Sec. III, the dynamics of a catalytic chemical reaction is generally characterized by the mean firstpassage time τp and the mean first transition path time τt . According to the Kramers theory, τp gives the time to overcome
an energy barrier, and the inverse of it is a chemical reaction
rate [23]. While most of the first passage time is spent by thewaiting time, the actual time required for a state transition canbe much smaller. Such a short timescale is characterized by τt[25–27,29,30]. For a nucleic acid folding, it was estimated tobe τt  10−5 s [31–33].
In this section, we obtain the analytical expressions of τpand τt in terms of the potential parameters A and F in Gr .
Within the decoupling approximation, Kd 2
A, we consider
the dynamics of θ by using Eq. (15). Then one can obtain bothτp and τt for a general potential [27].
A. Mean first passage time τp

The first passage time tp is a time for a reaction thatstarted from the initial value θ0 reaches the final value θf forthe first time. Notice that θ0 = 0 and θf = 2π in our model.
Mathematically, this is equivalent to consider a Brownian motion of a particle in a semi-infinite system with an absorbing
boundary condition at θ = θf , and to measure the time until
the particle is absorbed at θ = θf [see Fig. 7(a)]. However,
because the state θ = θ0 can be visited multiple times, mostof the first passage time is spent by a long waiting time. Since
tp is a stochastic quantity and has a broad distribution, it isuseful to consider its average value τp called the mean firstpassage time.
A formal derivation of τp is reviewed in Appendix B. For
an arbitrary periodic function Gp , τp can be given by
 2π
 2π
1
dx
dy
τp =
D(1 − e−2π f ) 00
× exp[gp (x) − gp (x − y) − f y],(22)

where gp = Gp /(kB T ) and f = F/(kB T ) are the dimensionless potential and nonequilibrium force, respectively.
Moreover, we have defined the diffusion constant D = MkB T ,where M is the Onsager coefficient used in Eqs. (6) and (15).
Note that D has the dimension of inverse time in the present
model.
The periodic potential Gp (θ ) in Eq. (3) should satisfy
Gp (θ + 2π ) = Gp (θ ) and its energy barrier is A. Among
various possibilities, the simplest form would be Gp (θ ) =−A cos θ that was used in our numerical simulation. However,to perform the integral in Eq. (22) analytically, we employ
here the following linear functions:

⎧ 
2
⎪
⎪A
θ −1
for 0  θ < π⎨πGp (θ ) =
,
(23)
⎪
2
⎪
⎩A − θ + 3
for π  θ < 2ππ

FIG. 7. (a) Schematic description of the first passage time tp . We
consider a semi-infinite system with an absorbing boundary condition at θ = 2π , and measure the time until a Brownian particle isabsorbed at θ = 2π . Here the state θ = 0 can be visited multiple
times, which gives rise to a long waiting time. (b) Schematic description of the first transition path time tt . We consider a finite systemwith absorbing boundary conditions at θ = 0 and 2π , and measure
the time until a Brownian particle is absorbed at θ = 2π . When theparticle is absorbed at θ = 0 [the red (gray) trajectory], such an event
is not counted.

as depicted in Fig. 2(a). With this periodic potential, one cananalytically obtain τp as
πf
1 + e−2π f
2π 2
8a2
τp =
−
D π 2 f 2 − 4a2(π 2 f 2 − 4a2 )2 1 − e−2π f+16a2 cosh(2a) e−π f,(π 2 f 2 − 4a2 )2 1 − e−2π f(24)where a = A/(kB T ) is the dimensionless energy barrier. Thenthe asymptotic expressions of τp are given as follows:
2π
2π
=
(a
1 or af ),DfMFπ (kB T )2 e2A/kB T
π e2a=(a  1 and fτp ≈2Da2 f
2MA2 F
τp ≈

(25)
1),
(26)

τp ≈

2 2a −π f

2A/kB T −πF/kB Te
π e e
π kB Te=Da2
MA2(a  f  1).
2(27)

Here we have recovered the dimension in the last expressions
for the clarity sake. Since Eq. (25) does not depend on thetemperature, thermal fluctuations are irrelevant in this limit.
This is not the case for Eqs. (26) and (27) which diverge when
the temperature vanishes due to the exponential factors.
In Fig. 8(a), we plot the scaled τp in Eq. (24) as a functionof f for a = 0.1 and 10. For a = 0.1 (black), the entire behavior is simply approximated by Eq. (25). For a = 10 (red062113-7

KENTO YASUDA AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 103, 062113 (2021)

FIG. 8. (a) The dimensionless mean first passage time τp as a function of the dimensionless nonequilibrium force f for a = 0.1 (black) and
10 (red or gray). The plot for a = 0.1 almost coincides with the asymptotic form in Eq. (25). The dotted and dashed lines are the asymptotic
expressions in Eqs. (26) and (27), respectively. (b) The dimensionless mean first transition path time τt as a function of the dimensionlessnonequilibrium force f for a = 0.1 (black) and 10 (red or gray). The dotted and dashed lines are the asymptotic expressions in Eqs. (35)and (36), respectively. The dash-dotted line is the asymptotic expression in Eq. (37). (c) The ratio τt /τp as a function of the dimensionlessnonequilibrium force f for a = 0.1 and 10. The asymptotic value of τt /τp approaches unity for large f , while it strongly depends on a forsmall f .
3
2π f
) + 3(2a)2 (1 − e2π f ),1 = (2a) (1 − eor gray), on the other hand, we have plotted Eqs. (26) (dottedline) and (27) (dashed line), which are in good agreement withthe full expression of τp .

22a π f− 3e2π f )2 = (2a) (−1 + 2e e

B. Mean first transition path time τt

The first transition path time tt is a time for a reactionthat started from the initial value θ0 reaches the final value
θf without returning to θ0 . Mathematically, this is equivalent
to consider a Brownian motion of a particle in a finite systemwith absorbing boundary conditions both at θ = θ0 and θ = θfand to measure the time until the particle is absorbed at θ = θf[the black trajectory in Fig. 7(b)]. When the particle is absorbed at θ = θ0 , such an event is not counted [the red (gray)trajectory in Fig. 7(b)]. Since tt is also a random quantity, we
consider its average τt called the mean first transition pathtime.
A formal derivation of τt is explained in Appendix C, andthe result is given by
 2π
 2π  x−1  2π1τt =dw exp[gr (w)]
dx
dy
dz
D 0
0
x
0
× exp[−gr (x) + gr (y) + gr (z)],(28)

where gr = Gr /(kB T ). Using Eqs. (3) and (23), we can analytically obtain τt as
2π 2

τt =

D[2a(1 − 2e2a eπ f + e2π f ) + π f (e2π f − 1)]×0+

1 (π f ) +23
2 (π f ) + 3 (π f ) +(π 2 f 2 − 4a2 )24 (π f )

4

,

(29)
where
4 2a π f
0 = −2(2a) e e+ (2a)3 (e−2a eπ f + 3e2a eπ f − 2 − 2e2π f ),

(30)

(31)

+ (2a)(−e−2a eπ f + 5e2a eπ f − 2e2π f − 2),

(32)

2π f
− 1) − (e2π f − 1),
3 = 2a(e

(33)

2π f
+ 1.
4 =e

(34)

Then the asymptotic expressions of τt are given as follows:
2π 2
2π 2
=
(a
1 and f
1),
(35)
3D
3MkB T
2π
2π=( f  1 and a
f ),
(36)
τt ≈
Df
MF
π2
π2
=
(a  1 and a  f ).
(37)
τt ≈
Da
MA
In the limit of Eq. (35), the transition process is dominated bythermal fluctuations. On the other hand, Eqs. (36) and (37) are
independent of the temperature, and hence the transitions occur deterministically. The scaling relation τt ∼ 1/a in Eq. (37)
was also obtained before for a quadratic potential [34].
In Fig. 8(b), we plot the scaled τt as a function of f for a =0.1 and 10. For a = 0.1 (black), τt is constant for f < 1 and it
decreases for f > 1. This behavior is in accordance with theasymptotic expressions in Eq. (35) (dotted line) and Eq. (36)(dashed line). For a = 10 (red or gray), on the other hand, τttakes a maximum value around f ≈ a. The dash-dotted line isthe asymptotic expression in Eq. (37).
In Fig. 8(c), we plot the ratio τt /τp as a function of f for
a = 0.1 and 10. For a = 0.1 (black), a power law dependence
is seen for f < 1, and τt /τp approaches unity for f > 1. For
a = 10 (red or gray), on the other hand, τt /τp is vanishingly
small for f < 1, and it grows exponentially for 1 < f < a.
It is worth mentioning here the characteristic difference
between τp and τt . In the limit of f → 0, τp diverges [see062113-8

τt ≈NONRECIPROCALITY OF A MICROMACHINE DRIVEN …

PHYSICAL REVIEW E 103, 062113 (2021)

Eqs. (25) and (26)] while τt remains finite [see Eqs. (35) and(36)]. This is because a nonequilibrium driving force is always
required for the net chemical reaction with a finite τp . On
the other hand, τt can be evaluated even in the equilibrium
situation.
VI. SUMMARY AND DISCUSSION

In this paper, we have discussed cyclic state transitionsof a micromachine driven by a catalytic chemical reaction.
We have proposed a minimum model of a generic micromachine and calculated the nonreciprocality in Eq. (1) to quantify
the functionality. Our model uses the reaction variable θ andthe state variables si which are coupled to each other throughthe mechanochemical coupling mechanism. The tilted periodic potential Gr (θ ) for catalytic reaction is characterized by
the energy barrier A and the nonequilibrium force F . Importantly, the state variables si are required to change periodicallyin time for a catalytic reaction.
To investigate our model analytically, we have employed
the decoupling approximation for the Onsager’s equations.
Furthermore, we have assumed that the reaction variable θobeys a deterministic stepwise dynamics characterized by themean first passage time τp and the mean first transition pathtime τt . Under these assumptions, we have first obtained thetime dependencies of the state variables s1 and s2 in termsof τp and τt [see Eqs. (17) and (18)]. We find that the magnitudes of s1 and s2 become smaller when τt /τp or γ τp isdecreased. Then we have analytically obtained the nonreciprocality R12 within the decoupling approximation [see Eq. (19)].
One of the important results is the asymptotic scaling relation
R12 ∼ (γ τt )2 for τt /τp
1 [see Eq. (21)]. Using Eq. (15) inthe small coupling limit, we have further obtained τp [seeEq. (24)] and τt [see Eq. (29)] in terms of the potential parameters A and F .
So far, the nonreciprocality R12 has been obtained in termsof τp and τt , whereas they are further expressed in termsof A and F . For a realistic chemical reaction such as ATPhydrolysis, τt /τp is typically small, and we expect that R12 canbe well approximated by Eq. (21). On the other hand, the limitof Eq. (37) is appropriate for τt when a = A/(kB T ) is large
enough. Using the corresponding asymptotic expressions, therelevant scaling for the nonreciprocality R12 turns out to be

 2

dγ
dμK 2
R12 ∼
sin(φ2 − φ1 ) ∼sin(φ2 − φ1 ).
Da
MA
(38)
In particular, the relation R12 ∼ 1/A2 implies that the higher
the energy barrier is, the smaller the nonreciprocality becomes. This scaling relation is another important result of thepresent model.
Next we discuss typical values of the model parameters.
From the experiment measuring the enhanced diffusion of a
motor protein, the energy barrier of ATP was estimated to
be A ∼ 10 kB T [35]. When a single ATP molecule is converted into ADP, the produced energy is roughly 20 kB T [6].
Then we estimate the nonequilibrium chemical force as F ∼20/(2π ) kB T ∼ 3 kB T (notice again that the dimension of F isenergy). Hence both A/kB T > 1 (a > 1) and A > F (a > f )are satisfied. Moreover, one can estimate from Eqs. (26) and
1, which justifies the assumption in(37) that τt /τp ∼ 10−6
Eq. (38).
Recent experiments reported the diffusion enhancement of
enzymes due to catalytic chemical reactions [8,11]. Whena self-propelled particle undergoes a rotational diffusion, its
translational diffusion coefficient increases by D = V 2 τrot ,where V is the propulsion velocity and τrot is the rotational diffusion time [36]. Since the propulsion velocity is proportional
to the nonreciprocality, V ∼ R, the diffusion enhancement canbe estimated as D ∼ R2 τrot ∼ τrot /A4 . So far, the relation
between the energy barrier and the functionality of a micromachine has not yet been investigated. We predict that thechange in the energy barrier can be reflected in the diffusion
enhancement of enzymes.
In the present study, we have mainly discussed the casewhen there are only two degrees of freedom (s1 and s2 ) ofa micromachine. Although this is a minimum and sufficientnumber to discuss the nonreciprocality R12 , one needs to take
into account a large number of state variables to describe
the dynamics of realistic enzymes. As mentioned in the introduction, the total functionality
of a micromachine can be
expressed by the weighted sum i, j Wi j Ri j , where Wi j is theweighting tensor that depends on the properties of a micromachine. Hence, it is necessary to evaluate the nonreciprocalities
Ri j for all the binary combinations of the state variables.
Although the estimation of the weight tensor Wi j is beyond
the scope of this work, such a study will be important in the
future.
In the present work, the nonreciprocality R12 has beenobtained only in the weak coupling limit, Kd 2
A. Theinvestigation of the opposite strong coupling limit is also left
as a future work such as performing more extended numerical simulations. It would be also interesting to see the casewhen the off-diagonal elements of the Onsager coefficient arenonzero, i.e., μi j = 0 for i = j.
ACKNOWLEDGMENTS

We thank I. Sou and Y. Hosaka for useful discussions.
K.Y. acknowledges support by a Grant-in-Aid for JSPS Fellows (Grant No. 18J21231) from the Japan Society for the
Promotion of Science (JSPS). S.K. acknowledges support
by Grants-in-Aid for Scientific Research (C) (Grants No.
18K03567 and No. 19K03765) 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, and Technology of Japan.
APPENDIX A: DERIVATION OF EQS. (17) AND (18)

In this Appendix, we show the derivation of Eqs. (17) and(18). In the absence of the noise terms, Eqs. (13) and (14) can
be formally solved as t
δi (t ) = −e−γ tdt  eγ t d cos[θ (t  ) + φi ]θ̇ (t  ),(A1)−∞

where we have assumed 0  t  τp and ignored a term that
depends on the initial condition. Using Eq. (16) for θ and the062113-9

KENTO YASUDA AND SHIGEYUKI KOMURAPHYSICAL REVIEW E 103, 062113 (2021)

condition θ (t + τp ) = θ (t ) + 2π , we obtaint
2π −γ t2πt (τt − t  )dt  eγ t cos+ φi
δi (t ) = − de
τt
τt
0


∞  −nτp +τt
2π 
 γ t
−dt e cos(t + nτp ) + φi ,τtn=1 −nτp(A2)where (t ) is the Heaviside step function. Changing the variable to t  = t  + nτp in the second integral, we obtaint
2π −γ t2πt  γ t
(τt − t  )δi (t ) = − dedt e cos
+ φi
τt
τt
0


 τt
∞

2πt 
−γ nτp γ t −e
dt e cos+ φi , (A3)τt0n=1which results in Eqs. (17) and (18) after the integration.
APPENDIX B: DERIVATION OF EQ. (22)In this Appendix, we show the derivation of the mean firstpassage time τp in Eq. (22) [23,37,38]. For this purpose, we
consider a conditional probability distribution P(θ , t|θ0 ) for
which θ (0) = θ0 is imposed as the initial condition. ThenP(θ , t|θ0 ) satisfies the following Fokker-Planck equation:∂t P(θ , t|θ0 ) = L(θ )P(θ , t|θ0 ),(B1)

L(θ ) = D∂θ exp[−gr (θ )]∂θ exp[gr (θ )].

(B2)

Similarly, P(θ , t|θ0 ) also satisfies the following backward
Fokker-Planck equation:∂t P(θ , t|θ0 ) = L† (θ0 )P(θ , t|θ0 ),L† (θ0 ) = D exp[gr (θ0 )]∂θ0 exp[−gr (θ0 )]∂θ0 .

(B3)
(B4)

We employ the reflective boundary condition at θ → −∞and the absorbing boundary condition at θ = 2π . Then the
total probability distribution decays due to the latter boundary
condition. Here we introduce the survival probability definedas 2πS(t, θ0 ) =dθ P(θ , t|θ0 ).
(B5)

Next, one can show from Eqs. (B3), (B5), (B6), and (B8)that
−1 = L† (θ0 )τp (θ0 ),where we have used the conditions limt→0 tKp (t ) = 0 andlimt→∞ tKp (t ) = 0. Using the reflective boundary condition(∂θ0 τp = 0 at θ0 → −∞) and the absorbing boundary conditions (τp = 0 at θ0 = 2π ), one can solve the above equation to
obtain [39]
 x

1 2π
τp (θ0 ) =dx exp[gr (x)]dY exp[−gr (Y )].
D θ0
−∞
(B10)
Setting θ0 = 0 and using Eq. (3), we get
 ∞

1 2π
dx
dy exp[gp (x) − gp (x − y ) − f y ],
τp =
D 0
0
(B11)
where y = x − Y . Because the periodicity of Gp is 2π , weobtain 2π
∞1  −2nπF 2π
τp =
e
dx
dy
D n=0
0
0
× exp[gp (x) − gp (x − y) − f y],

dS
.
(B6)
dt
From the condition S(0) = 1, the following normalization
condition holds:
 ∞dt Kp (t, θ0 ) = 1.
(B7)
Kp (t, θ0 ) = −0

The mean first passage time τp is defined as the first momentof the distribution function ∞τp (θ0 ) =dt tKp (t, θ0 ).
(B8)

(B12)

where y = y − 2nπ . Since F > 0, we can easily evaluate theinfinite series and obtain Eq. (22).
APPENDIX C: DERIVATION OF EQ. (28)In this Appendix, we show the derivation of the mean
first transition path time τt in Eq. (28) [25–27]. In this case,the absorbing boundary condition is imposed both at θ = 0
and θ = 2π . Hence, unlike τp , the probability is absorbed
from both of the boundaries, although the probability of being
absorbed at θ = 2π determines τt .
The distribution function of the first transition path time isgiven by1
J (θ = 2π , t|θ0 ),N
where J (θ , t|θ0 ) is a probability flux
Kt (t, θ0 ) =J (θ , t|θ0 ) = −D exp[−gr (θ )]−∞Then the distribution function of the first passage time is given
by

(B9)

(C1)

∂
[exp[gr (θ )]P(θ , t|θ0 )].
∂θ
(C2)

In Eq. (C1), N is the normalization constant that is fixed by
the condition ∞dt Kt (t, θ0 ) = 1.
(C3)
0

The mean first transition path time τt is defined as the first
moment of the distribution function
 ∞dt tKt (t, θ0 ).
(C4)
τt (θ0 ) =0The backward Fokker-Planck equation for the probability
flux J is given by0062113-10∂t J (θ , t|θ0 ) = L† (θ0 )J (θ , t|θ0 ).

(C5)

NONRECIPROCALITY OF A MICROMACHINE DRIVEN …

PHYSICAL REVIEW E 103, 062113 (2021)

From Eqs. (C1), (C3), and (C5), one can show that
0 = L (θ0 )N (θ0 ),†(C6)

where we have used the conditions J (2π , 0|θ0 ) = 0 andJ (2π , ∞|θ0 ) = 0. Solving this equation with the boundary
conditions N (0) = 0 and N (2π ) = 1, we obtain
 2π
 θ0
−1
N (θ0 ) =dy0 exp[gr (y0 )]dy0 exp[gr (y0 )] .
0of the absorbing boundary conditions ψ (0) = 0 and ψ (2π ) =0, we can solve the above equation to obtain1
ψ (θ0 ) =
D

 2π
0


 θ0
× [1 − N (θ0 )]dx exp[−gr (x)]N 2 (x)

0

(C7)

+ N (θ0 )From Eqs. (C1), (C4), and (C5), we obtain−N (θ0 ) = L† (θ0 )ψ (θ0 ),−1

dw exp[gr (w)]

 2π
θ0

0

dx exp[−gr (x)][1 − N (x)]N (x) .
(C9)

(C8)

where ψ (θ0 ) = τt (θ0 )N (θ0 ) and we have used the conditions
limt→0 tKt (t, θ0 ) = 0 and limt→∞ tKt (t, θ0 ) = 0. With the use

Since τt = limθ0 →0 [ψ (θ0 )/N (θ0 )], only the second term remains and we obtain Eq. (28).

[1] S. Toyabe and M. Sano, J. Phys. Soc. Jpn. 84, 102001 (2015).
[2] A. I. Brown and D. A. Sivak, Chem. Rev. 120, 434 (2020).
[3] Y. Togashi, T. Yanagida, and A. S. Mikhailov, PLoS Comput.
Biol. 6, e1000814 (2010).
[4] M. L. Mugnai, C. Hyeon, M. Hinczewski, and D. Thirumalai,Rev. Mod. Phys. 92, 025001 (2020).
[5] T. Harada and S.-i. Sasa, Phys. Rev. Lett. 95, 130602 (2005).
[6] S. Toyabe, T. Okamoto, T. Watanabe-Nakayama, H. Taketani,
S. Kudo, and E. Muneyuki, Phys. Rev. Lett. 104, 198103(2010).
[7] T. Ariga, M. Tomishige, and D. Mizuno, Phys. Rev. Lett. 121,
218101 (2018).
[8] H. S. Muddana, S. Sengupta, T. E. Mallouk, A. Sen, and P. J.
Butler, J. Am. Chem. Soc. 132, 2110 (2010).
[9] R. Golestanian, Phys. Rev. Lett. 115, 108102 (2015).
[10] S. Ghosh, A. Somasundar, and A. Sen, Annu. Rev. Condens.
Matter Phys. 12, 177 (2021).
[11] A.-Y. Jee, S. Dutta, Y.-K. Cho, T. Tlusty, and S. Granick,Proc. Natl. Acad. Sci. USA 115, 14 (2018).
[12] T. Sakaue, R. Kapral, and A. S. Mikhailov, Eur. Phys. J. B 75,
381 (2010).
[13] J. E. Avron, O. Kenneth, and D. H. Oaknin, New J. Phys. 7, 234(2005).
[14] R. Golestanian and A. Ajdari, Phys. Rev. E 77, 036308 (2008).
[15] E. Lauga, The Fluid Dynamics of Cell Motility (Cambridge
University Press, Cambridge, UK, 2020).
[16] G. Grosjean, M. Hubert, G. Lagubeau, and N. Vandewalle,
Phys. Rev. E 94, 021101(R) (2016).
[17] I. Sou, Y. Hosaka, K. Yasuda, and S. Komura, Phys. Rev. E 100,
022607 (2019).
[18] I. Sou, Y. Hosaka, K. Yasuda, and S. Komura, Physica A(Amsterdam, Neth.) 562, 125277 (2021).
[19] M. Leoni and P. Sens, Phys. Rev. Lett. 118, 228101 (2017).
[20] M. Tarama and R. Yamamoto, J. Phys. Soc. Jpn. 87, 044803(2018).

[21] A. Shapere and F. Wilczek, J. Fluid Mech. 198, 557 (1989).
[22] K. Dill and S. Bromberg, Molecular Driving Forces: Statistical Thermodynamics in Biology, Chemistry, Physics, andNanoscience (Garland Science, London, 2010).
[23] P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62,
251 (1990).
[24] M. Doi, Soft Matter Physics (Oxford University, Oxford, UK,2013).
[25] G. Hummer, J. Chem. Phys. 120, 516 (2004).
[26] S. Chaudhury and D. E. Makarov, J. Chem. Phys. 133, 034118(2010).
[27] W. K. Kim and R. R. Netz, J. Chem. Phys. 143, 224108 (2015).
[28] A. Ghanta, J. C. Neu, and S. Teitsworth, Phys. Rev. E 95,
032128 (2017).
[29] M. Laleman, E. Carlon, and H. Orland, J. Chem. Phys. 147,
214103 (2017).
[30] M. Caraglio, T. Sakaue, and E. Carlon, Phys. Chem. Chem.
Phys. 22, 3512 (2020).
[31] K. Neupane, D. B. Ritchie, H. Yu, D. A. N. Foster, F. Wang, and
M. T. Woodside, Phys. Rev. Lett. 109, 068102 (2012).
[32] H. S. Chung and W. A. Eaton, Nature (London) 502, 685(2013).
[33] K. Neupane, F. Wang, and M. T. Woodside, Proc. Natl. Acad.
Sci. USA 114, 1329 (2017).
[34] H. S. Chung, J. M. Louis, and W. A. Eaton, Proc. Natl. Acad.
Sci. USA 106, 11837 (2009).
[35] R. Hayashi, K. Sasaki, S. Nakamura, S. Kudo, Y. Inoue, H. Noji,and K. Hayashi, Phys. Rev. Lett. 114, 248101 (2015).
[36] H.-R. Jiang, N. Yoshinaga, and M. Sano, Phys. Rev. Lett. 105,
268302 (2010).
[37] P. Reimann, C. Van den Broeck, H. Linke, P. Hänggi, J. M.
Rubi, and A. Pérez-Madrid, Phys. Rev. Lett. 87, 010602 (2001).
[38] K. Hayashi and S.-i. Sasa, Phys. Rev. E 69, 066119 (2004).
[39] N. Goel and N. Richter-Dyn, Stochastic Models in Biology(Academic Press, New York, 1974).

062113-11