Odd elasticity of a catalytic micromachine

We perform numerical simulations of a model micromachine driven by catalytic chemical reactions. Our model includes a mechano-chemical coupling between the structural variables and the nonequilibrium variable describing the catalytic reactions. The time-correlation functions of the structural variables are calculated and further analyzed in terms of odd Langevin dynamics. We obtain the effective odd elastic constant that manifests the broken time-reversal symmetry of a catalytic micromachine. Within the simulation, we separately estimate the quantity called nonreciprocality and show that its behavior is similar to that of the odd elasticity.Our approach suggests a new method to extract the nonequilibrium properties of a micromachine only by measuring its structural dynamics.


I. INTRODUCTION
In recent years, the physics of micromachines such as bacteria, motor proteins, and artificial molecular machines has been intensively studied [1][2][3].A micromachine can be defined as a small object that extracts energy from chemical substances in the system and further exhibits mechanical functions [4,5].The interplay between the structural dynamics of such a small object and associated chemical reaction is crucial for biological functions of a micromachine [6][7][8].Owing to the developments in nonequilibrium statistical mechanics and experimental techniques, various research has been conducted to reveal the energetics of a single micromachine [9][10][11][12].
In our previous work, we proposed a model that describes cyclic state transitions of a micromachine driven by catalytic chemical reactions [13].As we shall explain later, we considered a mechano-chemical coupling of the variables representing the chemical reaction and the internal structural state of a micromachine.To estimate the functionality of a micromachine, we focused on the physical quantity called "nonreciprocality" that represents the area enclosed by a trajectory in the conformational space.Such a quantity directly determines the average velocity of a three-sphere microswimmer [14][15][16] or the crawling speed of a cell on a substrate [17,18].
Recently, Scheibner et al. introduced the concept of odd elasticity that is useful to characterize nonequilibrium active systems [19,20].Odd elasticity arises from antisymmetric (odd) components of the elastic modulus tensor that violates the energy conservation law and thus can exist only in active materials [21,22] or biological systems [23].Recently, some of the present authors proposed a thermally driven microswimmer with odd elastic- * komura@wiucas.ac.cn ity and demonstrated that it can exhibit a directional locomotion [24,25].The concept of odd elasticity was also applied to Purcell's three-link swimmer model [26] and to robotic systems [27].We also showed that antisymmetric parts of the time-correlation functions in odd Langevin systems are proportional to the odd elasticity [28,29].Whilst many theoretical models are proposed to explain microscopic physical origins of odd elasticity, quantitative measurements of the odd elastic modulus are still lacking in these stochastic and nonequilibrium systems.
In this paper, we perform extensive numerical simulations of a model micromachine driven by catalytic chemical reactions [13].We calculate the time-correlation functions of the structural variables and analyze them in terms of effective even and odd elasticity of coupled Langevin equations [28].Importantly, the presence of odd elasticity reflects the broken time-reversal symmetry of a catalytic micromachine.We also show that the nonreciprocality is proportional to the odd elasticity in the current stochastic model.Our approach provides us with a new method to extract the nonequilibrium properties of a micromachine only by measuring its structural dynamics without invoking the dynamics of the chemical reaction [30].
The current work is based on our previous two papers in Refs.[13] and [28], and further completes our general view on stochastic micromachines.In these works, we explicitly discuss the effects of thermal fluctuations which provide driving forces for micromachines.One of the results specific to our stochastic system is that the average work is proportional to the square of the odd elastic constant (see Sec. III D) [24].It should be also emphasized that the concept of odd elasticity is not limited to elastic materials.Although we focus on catalytic micromachines, most of the results also apply to general odd stochastic systems such as the odd microswimmer [28,29].In our previous paper in Ref. [13], we solved our model almost analytically by using several approxi-mations.In this work, on the other hand, we numerically solve the same equations in the presence of noise without using any crude approximations.Such a study is necessary for the complete understanding of the proposed model.
In Sec.II, we briefly review our model of a micromachine driven by a catalytic chemical reaction and explain the numerical method [13].In Sec.III, we show the results of the numerical simulations and discuss the properties of the time-correlation functions [28].Then we estimate effective even and odd elastic constants by using the result of the overdamped odd Langevin equations.A summary of our work and some discussion are given in Sec.IV.

A. Model
We first explain our model of a catalytic micromachine which was introduced in Ref. [13] and is schematically shown in Fig. 1(a).Consider a system which contains one enzyme molecule (E) that acts as a micromachine, substrate molecules (S), and product molecules (P).The enzyme molecule plays the role of a catalyst and the corresponding chemical reaction is written as [31] S + E ⇄ ES → P + E where ES indicates a complex molecule.In our model, we introduce a dimensionless reaction variable θ(t) to quantify the extent of the catalytic reaction.The reaction variable θ is a continuous number and increases 2π for each reaction.
According to the Kramers theory, the free energy G r describing a chemical reaction is given by a tilted periodic potential [32] as schematically represented in Fig. 1(b).The first term is a periodic potential with a period of 2π and A is the energy barrier.This is because θ increases by 2π for one cycle of chemical reaction and should experience the same potential.On the other hand, F in Eq. ( 2) represents the chemical potential difference that drives catalytic reaction.For example, F represents the chemical potential change before and after ATP hydrolysis.Next, we introduce the dimensionless state variables s i (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.To introduce the mechano-chemical coupling mechanism, we assume that each state variable s i experiences a harmonic potential, (C i /2)[s i − ℓ i (θ)] 2 , where C i is the coupling parameter and ℓ i (θ) is the dimensionless natural state that depends on the reaction variable θ.We consider that the natural state ℓ i (θ) micromachine substrate < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 r q h / K Schematic picture of a micromachine characterized by the conformational state variables s1, s2, and s3.Moreover, the adhesion between the domains and the substrate is described by the variables s4 and s5.(b) The tilted periodic potential Gr(θ), where θ is the catalytic reaction variable.As expressed in Eq. ( 2), A is the energy barrier and F corresponds to the nonequilibrium force.A possible trajectory of θ is schematically shown by the blue arrow.The value of θ fluctuates around the minimum of the potential and the transition to the next minimum takes place occasionally.Reprinted figure with permission from Ref. [13].Copyright (2021) by the American Physical Society.
changes periodically and assume the simplest periodic form ℓ i (θ) = sin(θ + ϕ i ), where ϕ i is the constant phase difference relative to the reaction phase θ [33,34].Under these assumptions, we consider the following mechanochemical coupling energy G c between θ and s i : Then the total free energy G t in our model is simply given by For the time evolution of θ and s i , we employ the Onsager's phenomenological equations in the presence of noises [35,36] where dot indicates the time derivative.In the above equations, the time-evolutions of the variables are proportional to the respective thermodynamic forces, and they relax to the thermal equilibrium state for which the total free energy is minimized.Here, µ θ and µ ij are the mobility coefficients, and ξ θ and ξ i represent thermal fluctuations which satisfy the following fluctuationdissipation theorem [35,36] where k B is the Boltzmann constant and T is the temperature.
Although our model is general and includes many degrees of freedom, we make several simplifications to capture the physical insight of a catalytic micromachine.First, we only consider two degrees of freedom, i.e., s 1 and s 2 .Second, the mobility coefficients µ ij (i, j = 1, 2) is assumed to have the form µ ij = µ s δ ij , where δ ij is the Kronecker delta.Third, the coupling free energy is symmetric between the two degrees of freedom, i.e., C 1 = C 2 = C.In this work, we investigate the case when ϕ 1 = 0 because different choices of ϕ 1 give qualitatively the same ϕ 2 -dependence.Then Eqs. ( 5) and ( 6) reduce to the following simplified set of equations [13] Notice that the above equations are difficult to solve analytically because they are highly nonlinear although some attempts were made in Ref. [13].

B. Numerical method
In the above model, the variables θ and s i are dimensionless, whereas A, F , and C have the dimension of energy.The latter quantities are all scaled by the thermal energy k B T and the corresponding dimensionless quantities are defined by Â = A/(k B T ), F = F/(k B T ), and Ĉ = C/(k B T ), respectively.The dimensionless time is defined by t = 2µ θ k B T t and the ratio between the two mobilities is fixed to µ s /µ θ = 1 in our simulation.
We numerically solve the coupled stochastic differential equations in Eqs. ( 10)-( 12) by integrating the following dimensionless quantities where dW θ , dW 1 , dW 2 ∼ √ dt are the increment of Wiener processes.
When we numerically calculate the statistical average ⟨X⟩ of a quantity X, we use the long-time average.This is justified because we have checked that the ensemble average and the time average give the same result in our simulation after the system reaches the steady state.

A. Time evolution
As an example of the simulation result, we show in Fig. 2(a) the time evolution of θ (top), s 1 (middle), and s 2 (bottom) for certain set of parameters.When the catalytic reaction occurs, we see that the reaction variable θ increases in a stepwise manner by 2π, whereas the state variables s 1 and s 2 undergo almost random fluctuations.Corresponding to the occasional increase of θ, both s 1 and s 2 tend to show peaks as indicated by the vertical red lines.Although the result in Fig. 2(a) demonstrates that the present model micromachine is indeed driven by thermal fluctuations, it is difficult to isolate the peaks of s i because they are almost comparable to the background fluctuations.Moreover, the state transitions become very rare when Â ≫ 1.
In Fig. 2(b), we plot a typical trajectory of s 1 and s 2 for one typical cycle of chemical reaction.While θ fluctuates around the local minimum, s 1 and s 2 also fluctuate (black trajectory) around the average structure specified by the white circle.When θ increases by 2π, the trajectory of s 1 and s 2 follows a directional closed loop (red trajectory).The direction is clockwise when 0 < ϕ 2 < π (as shown by the red arrow) or counterclockwise when π < ϕ 2 < 2π.The area enclosed by the trajectory can be obtained by the following quantity called "nonreciprocality" [13] where the integral is taken over one cycle.The statistical property of R 12 in our stochastic simulation will be discussed later in Sec.III D. (We avoid to use the terminology "nonreciprocity" that has broader and general meanings [37][38][39].)

B. Time-correlation functions
From the time series of s 1 and s 2 in Fig. 2(a), we calculate the structural time-correlation functions of a catalytic micromachine.We first define the time-dependent fluctuation by x i (t) = s i (t) − ⟨s i ⟩, where ⟨s i ⟩ is the average of s i .Then the matrix of time-correlation functions ϕ ij (t) is defined as In the above, ϕ ij (t) is decomposed into the symmetric and antisymmetric parts which satisfy ϕ S ij (t) = ϕ S ji (t) and , respectively [28].For our later purpose, we also define the equal-time-correlation functions with a bar as φij = ϕ ij (0).
When the system is in equilibrium for which the timereversal symmetry holds, the correlation functions must satisfy the reciprocal relation ϕ eq ij (t) = ϕ eq ji (t) for i ̸ = j [35,36].This is equivalent to saying that the antisymmetric part of the correlation function ϕ A ij (t) should vanish in equilibrium [28].In nonequilibrium situations, however, the antisymmetric parts can exist because the time-reversal symmetry can be generally violated [40][41][42].
In Fig. 3(a), we plot both the self-correlation functions ϕ 11 (t) (black) and ϕ 22 (t) (red) as a function of the dimensionless time t.The overall behavior of the selfcorrelation functions is nearly described by an exponential function with a slight modification.Although the equal-time-correlation functions φ11 and φ22 take different values, the slopes of the initial decay are the same for these self-correlation functions.Notice that the initial decay of the correlation function reflects the mobility of the system.In Fig. 3(b), we plot both the symmetric part of the cross correlation function ϕ S 12 (t) = [ϕ 12 (t) + ϕ 21 (t)]/2 (black) and the antisymmetric part of the cross correlation function ϕ A 12 (t) = [ϕ 12 (t) − ϕ 21 (t)]/2 (red) as a function of t.We clearly see that ϕ A 12 (t) is nonzero, namely, ϕ 12 (t) ̸ = ϕ 21 (t) for 0 < t < 2. This means that the time-reversal symmetry is explicitly broken and hence the antisymmetric part of the cross-correlation function exists in the present catalytic system.As we discuss below, the antisymmetric parts of the correlation functions can be quantitatively characterized by effective odd elasticity [28].

C. Effective even and odd elastic constants
The short-time behaviors of the obtained correlation functions have the following characteristic features; (i) the initial slopes of the self-correlation functions ϕ 11 (t) and ϕ 22 (t) are the same, (ii) the initial slope of the symmetric part of the cross correlation function ϕ S 12 (t) vanishes, i.e., φS 12 (0) = 0, (iii) the initial value of the antisymmetric part of the cross correlation function φA 12 vanishes, i.e., φ12 = φ21 .
Under these conditions, the obtained correlation functions can be interpreted in terms of the following coupled odd Langevin equations: where µ is the effective mobility and ζ α (α = 1, 2) is Gaussian white noise that satisfies the fluctuation dissipation relation.More general situations and the statistical property of the noises were discussed in Ref. [28].
The effective elastic constant matrix K αβ can be written in terms of even and odd elastic constants as follows [see
In Appendix A, we show that the short-time behaviors of the symmetric and antisymmetric parts of the correlation functions obtained from Eqs. ( 18) and ( 19) are given by [see Eqs.(A11) and (A12)] where tr[K e ] is the trace of the matrix K e αβ and ϵ αβ is the 2D Levi-Civita tensor.The equal-time-correlation function φαβ in Eq. ( 21) is now given by [see Eq. (A10)] where det[K e ] is the determinant of K e αβ , ((K e ) −1 ) αβ is the inverse of K e αβ , and Using these expressions, we extract the effective even and odd elastic constants of the catalytic micromachine.We first checked that µ/µ θ ≈ 1.0 holds in the simulation, which is consistent with the choice of the parameter µ s /µ θ = 1.Hence we can identify µ with µ s and obtain the effective elastic constants.
In Fig. 4(a), we plot the three dimensionless even elastic constants Ke αβ = K e αβ /(k B T ) as a function of the phase difference ϕ 2 .We see that they are periodic function of ϕ 2 .The period of K e 11 and K e 22 is π, whereas that of K e 12 is 2π.In Fig. 4(b), on the other hand, we plot ko = k o /(k B T ) as a function of ϕ 2 .Notice that k o can take negative values and roughly approximated as k o ∼ − sin ϕ 2 .The present analysis implies that the effective elasticity (both even and odd) can be obtained only by measuring the structural dynamics without knowing any detailed dynamics of the chemical reaction variable θ.
In Figs.5(a) and (b), on the other hand, we plot Ke αβ and ko , respectively, as a function of the coupling parameter Ĉ while the phase difference is fixed to ϕ 2 = π/2.We see in Fig. 5(a) that the even elastic constants K e αβ increase with C.Although the sign of the odd elastic constant k o is negative in Fig. 5(b), its magnitude also increases when C is made larger.

D. Nonreciprocality and power efficiency
To characterize the nonequilibrium degree of a micromachine, it is useful to consider the nonreciprocality defined in Eq. ( 16).For the current stochastic model, we first discuss the average number of reaction cycle ⟨N ⟩ per unit time that is obtained by dividing the total number of cycles by the total simulation time.In Fig. 6 plot ⟨N ⟩ as a function of the phase difference ϕ 2 for two combinations of A and F values.We find that ⟨N ⟩ is a periodic function of ϕ 2 and the period is π.
Next, we calculate the statistical average of nonreciprocality R 12 over many cycles of chemical reaction.To do so, we accumulate the quantity ṡ1 s 2 over long time and divide it by the total number of cycles during the simulation to obtain ⟨R 12 ⟩ per one cycle which is plotted in Fig. 7 as a function of ϕ 2 .We see that the behavior of ⟨R 12 ⟩ is similar to that of k o in Fig. 4(b) and can be approximated as ⟨R 12 ⟩ ∼ sin ϕ 2 .This result is reasonable because it was shown before that ⟨R 12 ⟩ is given by the time derivative of Eq. ( 22) and can be written as [26] where ⟨N ⟩ −1 is the average period of the reaction cycle.This expression confirms that ⟨R 12 ⟩ is proportional to k o .Hence the average nonreciprocality ⟨R 12 ⟩ becomes nonzero only if there is finite odd elasticity.We note again that such dynamics occur when the elastic constants are nonreciprocal, i.e., K 12 ̸ = K 21 in Eq. (20).The presence of odd elasticity indicates that mechanical work is done by the micromachine through its structural change and the average work per unit cycle is given by [19] In Fig. 8, we plot the average power efficiency defined by η = ⟨W ⟩/(2πF ) as a function of the phase difference ϕ 2 by using the simulation result.Notice that the energy supplied during one cycle of the catalytic reaction is 2πF .The power efficiency vanishes when ϕ 2 = 0, π, and 2π for which the deformation of the micromachine is reciprocal.This is also reasonable because if we substitute Eq. ( 24) into Eq.( 25) the average work should approximately behave as [24] ⟨W We also find that the power efficiency increases when the activation energy A and the nonequilibrium driving force F become smaller.Hence the micromachine can have a higher η value when thermal fluctuations are larger compared with the potential barrier.

IV. SUMMARY AND DISCUSSION
In this paper, we have performed numerical simulations of a model micromachine driven by catalytic chemical reactions as previously proposed by some of the present authors [13].We have obtained the structural time-correlation functions and also extracted their symmetric and antisymmetric parts (see Fig. 3).These results were further analyzed in terms of odd Langevin dynamics and we have obtained the effective even and odd elastic constants of the catalytic system (see Figs. 4 and  5).The presence of the odd elasticity demonstrates the broken time-reversal symmetry and we have confirmed that the nonreciprocality ⟨R 12 ⟩ is proportional to the odd elastic constant k o (see Fig. 7).We have also calculated the power efficiency of a micromachine and found that it increases when the activation energy A becomes smaller (see Fig. 8).
Let us summarize here our perspectives on stochastic micromachines, which have been discussed in our previous papers [13,28] as well as in this work that completes our view.First, the concept of odd elasticity is not necessarily limited to elastic materials as in the original work [19] but can also be applied to general nonreciprocal systems with odd interactions.Second, we have explicitly discussed the role of thermal fluctuations in catalytic micromachines and showed that the average work is proportional to the square of the odd elastic constant [see Eq. ( 26)].Third, we have suggested how to extract nonequilibrium properties of a micromachine only by measuring its structural dynamics even when the chemical reaction variable is a hidden nonequilibrium variable.
Recently, we have derived dynamical equations for a nonequilibrium active system with odd elasticity by using Onsager's variational principle [30].We showed that the elimination of an extra variable that is coupled to the nonequilibrium driving force leads to the nonreciprocal set of equations for the structural coordinates [43].The obtained nonreciprocal equations manifest the physical origin of the odd elastic constants that are proportional to the nonequilibrium force and the friction coefficients [19].In the present work, the reaction variable θ(t) is driven by the nonequilibrium force F and coupled to the structural variables s i .Similar to Ref. [30], the elimination of the nonequilibrium variable results in the nonreciprocal deformation of a micromachine characterized by R 12 in Eq. ( 16).This explains how a micromachine converts chemical energy to mechanical work and such a mechanism can be common for proteins, enzymes, and even microswimmers [24,25].
We have analyzed the time-correlation functions of a catalytic micromachine in terms of overdamped Langevin equations with odd elasticity as shown in Eq. (A1).In Appendix A, we have assumed that the corresponding mobility matrix M αβ is symmetric, i.e., M αβ = M βα , according to Onsager's reciprocal relation [36].If we remove this assumption and if M αβ is not symmetric, one also needs to consider the odd viscosity [44].Odd viscosity accounts for the fluid flow perpendicular to the velocity gradient and does not contribute to energy dissipation [45][46][47][48].Hence, odd viscosity needs to be taken into account when the surrounding environment of a micromachine is in out-of-equilibrium situations.In such cases, however, we cannot rely on the fluctuation dissipation relation that we have used in the derivation of the correlation functions [28].
Finally, we discuss some typical values of the parameters used in our numerical simulations and evaluate the effective odd elastic constant.According to the experiments on F 1 -ATPase, the activation energy and the chemical potential difference can be roughly estimated as A ≈ 10 k B T [11] and F ≈ 3 k B T [10].Using the relation for work per cycle ⟨W ⟩ = 2πF η as in the previous section, we estimate ⟨W ⟩ ≈ 2 k B T when η ≈ 0.1 (see Fig. 8).We shall further identify this work as ⟨W ⟩ ∼ κ o d 2 , where κ o is the odd elastic constant in unit of J/m 2 (k o in this paper has the dimension of energy) and d ≈ 10 −8 m is the typical enzyme size [28].Then one can roughly estimate as κ o ≈ 10 −4 J/m 2 , which is also consistent with our previous estimate for a kinesin molecule [28].(see Eq. (C6) in Ref. [28]).

FIG. 2 .
FIG. 2. (Color online) (a) Time evolution of θ, s1, and s2 as a function of dimensionless time t = 2µ θ kBT t when Â = 20, F = 16, Ĉ = 20, ϕ2 = π/2, and µs/µ θ = 1.The red vertical lines indicate the moments when the catalytic chemical reactions take place.(b) The trajectory of s1 and s2 corresponding to one cycle of chemical reaction that occurs in (a).When θ fluctuates around the local minimum of the potential energy Gr(θ), s1 and s2 also fluctuate (black trajectory for 115 ≤ t < 216.5) around the average position (white circle).When the chemical reaction takes place and θ changes by 2π, the trajectory forms a closed loop (red trajectory for 112.5 ≤ t < 115) in a clockwise manner as shown by the red circular arrow.

12 FIG. 3 .
FIG. 3. (Color online) (a) The self-correlation functions ϕ11(t) (black) and ϕ22(t) (red) as a function of the dimensionless time t.(b) The symmetric part of the cross correlation function ϕ S 12 (t) = [ϕ12(t) + ϕ21(t)]/2 (black) and the antisymmetric part of the cross correlation function ϕ A 12 (t) = [ϕ12(t) − ϕ21(t)]/2 (red) as a function of the dimensionless time t.For both (a) and (b), the parameters are the same as in Fig. 2.

FIG. 5 .
FIG. 5. (Color online) (a) The dimensionless effective even elastic constants Ke 11 (black), Ke 22 (red), and Ke 12 (blue) as a function of the dimensionless coupling parameter Ĉ = C/(kBT ).(b) The dimensionless effective odd elastic constant ko as a function of Ĉ.For both (a) and (b), the other parameters are the same as in Fig. 2.

16 FIG. 6 .
FIG. 6. (Color online) The average number of reaction cycle per unit time ⟨N ⟩ as a function of the phase difference ϕ2 when Â = 20 and F = 16 (black) and when Â = 10 and F = 8 (red).The other parameters are the same as in Fig. 2.
edges the JSPS, KAKENHI for Transformative Research Areas A (No. 21H05309) and the Japan Science and Technology Agency (JST), PRESTO Grant (No. JP-MJPR1921).K.Y. and K.I. were supported by the Research Institute for Mathematical Sciences, an International Joint Usage/Research Center located in Ky-oto University.S.K. acknowledges the support by the National Natural Science Foundation of China (Nos.12274098 and 12250710127) and the startup grant of Wenzhou Institute, University of Chinese Academy of Sciences (No. WIUCASQD2021041).A.K. and K.Y. contributed equally to this work.