Mechanochemical enzymes and protein machines as hydrodynamic force dipoles: The active dimer model

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 non-equilibrium fluctuating flows collectively generated by the enzymes, mixing in the solution and diffusion of passive particles within it are expected to get enhanced. Here, we investigate the intensity and statistical properties of such force dipoles in the minimal active dimer model of a mechanochemical enzyme. In the framework of this model, novel estimates for hydrodynamic collective effects in solution and in lipid bilayers under rapid rotational diffusion are derived, and available experimental and computational data is examined.


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 the ligand-bound state, are often accompanied by conformational transitions in them. Thus, these macromolecules would repeatedly change their shapes in each next turnover cycle. The primary role of 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, such motions are moreover employed to bring about the required machine function or to generate work.
Since enzymes are in solution, their active conformational changes are accompanied by flows in the fluid around them. Such non-equilibrium flows can affect internal mechanical motions in the enzymes and also influence translational and rotational diffusion of such proteins, as demonstrated by MD simulations for a model protein 1 and adenylate kinase. 2 It has been discussed whether hydrodynamic self-propulsion of enzymes could furthermore occur, in the models where either instantaneous transitions 3,4 or ligand-induced continuous conformational motions take place. [5][6][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 or transporters. Essentially, they represent protein machines powered by ATP hydrolysis or other catalytic reactions in them. Within each operation cycle, the shapes of their membrane domains typically change, 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 fluctuating flows 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 and diffusion enhancement should therefore arise. Additionally, chemotaxislike effects in the presence of spatial gradients in the 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 self-propulsion for proteins, which is predicted to be weak. [5][6][7] Following the original publication, 12 extensive further research has been performed. [13][14][15][16][17][18][19][20] The effects of rotational diffusion and of possible nematic ordering for enzymes were considered, 14 the phenomena in biomembranes were extensively analyzed, 15,16 and the 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 and principal theoretical predictions could thus be verified. 20 At low Reynolds numbers, the flow distribution produced by an object, changing the shape due to internal forces within it, can be characterized in the far field as that corresponding to a hydrodynamical force dipole. If the time-dependent stochastic force dipole of an enzyme is known, the collective hydrodynamic ef-fects in solution of such enzymes are predicted by the mean-field theory. 12 The difficulty, 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 effects 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, the active dimer model is formulated. The active dimer represents a minimal model where ligand-induced mechanochemical motions are reproduced. 12,[20][21][22] After presenting the model, we undertake an approximate analytical investigation of statistical properties of the force dipoles corresponding to active dimers in subsection 2.2, followed by a numerical study in subsection 2.3. Quantitative estimates for the intensity of hydrodynamical force dipoles in real enzymes are obtained in subsection 2.4. Section 3 corresponds to the second part. Based on the active dimer results, we obtain in subsection 3.1 more precise analytical and numerical estimates for the maximal diffusion enhancement for passive particles in solutions of active enzymes, taking into account fast rotational diffusion of enzymes. Similar estimates for diffusion enhancement of passive particles in lipid bilayers are derived in subsection 3.2.
The results are discussed in Section 4. There, we analyze the available experimental and computational data for diffusion enhancement in, respectively, subsections 4.1 and 4.2. Conclusions and an outline for the perspectives of further research are provided in Section 5.

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 and 2 interacting via a potential u(r) that depends on the distance r = |r 1 − r 2 | between them. The forces acting on the particles are f 1 = −∂ u/∂ r 1 = f and f 2 = − f . If the dimer is immersed into a viscous fluid, the velocity V of the hydrodynamic flow far enough from the dimer is approximately given by 12 where G αβ (R) is the mobility tensor depending on the position R of the dimer with respect to the observation point, e = (r 1 − r 2 )/r is the unit orientation vector of the dimer, and m = f r is the magnitude of the force dipole. Summation over repeated indices is assumed. The force dipole is present only if there are nonvanishing 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 much larger than their size.
The minimal active dimer model has been proposed 12,21 (see also review 23 ) to imitate mechanochemical conformational motions accompanying a catalytic turnover cycle in an enzyme. Note 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 0 and stiffness k 0 . A substrate particle (red) arrives (A) and binds as a ligand to the dimer by forming an additional elastic link with stiffness κ that connects the two beads (B). The natural length c of this additional link is taken to be shorter than 0 . Therefore, it tends to contract the dimer until a new equilibrium conformation (C) with a certain distance 1 between the beads is reached. Once this has taken place, a chemical reaction, that converts the ligand from the substrate 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 1 that is shorter than the natural length 0 . Therefore, the spring expands and the domains move apart 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 and therefore we do not consider reverse product binding events. Moreover, possible dissociation events for the substrate are neglected assuming that its affinity is high. Note that, since the product is immediately released once it has been formed, the ligand inside our model enzyme is always only in the substrate form. 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 and where x is the distance between the beads and The overdamped dynamics of the dimer in the ligand state s is described by the Langevin equation where γ is the mobility coefficient. To account for thermal fluctuations, this equation includes thermal noise, where k B 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 study 22 of diffusion enhancement for a single dimer itself. In the Oseen approximation, such interaction terms are proportional to the small parameter a/ 0 , leading to corrections of the same order for the force dipoles, neglected by us.
Stochastic transitions between the two ligand states take place at constant rates v 0 and v 1 within narrow windows of width ρ near x = 0 and x = 1 . If probability distributions p s (x,t) are introduced, they obey a system of two coupled Fokker-Planck equations and where u 0 (x) = v 0 for 0 − ρ < x < 0 + ρ and vanishes outside of this interval; u 1 (x) = v 1 for 1 − ρ < x < 1 + ρ and zero outside the interval. Note that the rate v 0 of substrate binding is proportional to the substrate concentration.
If the transition windows are very narrow, i.e., ρ 0 and ρ 1 , one can use the approximation where ν 0 = 2v 0 ρ and ν 1 = 2v 1 ρ. Fig. 2 shows the energy diagram of the model. Within each cycle, the dimer dissipates in mechanochemical motions the energy ∆E = ∆E 0 + ∆E 1 which is furthermore equal to the difference E sub − E prod of the energy E sub = E 1 ( 0 ) − E 0 ( 0 ) supplied with the substrate and the energy E prod = E 1 ( 1 ) − E 0 ( 1 ) removed with the product. We have The energy difference ∆E is always positive and, hence, the considered active dimer represents an exothermic enzyme.
The force dipole of the active dimer is m = k 0 ( 0 − x)x for s = 0 and m = k 1 ( 1 − x)x for s = 1. Note that therefore m ≤ k 0 2 0 /4 for s = 0 and m ≤ k 1 2 1 /4 for s = 1. When the transition windows are narrow, the probability rate w 0 that substrate binding, i.e., a transition to state s = 1, occurs per unit time in the state s = 0 is approximately On the other hand, the probability rate w 1 that product release, i.e., a transition to state s = 0, occurs per unit time in the state s = 1 is then approximately given by These equations are derived in Appendix A. Moreover, the characteristic relaxation times of the dimer in the states s = 0 and s = 1 are, respectively, τ 0 = (γk 0 ) −1 and τ 1 = (γk 1 ) −1 . The parameter combinations w 0 τ 0 and w 1 τ 1 play an important role in determining the kinetic regimes. If the condition w 0 τ 0 1 is satisfied, equilibration to thermal distribution in the state s = 0 usually takes place before a transition to the state s = 1, i.e., binding of a substrate, occurs. If the opposite condition w 0 τ 0 1 holds, such transition takes place immediately after the transition window at x = 0 is reached. If w 1 τ 1 1, the equilibration takes place in the state s = 1 before a transition to the state s = 0, i.e., the reaction and the product release, occurs. In the opposite limit with w 1 τ 1 1, the reaction takes place and product becomes released immediately once the respective window at x = 1 is reached.
Note that, because the rate w 0 is proportional to substrate concentration, the condition w 0 τ 0 1 corresponds to the substrate saturation regime for the considered model enzyme. The condition w 1 τ 1 1 implies that the enzyme waits a long time before the product is released.

Approximate analytical results for force dipoles
At thermal equilibrium in the absence of substrate, p 1 (x) = 0 and Since m = k 0 ( 0 −x)x, one can easily find the equilibrium statistical distribution for force dipoles by using the condition P eq (m)dm = p 0 (x)dx. Using, for convenience, the dimensionless force dipole magnitude m = m/(k 0 2 0 ) and dimensionless temperature θ = k B T /(k 0 2 0 ), we get If θ 1, this distribution is approximately Gaussian and localized at m = 0, i.e., Using the distribution in eqn (14), one finds that the mean force dipole is The correlation function C(t) = ∆m(t)∆m(0) for variations ∆m = m − m of force dipoles is 20 where τ 0 = (γk 0 ) −1 is the characteristic relaxation time for the dimer in the state s = 0. As shown in Appendix B, the exact relation m = −k B T holds for the dimer in any steady state and, therefore, in any of these limits.
For an active dimer, approximate analytical estimates can be obtained in the four characteristic limits described below. The two of 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 release occur once an appropriate conformation (x = 1 ) has been reached. In regime D, the overall kinetic rate is, on the other hand, limited by the waiting time for product formation and release.
A. The limit of w 0 τ 0 1 and w 1 τ 1 1 If these conditions are satisfied, binding of the substrate and product release have large waiting times. In this limit, there are two almost independent equilibrium subpopulations of dimers in the states s = 0 and s = 1. The relative weights of the subpopulations are w 1 /(w 1 + w 0 ) and w 0 /(w 1 + w 0 ). Therefore, all statistical properties are given by the sums of contributions from different states taken with the respective weights. Particularly, the correlation function of force dipoles is We can use the above equation to determine the non-equilibrium part of the fluctuation intensity of force dipoles Because ∆m 2 = C(0), we have As follows from eqn (17), the equilibrium fluctuation intensity is Since the effective binding rate w 0 of the substrate is proportional to its concentration c, i.e., w 0 = ηc, eqn (20) yields the Michaelis-Menten form of the dependence of ∆m 2 A on the substrate concentration.
Remarkably, the catalytic activity of the model enzyme can thus lead not only to some enhancement, but also to reduction of fluctuations of the force dipoles. According to eqn (20), reduction should be observed if k 1 k B T are satisfied, thermal fluctuations can be neglected and the dimer essentially behaves as a deterministic oscillator. Then, the solution can be obtained by integrating eqn (5) with appropriate boundary conditions. This yields x(t) = 1 +( 0 − and The respective time-dependent force dipole is Hence, it is negative for s = 1 and positive for s = 0. The force dipole varies within the interval m min < m < m max , where the minimum value m min = −k 1 0 ( 0 − 1 ) is taken at t = 0, i.e., in the state s = 1 just after substrate binding, and the maximum value m max = k 0 1 ( 0 − 1 ) is reached at t = T 1 , in the 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 have also taken the values outside of this interval.
It can be checked by direct integration that the period-averaged force dipole for the deterministic active dimer is m(t) det = 0. The correlation function for the deterministic oscillating dimer is defined as the period average The explicit analytical form of this periodic correlation function is too complicated and we do not give it (analytical results for correlation functions are also omitted below in the limits C and D). The mean-square intensity of force dipoles is m(t) 2 det = C det (0). In the limit ρ → 0, we approximately have When k 1 ∼ k 0 , this equation yields the scaling m(t) 2 det ∼ k 2 0 .
C. The limit of w 0 τ 0 1 and w 1 τ 1 1 If these conditions are satisfied, the model enzyme waits a long time for binding of a substrate (because the substrate concentration is low), but then it performs a rapid reaction cycle. An approximate solution in this regime can be obtained if, additionally, the conditions k 0 2 0 , that thermal fluctuations are weak. Moreover, we shall assume that the transition window for substrate binding is narrow, i.e., the approximation in eqn (9) holds for u 0 (x).
In this case, the dependence x(t) consists of a sum of statistically independent rare pulses, each corresponding to one reaction cycle: where given by eqn (22). The pulses appear at random time moments t j and the probability of their appearance per unit time is w 0 . Moreover, we also have where Hence, this represents a random Poisson process. Its first two statistical moments are approximately m(t) = 0 and Taking into account eqn (11), we notice that, when k 1 ∼ k 0 , the scaling m 2 (t) ∼ k This situation corresponds to substrate saturation and a long waiting time for the reaction and product release in the ligand-bound state. A derivation, similar to that given above, shows that, if k 0 2 0 k B T and k 1 2 1 k B T , we approximately have m(t) = 0 and If we take into account eqn (12), it can be noticed that, when 0 is again obtained.

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 with the 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 k B T k 0 2 0 is satisfied, so that thermal fluctuations are relatively weak.
Before proceeding to simulations, the model was nondimensionalized. The dimensionless variables weret = t/τ 0 , (5) was numerically integrated, complemented by transitions between the ligand states.
In the simulations, we had 1 = 0.55 0 , k 1 = 2k 0 , and ρ = 0.01 0 . We have kept constantṽ 1 = 2, but varied the parameterṽ 0 . A relatively low dimensionless temperature θ = 0.0018 was chosen to satisfy the condition k B T k 0 2 0 . Under such choice, m(t) 2 det / ∆m 2 eq = 19.3 and w 1 τ 1 = 0.27. Note that, because of the last condition, there was a significant random variation in the waiting times for substrate conversion and product release. Moreover, waiting times for substrate binding, characterized by the rate w 0 , could also vary. These effects 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 taken place. 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), the substrate binding rate is increased. As a result, the dimer is frequently cycling, already resembling an oscillator. Nonetheless, the random 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 the absence of the substrate, given by eqn (14). It represents a narrow Gaussian peak at m = 0. The distribution atṽ 0 = v 0 τ 0 = 0.03 (red) is almost indistinguishable from it. The blue curve is the distribution for active dimers corresponding to Fig. 3(b). Now, the distribution is more broad and the central peak is smaller. 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 v 0 , 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 by eqn (25).
Normalized correlation functions of force dipoles at different substrate binding rates are shown in Fig. 6. In the absence of the substrate (for v 0 = 0) the dependence is monotonous (it is given by eqn (17)). As the substrate concentration is increased, 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 in Fig. 6) to the dependence depend on the substrate binding rate. The oscillation period under saturation conditions is still larger than T c /τ 0 = 5.7 for the deterministic dimer according to eqn (23). This is because of an additional waiting time for product release. The characteristic relaxation time is about 1/(Γτ 0 ) = 2.
It should be stressed that the form in eqn (30) of the correlation 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 nonharmonic, as seen in Fig. 3.
There are two effects that make the dimer model stochastic, i.e., the thermal noise in the dynamical equation (5) and random transitions between the ligand states s = 0 and s = 1. When θ → 0, the thermal noise vanishes, but random transitions between the states nonetheless remain. This second stochastic effect 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 solution of the master equations (7) and (8).

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 do this, the relationship between such a simple model and the actual proteins needs to be first discussed. Proteins fold into a definite conformation that however incorporates many different substates. Slow dynamics of proteins represents wandering over a Markov network of such metastable 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 cholesterol oxidase revealed that thermal conformational fluctuations in this enzyme, in the absence of the substrate, had correlations persisting even 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 described by a set of effective collective coordinates.
In a detailed study of adenylate kinase, 29 combining all-atom MD simulations with single-molecule fluorescence resonance energy 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 the final catalytically efficient 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 of protein machines and molecular motors, such as myosin V and F 1 -ATPase 30 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 mo-tions of the domains that can be often characterized by a single coordinate, such as a hinge angle or a distance between the centers of 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 to estimate the magnitude of self-propulsion effects in the enzymes. 7 In the framework of the active dimer model, statistical properties of force dipoles in different kinetic regimes can be analyzed and characteristic order-of-magnitude estimates for the intensity of such dipoles for typical enzymes and protein machines can be derived.
In Section 2.2, four kinetic regimes have been outlined. The two of 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 high catalytic activity and, thus, the strongest nonequilibrium force dipoles can be expected.
In regime B, mechanochemical motions are limiting the overall catalytic rate. In other words, product formation and its release occur once an appropriate conformation (x = 1 ) has been reached. Such regime is characteristic, for example, for adenylate 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, catalytic hydrolysis 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 the closed 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 stiffness constants k 0 and k 1 and inter-domain distances 0 and 1 in the open (s = 0) and closed (s = 1) conformations, respectively. The typical size of a protein is of the order of tens of nanometers and this would be also the characteristic distances 0 and 1 between the domains. Moreover, if the open and closed states are distinctly different, as, for example, in adenylate kinase or myosin, the change ∆ = 0 − 1 is comparable in magnitude to 0 and 1 . As characteristic values for order-of-magnitude estimates, one can, for example, choose 0 = 10 nm and 1 = 5 nm in the open and the closed states, respectively.
In the active dimer, two domains (beads) are connected by a spring. In real proteins, they can be, instead, connected by a hinge with the elastic energy equilibrium angle Θ 0 . This can also be approximately written as so that the hinge is described as an elastic spring with x = Θ, 0 = Θ 0 and the effective stiffness k = K/ 2 , where is the characteristic linear size of the domains connected by the hinge.
The stiffness of the converter hinge in myosin V was estimated in single-molecule experiments by Kinoshita with coworkers 33 to be about K = 5 k B T /rad 2 both in the open and the closed states. On the other hand, the data of high-speed AFM observations by Ando with coworkers 34 corresponds to a higher value of K = 23 k B T /rad 2 . The difference may be due to the fact that the hinge becomes softer for larger angles. In our estimates below, we take K = 10 k B T /rad 2 . Choosing = 10 nm, this leads to k = 0.1 k B T /nm 2 .
As noticed above, in adenylate kinase, the overall turnover rate under substrate saturation in an enzyme is limited by conformational transitions between the open and closed states (and hence the turnover rate is about 10 3 s −1 ). The maximum intensity of force dipoles can be estimated by using eqn (25). If the parameter values k 0 = k 1 = 0.1 k B T /nm 2 , 0 = 10 nm, 1 = 5 nm and ρ = 1 nm are chosen, the non-equilibrium mean-square fluctuation intensity ∆m 2 A of force dipoles in such enzymes is estimated as approximately 80 [pN·nm] 2 . This is similar to the previous estimate of 100 [pN·nm] 2 in Ref. 12 based on 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 product formation and its release. In this case, the intensity of force dipoles can be estimated using eqn (29). There, the rate w 1 of product formation and release is approximately the same as the overall turnover rate, whereas τ 1 corresponds to the conformational transition time. Choosing the turnover rate of 15 s −1 , as in myosin V, and the conformational transition time of 1 ms and keeping the same other parameters as above, the non-equilibrium mean-square fluctuation intensity of force dipoles can then be estimated as 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 µs) or urease (59 µs). 36 Moreover, transition times from open to close confirmations can be also very short in some enzymes. For example, for phosphoglycerate kinase (PGK), neutron spin-echo spectroscopy yields conformational transition 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.
Equation (10) relates the energy (generally, enthalpy) dissipated in mechanochemical motions within the turnover cycle of an active dimer to the stiffness 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, k 0 = k 1 = 0.1 k B T /nm 2 and 0 − 1 = 10 nm for myosin, we obtain ∆E = 10 k B T , which is in reasonable agreement with the energy of about 20 k B T supplied to this molecular motor with ATP (only half of this energy is used in the power stroke).
Note that, assuming for simplicity that k 0 = k 1 = k, eqn (10) can be also written as thus expressing the stiffness in terms of the energy ∆E dissipated in mechanochemical motions and the conformational change An enzyme is stiffer 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 conformational transitions 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 form where ζ 0 is a dimensionless factor of order unity that also includes the logarithmic term and we have taken k 0 = k 1 = k. Note that this estimate holds assuming that the force dipoles in the catalytically active enzyme are much stronger than those due to thermal fluctuations in the absence of substrate.
Substituting k from eqn (33), a simple order-of-magnitude estimate is obtained where ζ 1 is another dimensionless factor of order unity. Moreover, by using eqn (21) and (33), we furthermore get if the condition k 2 0 k B T holds. These results show that the intensity of force dipoles is strongly sensitive to the magnitude of mechanochemical motions within the turnover cycle. Moreover, they show that, in strongly exothermic enzymes, force dipoles are greatly enhanced when catalytic activity takes place.
The above-mentioned catalase and urease enzymes are not only exceptionally fast, but also highly exothermic, with ∆H = 100 kJ/mol for catalase and ∆H = 59.6 kJ/mol for urease. 36 Hence, large energies of 42 k B T or 25 k B T are released in them and 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 proposed 36 that chemoacoustic intramolecular effects caused by strong heat release 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.

Diffusion 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 and numerical estimates for diffusion enhancement of passive particle in solutions of active enzymes. In the previous studies, 12,14 the magnitude of diffusion enhancement has been expressed in terms 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 difficulties in the interpretation of experimental results and in the analysis of the computational data. In this section, we use the statistical properties of force dipoles for active dimers, determined above in Section 2, to estimate the diffusion enhancement effects for enzymes in water solution and for active protein inclusions in biomembranes.

Diffusion effects of enzymes in water solutions
As previously shown 12,14 , the change D A in the diffusion coefficient of passive tracer particles in a three-dimensional (3D) solution is given by where n is the concentration of active enzymes, µ is viscosity, and cut is a microscopic cut-off length of the order of a protein size. Moreover, we have where C(t) is the correlation function of force dipoles corresponding to the enzymes and σ (t) is the orientational correlation function for them; χ eq is given by the same equation, but with C(t) replaced by C eq (t).
The orientational correlation function has the form where τ rot is the orientational correlation time. As seen from eqn (37), (38) and (39), the magnitude of diffusion enhancement is sensitive to the relationship between the correlation time of force dipole and the orientational correlation time.
According to the Stokes equation, rotational diffusion coefficient for a spherical particle of radius R is The orientational correlation time is defined as τ rot = 1/D rot . Since proteins are not spheres, their orientational correlation times are shorter than given by the Stokes estimate. Even in crowded solutions, they do not exceed a microsecond. [40][41][42] For active dimers, the orientational correlation times can be estimated by using eqn (40) with R = 0 .
For the active dimer model, the correlation functions of force dipoles and, therefore, their correlation times were analytically determined at equilibrium [see eqn (17)] and in the limit A [see eqn (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 τ 0 and τ 1 . In the discussion of conformational relaxation phenomena in proteins in subsection 2.4, we have noticed that slow conformational relaxation processes, involving relative domain motions in real enzymes, would usually lie in the microsecond to millisecond range. Hence, correlation times for force dipoles of the enzymes would be typically longer than their orientational correlation times.
Below, we assume that the orientational correlation time is much shorter than the correlation time for force dipoles. By using eqn (39) and putting C(t) ≈ C(0) = ∆m 2 in eqn (38), we approximately find Hence, the change in the diffusion coefficient can be estimated as Substituting approximate analytical expressions for ∆m 2 A in different limiting regimes, obtained in subsection 2.2, or using the numerical data from subsection 2.3, we can now estimate and analyze diffusion enhancement.
Particularly, it was found in subsection 2.3 that ∆m 2 A for active dimers has a Michaelis-Menten dependence on the substrate concentration. As follows from the above equation, the same dependence should hold for D A . 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 subsections 2.2 and 2.3, on the properties of the turnover cycle of an enzyme. The highest intensity is found in regime B, i.e., when there is no long waiting time for product release, and the turnover rates are limited by conformational transitions within the catalytic cycle (adenylate kinase is an example of such an enzyme).
In such asymptotic regime, ∆m 2 A is given by eqn (35) provided that the condition k 2 0 k B T holds. Substituting this expression into eqn (42), we obtain Here, R 0 is the radius of a tracer particle, 0 is the characteristic size of an enzyme (i.e., the dimer length), ∆ specifies the magnitude of the conformational change, ∆E is the free energy supplied with the substrate and dissipated within each cycle. For the equi-librium diffusion constant D T , the Stokes equation has been used. Moreover, the microscopic cut-off length is chosen as cut = 0 + R 0 and ν = 4πζ 1 /5 is a numerical factor of order unity.
Equation (43) can be employed to estimate the maximum relative diffusion 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 ∆E = 10 k B T and ∆ = 0.1 0 . As the enzyme concentration, we take n = 1 µM. This corresponds to a noncrowded solution where the mean distance between the enzymes is about ten times larger than their size ( ∼ 10 0 ). Moreover, we consider passive particles with the sizes comparable to that of an enzyme (R 0 ∼ 0 ). Under these conditions, we have D A ∼ 10D T , i.e., diffusion of tracer particles is ten times faster in the solution of catalytically active enzymes.
Dependence of diffusion enhancement on the orientational correlation time is additionally discussed in Appendix D.

Diffusion effects of active protein inclusions in biomembranes
It is known that, on the length scales shorter than the Saffman-Delbrü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 cyclically change 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 difference to water solution is that, for the biomembranes as 2D fluids, hydrodynamic diffusion enhancement effects are non-local. For such systems, eqn (37) is replaced by 12,14 Here, χ is again given by eqn (38) with σ (t) being the planar orientational correlation function for protein inclusions. Moreover, µ 2D is the 2D viscosity of the lipid bilayer, related as µ 2D = hµ 3D to its 3D viscosity µ 3D (where h is the bilayer thickness); n 2D 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 R m (shorter than the Saffman-Delbrück length) within a membrane. Then, diffusion enhancement for a passive particle of radius R 0 located in the center of the disc is 12,14 where ζ m = (1/32π) ln(R m / cut ), cut = R 0 + 0 , and χ is given by the integral in eqn (38) where, however, σ (t) is the planar orien-tational correlation function for proteins inside a membrane. The viscosity µ 3D of lipid bilayers is about 10 3 times higher than that of water and, therefore, both translational and rotational diffusion is much slower in them. From experiments, it is known that diffusion constants for proteins in lipid bilayers are about D T = 10 −10 cm 2 /s, i.e., about 10 3 times smaller than in water for similar proteins. One can therefore expect that rotational diffusion of proteins in lipid bilayers would be slowed by about a factor of 10 3 too, yielding orientational correlation times τ rot that might approach a millisecond, still being shorter than the turnover time of an enzyme.
The magnitude of diffusion enhancement in eqn (46) can be determined by modeling protein inclusions as active dimers that lie flat in the membrane. Then, the same estimate (35) for ∆m 2 A can be used. Combining all terms, diffusion enhancement in eqn (46) for a passive particle in the center of a protein raft approximately is where the dimensionless prefactor is ν m = ζ 1 ζ m and 2D = n is the mean distance between inclusions in the membrane.
To obtain a characteristic order-of-magnitude estimate, the 3D viscosity of the lipid bilayer is chosen as µ 3D = 1 Pa·s and the thickness of the bilayer as h = 1 nm. For protein inclusions, we assume that ∆E = 10 k B T and ∆ ∼ 0 . The orientational correlation time is taken to be τ rot = 100 µs and the mean lateral distance between the proteins is 2D = 10 nm. For such parameter values, the maximal possible diffusion enhancement under substrate saturation conditions is about D A = 10 −9 cm 2 /s. For comparison, Brownian diffusion constants for proteins in lipid bilayers are of the order of 10 −10 cm 2 /s and diffusion constants for lipids are about 10 −8 cm 2 /s.

Discussion
Using the results of our study, available experimental and computational data on diffusion enhancement in solutions of catalytically active enzymes can be discussed.

Experimental data
Diffusion enhancement for the enzymes has been reported in solutions of several catalytically active enzymes, at the concentrations varying between 1 nM and 10 nM. 36,[43][44][45][46] With the exception of aldolase 43 (for which, however, the enhancement could not be independently confirmed 47 ), all these enzymes were exothermic and had high turnover rates of about 10 4 s −1 . 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 can be understood in the framework of the original theory 12 and its subsequent extensions, including the present work. The fact that a significant diffusion enhancement (by tens of percent) was observed already at low nanomole concentrations can still be per-haps 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 the experimentally observed diffusion enhancement on the enzyme concentration clearly contradicts the theory 12 where diffusion enhancement arises as a collective hydrodynamic effect. Effectively, diffusion enhancement was observed in the experiments 36,43-46 already for single molecules of enzymes.
When our study was completed, an interesting publication 48 has appeared where diffusion enhancement was demonstrated by a different method for several other reactions. Since catalyst's diffusion was not affected by its concentration, this was again a single-particle property not covered by the theory. 12 Experiments on optical tracking of particles in animal cells 49 and in bacteria or yeast 50 have been furthermore performed. They have shown that, when metabolism was suppressed (by depletion of ATP), diffusion dropped to undetectable levels 49,51 or it was much slowed down and replaced by subdiffusion characteristic for a colloidal glass. 50 Strong reduction of diffusion under metabolism suppression was moreover found in various cytoplasm extracts. 52 It should be also noted that diffusion enhancement has been experimentally observed within chromatin in a living biological cell. 53 This was explained by active operation of molecular machines involved in transcription and translation of DNA. 54 The cytoplasm of a living cell represents a crowded solution of proteins. In bacteria, the volume fraction of proteins in cytosol is about 30 percent, 55 with the highest concentrations of the order of 100 µM reached for glycolysis enzymes. Most of the enzymes in the cell are mechanochemical, i.e., they exhibit conformational changes in their catalytic cycles. Typical turnover times of enzymes in a biological cell are of the order of 10 ms.
According to the previous 12 and current estimates, substantial diffusion enhancement due to hydrodynamic collective effects should thus be expected under metabolism in the cytoplasm. There are, however, also other mechanisms that can contribute to diffusion enhancement in the cells.
The cytoskeleton of animal cells represents an active gel, with numerous myosin molecular motors operating within it. It is known that the activity of the motors can lead to development of non-equilibrium fluctuations in the cytoskeleton which induce in turn fluctuations and diffusion enhancement in the cytosol. 17,51,56 The skeleton of bacteria and yeast is however passive; moreover, metabolic diffusion enhancement in such cells could also be observed when their skeleton was chemically resolved. 50 Therefore, the active gel mechanism 56 cannot account for the effects 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 and relaxation phenomena strongly slowed down in them. 57 Indeed, such behavior could be observed both in the cells 50 and in the extracts 52 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 k B T per a protein molecule per a cycle, can lead to diffusion enhancement by one order of magnitude. 58 This provides an additional, non-hydrodynamic, mechanism that can contribute to the experimentally observed diffusion enhancement in living biological cells.
In summary, the analysis of the available experimental data reveals that the predicted diffusion enhancement 12 for passive particles caused by collective catalytic activity of enzymes could not so far been reliably confirmed.

Computational data
Large-scale computer simulations for colloids of active dimers have been performed by Dennison,Kapral and Stark. 20 In these simulations, the solvent was explicitly included and the multiparticle collision dynamics (MPCD) approximation 60 was employed, thus allowing to fully account for hydrodynamic effects.
To facilitate the comparison, we first give a summary of the essential parameter values in the study, 20 using the current notations employed by us. The natural lengths of the dimer in two ligand states were 0 and 1 = 0 /2, and the spring constants were k 0 and k 1 = 2k 0 . The dimensionless spring constant k 2 0 /(k B T ), characterizing stiffness of the dimer, was varied between 144 and 1440. The energy ∆E = (1/2)(k 0 + k 1 )( 0 − 1 ) 2 , supplied to a dimer and dissipated by it as heat within a single cycle, was changing therefore between 121.5 k B T and 1215 k B T . The simulations were performed under substrate saturation conditions. Product formation and release were possible within a window of half-width ρ = 0.025 0 near x = 1 . The rate v 1 of this transition could be varied in the simulations by a factor of 5.
The Langevin equation (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 C eq (t) was computed yielding the correlation time for fluctuations of its force dipole; this function could be well fitted to the theoretical dependence in eqn (17). Note that, when k 0 1, the relaxation time τ 0 = (γk 0 ) −1 of the dimer should be close to this correlation time. Moreover, we have τ 1 = (γk 1 ) −1 = τ 0 /2. Using such estimates, it can be shown that w 1 τ 1 varied between 0.001 and 0.1 in the simulations. 20 Because substrate saturation was assumed, conditions w 0 τ 0 1 and w 1 τ 1 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 v 0 τ 0 = 3 in Fig. 6. The correlation times varied, but remained of the same order of magnitude as the correlation time of the passive dimer. The forcedipole intensity ∆m 2 of active dimers was by about an order of magnitude larger than ∆m 2 eq for the passive ones. Depending on the parameters, it scaled as k α 0 with the exponent α in the range between 1.2 and 1.6, comparable with the exponent of 1.5 in eqn (29).
Orientational correlation functions σ (t) were furthermore computed for single dimers. 20 Remarkably, it was found that the orientational correlation time τ rot was sensitive to the conformational activity of the dimer, getting shorter by about an order of magnitude when such activity was switched on. Nonetheless, in all simulations τ rot was larger than the force dipole correlation time.
Multiparticle 3D computer simulations of colloids formed by active dimers were further performed. 20 In the simulations, the truncated potential for r < 2 1/24 (2r 0 ) and zero otherwise, with ε = 2.5 k B T and r 0 = 1.075 0 , was used to describe steep repulsive interactions between the beads belonging to different dimers. The interaction radius r 0 was chosen as defining the radius of a bead.
Since distances 0 and 1 = 0.5 0 in the open and closed dimer conformations were both smaller than 2r 0 = 2.15 0 , large overlaps between the beads in a dimer were present in the simulations. However, this did not affect 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.5 0 .
The volume fraction φ occupied by dimers was determined by taking into account the overlaps, but assuming that all dimers were in the equilibrium open state with the length of 0 . 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 effects, diffusion of a passive particle in the system of inactive dimers decreased with the volume fraction of them. The diffusion reduction at the highest taken volume fraction φ = 0.266 was less than ten percent, indicating that this colloidal system was still far from the glass transition threshold. 57 When the dimers were active, diffusion of tracers was increasing instead with the dimer volume fraction φ . For the most stiff active dimers with k Thus, collective hydrodynamic effects of active enzymes on diffusion of passive particles could be computationally confirmed. To speed up the calculations, model enzymes in the study 20 were chosen however to be unusually rapid (with the turnover times shorter than the rotational diffusion time) and unusually exothermic (with the heat release of hundreds of k B T per a turnover cycle). It would be therefore important to undertake such simulations also for the parameters closer approaching those of the real enzymes.

Conclusions and outlook
To our knowledge, the present work is the first study where hydrodynamic force dipoles of mechanochemical enzymes have been 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, such as adenylate kinase, and for protein machines, such as myosin, have been obtained.
We have also examined for what kinds of enzymes strong hydrodynamic effects may be expected. Our analysis reveals that, in the framework of the investigated model, these should be very rapid, highly exothermic and stiff enzymes, where the energy is dissipated in mechanical motions of a small amplitude. It is interesting to note that these general conditions are indeed satisfied, for example, for catalase or urease.
Using the derived statistical properties of force dipoles in the dimer model, more accurate estimates for diffusion enhancement for surrounding passive particles in solutions of active enzymes were obtained.
Based on these results, currently available experimental and computational data has been examined. We have concluded that, while the collective hydrodynamic effects of diffusion 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 data on diffusion enhancement in weak nanomole solutions of several fast exothermic enzymes cannot be explained in the framework of the theory 12 and alternative explanations for them should be sought. In experimental studies of diffusion phenomena in living cells and in cellular extracts, additional work is needed to distinguish possible hydrodynamic contributions from the effects of direct collisions between active proteins and the resulting kinetic crowding effects. Large-scale numerical simulations of crowded active colloids including hydrodynamic interactions between the particles are to be performed. It should be also pointed out that, although the effects of diffusion 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 diffusion enhancement effects. Especially, we have not considered possible origins of diffusion enhancement for single catalytically active enzymes, even though this question attracts much attention in view of the recent research. 45,46,48 Our focus was on diffusion enhancement for passive particles caused by hydrodynamic collective stirring of the solution by a population of active particles that cyclically 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 effects, accompanying functional conformational transitions, in allatom or coarse-grained molecular dynamics simulations for spe-cific enzymes and protein machines.

Conflicts of interest
There are no conflicts to declare.

A Transition probabilities
When transitions between the states s = 0 and s = 1 are rare, the solution of the master equations (7) and (8) can be approximately sought in the form where π s (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.
Substituting these expressions into eqn (7) and (8) and integrating over the variable x, one finds that the probabilities π s obey classical master equations for a two-level system, and Here w 0 and w 1 are effective rates of transitions between the states given by and The involved probability distributions in the statistically stationary case are and If the transition windows are narrow, approximations in eqn (9) can furthermore be used, so that we obtain w 0 = ν 0 p(x = 0 |s = 0), w 1 = ν 1 p(x = 1 |s = 1).
Thus, using the above expressions for distance distributions, we finally get and In the steady state, the probabilities are

B Average force dipole
Let us consider the second statistical moment x 2 . In a steady state, its time derivative is zero. On the other hand, by using eqn (7)-(8) and integrating by parts, we find Thus, we straightforwardly obtain that, for an active dimer in any statistically steady state, m = −k B T .
Note that here and also in the equations below, the integration limits over x are taken as +∞ and −∞. The actual limits are automatically selected by probability distributions p 0 (x) and p 1 (x).
The general solution of eqn (62) is s (x)e −λ n t + c.c.
where λ n and q (n) are eigenvalues and eigenvectors of the linear operatorL,L q (n) = λ n q (n) , and decomposition coefficients A n are determined by initial conditions.
Because the master equation must have a stable stationary solution, the operatorL should always possess a zero eigenvalue λ 0 = 0 and, furthermore, condition Re λ n > 0 should hold for all other eigenvalues n. 61 Generally, the eigenvectors can be ordered according to the increase of Re λ n (and therefore we can enumerate the eigenvalues in such a way that 0 < Re λ 1 ≤ Re λ 2 ≤ Re λ 3 ≤ ...). The stationary probability distributionp(x) coincides with the eigenvector q (0) (x).
The conditional probability G(x, s,t|x 0 , s 0 ) gives the probability to find the dimer in various states (x, s) at time t provided that it was in the state (x 0 , s 0 ) at time t = 0. It represents a special solution of the master equation (62) given by G(x, s,t|x 0 , s 0 ) = ∞ ∑ n=0 a n (x 0 , s 0 )q (n) s (x)e −λ n t + c.c.
where a n (x 0 , s 0 ) are the coefficients of decomposition of this initial condition over eigenvectors q (n) .
The force dipole m depends on the distance x between the domains and on the dimer state s, i.e., m(t) = m(x(t), s(t)). Therefore, in the statistically stationary state we have By using eqn (69) and (70), we find that, in the statistically sta-tionary state, the correlation function of force dipoles is If we retain in this decomposition only the first, most slowly decaying term, this yields C(t) ≈ B 1 e −λ 1 |t| + c.c. = C(0) cos α e −Γ|t| cos(Ω|t| − α).
Our numerical simulations, described in Section 2.3, have shown that, in the regimes approaching a deterministic oscillatory dimer, the correlation functions of force dipoles could be well fitted to the above dependence. This suggests that contributions from the higher, more rapidly decaying relaxation modes n > 1 have been indeed relatively small. As generally known, 62 noisy oscillators possess a slowly relaxing mode that corresponds to diffusion of the oscillation phase. It can be expected that, under chosen conditions, such a mode has been dominating the correlation functions for oscillatory dimers.
This yields a non-monotonous dependence of χ on the orientational correlation time. If the phase shift α is small and can be neglected (cf. Fig. 7), the maximum value χ max is reached at τ rot = (Ω − Γ) −1 and we have where is the limit of χ when τ rot Γ −1 and τ rot Ω −1 . These results allow us to discuss how the diffusion enhancement would depend on the orientational correlational time τ rot , not assuming that it is much shorter than the correlation time for force dipoles. If the approximation in eqn (30) holds, diffusion enhancement is determined by eqn (37) where χ is given by eqn (75). The diffusion enhancement depends nonmonotonously on the orientational correlation time. It increases linearly with τ rot at short times, then reaches a maximum at τ rot = (Ω − Γ) −1 and finally saturates at large orientational correlation times. For example, if we take the values Γ ≈ 1/(2τ 0 ) and Ω ≈ π/(3τ 0 ) corresponding to substrate saturation in Fig. 7, the maximum diffusion enhancement would be reached at τ rot = 1.8τ 0 and, at the maximum, it will be larger by about 30 percent than in the limit τ rot τ 0 .