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

Abstract

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.

---

1 Introduction

Ligand-induced mechanochemical motions are typical for enzymes. Binding or dissociation of a ligand (i.e., substrate or product) to such proteins, as well as chemical reactions within 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¹ and adenylate kinase.² 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–7

Lipid bilayers, forming biological membranes, behave as two-dimensional (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.¹⁰ As a result, active protein inclusions might even propel themselves through biomembranes.¹¹

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,¹² increased mixing in such systems and diffusion enhancement should therefore arise. Additionally, chemotaxis-like effects in the presence of spatial gradients in the concentration or the activity of enzymes can take place.¹² 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–7

Following the original publication,¹² extensive further research has been performed.^13–20 The effects of rotational diffusion and of possible nematic ordering for enzymes were considered,¹⁴ 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.¹⁹ Multiparticle numerical simulations of active oscillatory colloids, explicitly including hydrodynamic effects, were furthermore undertaken and principal theoretical predictions could thus be verified.²⁰

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 effects in solution of such enzymes are predicted by the mean-field theory.¹² 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–22 After presenting the model, we undertake an approximate analytical investigation of statistical properties of the force dipoles corresponding to active dimers in Section 2.2, followed by a numerical study in Section 2.3. Quantitative estimates for the intensity of hydrodynamical force dipoles in real enzymes are obtained in Section 2.4.

Section 3 corresponds to the second part. Based on the active dimer results, we obtain in Section 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 Section 3.2.

The results are discussed in Section 4. There, we analyze the available experimental and computational data for diffusion enhancement in, respectively, Sections 4.1 and 4.2. Conclusions and an outline for the perspectives of further research are provided in Section 5.

2 Statistical properties of force dipoles

2.1 The active dimer model

The simplest mechanical system that gives rise to a hydrodynamical force dipole is a dimer. It consists of two beads 1 and 2 interacting via a potential u(r) that depends on the distance r = |r₁ − r₂| between them. The forces acting on the particles are f₁ = −∂u/∂r₁ = f and f₂ = −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¹²where G_αβ(R) is the mobility tensor depending on the position R of the dimer with respect to the observation point, e = (r₁ − r₂)/r is the unit orientation vector of the dimer, and m = fr is the magnitude of the force dipole. Summation over repeated indices is assumed. The force dipole is present only if there are non-vanishing net interaction forces, i.e., if the distance between the particles in a dimer continues to change. As in the study,¹² 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 ref. 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.²²

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 [small script l]₀ and stiffness k₀. 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 [small script l]_c of this additional link is taken to be shorter than [small script l]₀. Therefore, it tends to contract the dimer until a new equilibrium conformation (C) with a certain distance [small script l]₁ 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 [small script l]₁ that is shorter than the natural length [small script l]₀. 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.

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

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

andwhere 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,

〈ξ(t₁)ξ(t₂)〉 = 2γk_BTδ(t₁ − t₂), (6)

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²² of diffusion enhancement for a single dimer itself. In the Oseen approximation, such interaction terms are proportional to the small parameter a/[small script l]₀, 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₀ and v₁ within narrow windows of width ρ near x = [small script l]₀ and x = [small script l]₁. If probability distributions p_s(x,t) are introduced, they obey a system of two coupled Fokker–Planck equations

andwhere u₀(x) = v₀ for [small script l]₀ − ρ < x < [small script l]₀ + ρ and vanishes outside of this interval; u₁(x) = v₁ for [small script l]₁ − ρ < x < [small script l]₁ + ρ and zero outside the interval. Note that the rate v₀ of substrate binding is proportional to the substrate concentration.

If the transition windows are very narrow, i.e., ρ ≪ [small script l]₀ and ρ ≪ [small script l]₁, one can use the approximation

u₀(x) = ν₀δ(x − [small script l]₀), u₁(x) = ν₁δ(x − [small script l]₁), (9)

where ν₀ = 2v₀ρ and ν₁ = 2v₁ρ.

Fig. 2 shows the energy diagram of the model. Within each cycle, the dimer dissipates in mechanochemical motions the energy ΔE = ΔE₀ + ΔE₁ which is furthermore equal to the difference E_sub − E_prod of the energy E_sub = E₁([small script l]₀) − E₀([small script l]₀) supplied with the substrate and the energy E_prod = E₁([small script l]₁) − E₀([small script l]₁) removed with the product. We have

The energy difference ΔE is always positive and, hence, the considered active dimer represents an exothermic enzyme.

Fig. 2 The energy diagram of the active dimer.

The force dipole of the active dimer is m = k₀([small script l]₀ − x)x for s = 0 and m = k₁([small script l]₁ − x)x for s = 1. Note that therefore m ≤ k₀[small script l]₀²/4 for s = 0 and m ≤ k₁[small script l]₁²/4 for s = 1.

When the transition windows are narrow, the probability rate w₀ 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₁ that product release, i.e., a transition to state s = 0, occurs per unit time in the state s = 1 is then approximately given byThese equations are derived in Appendix A. Moreover, the characteristic relaxation times of the dimer in the states s = 0 and s = 1 are, respectively, τ₀ = (γk₀)⁻¹ and τ₁ = (γk₁)⁻¹.

The parameter combinations w₀τ₀ and w₁τ₁ play an important role in determining the kinetic regimes. If the condition w₀τ₀ ≪ 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₀τ₀ ≫ 1 holds, such transition takes place immediately after the transition window at x = [small script l]₀ is reached. If w₁τ₁ ≪ 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, the reaction takes place and product becomes released immediately once the respective window at x = [small script l]₁ is reached.

Note that, because the rate w₀ is proportional to substrate concentration, the condition w₀τ₀ ≫ 1 corresponds to the substrate saturation regime for the considered model enzyme. The condition w₁τ₁ ≪ 1 implies that the enzyme waits a long time before the product is released.

2.2 Approximate analytical results for force dipoles

At thermal equilibrium in the absence of substrate, p₁(x) = 0 andSince m = k₀([small script l]₀ − x)x, one can easily find the equilibrium statistical distribution for force dipoles by using the condition P_eq(m)dm = p₀(x)dx. Using, for convenience, the dimensionless force dipole magnitude m̃ = m/(k₀[small script l]₀²) and dimensionless temperature θ = k_BT/(k₀[small script l]₀²), we getIf θ ≪ 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

〈m〉_eq = −k_BT. (16)

The correlation function C(t) = 〈Δm(t)Δm(0)〉 for variations Δm = m − 〈m〉 of force dipoles is²⁰

C_eq(t) = k₀[small script l]₀²k_BTe^−|t|/τ₀ + 2(k_BT)²e^{−2|t|/τ₀}, (17)

where τ₀ = (γk₀)⁻¹ is the characteristic relaxation time for the dimer in the state s = 0. As shown in Appendix B, the exact relation 〈m〉 = −k_BT holds for the dimer in any steady state.

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 = [small script l]₁) 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₀τ₀ ≪ 1 and w₁τ₁ ≪ 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₁/(w₁ + w₀) and w₀/(w₁ + w₀). 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

〈Δm²〉_A = 〈Δm²〉 − 〈Δm²〉_eq. (19)

Because 〈Δm²〉 = C(0), we haveAs follows from eqn (17), the equilibrium fluctuation intensity is

〈Δm²〉_eq = k₀[small script l]₀²k_BT + 2(k_BT)². (21)

Since the effective binding rate w₀ of the substrate is proportional to its concentration c, i.e., w₀ = ηc, eqn (20) yields the Michaelis–Menten form of the dependence of 〈Δm²〉_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₁[small script l]₁² < k₀[small script l]₀². Under this condition, the ligand-bound dimer (s = 1) is characterized by a lower fluctuation intensity of force dipoles than the free dimer (s = 0).

B. The limit of w₀τ₀ ≫ 1 and w₁τ₁ ≫ 1. In this limit, transitions take place once the respective transitions windows are entered. If additionally the conditions k₀[small script l]₀² ≫ k_BT and k₁[small script l]₁² ≫ k_BT 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) = [small script l]₁ + ([small script l]₀ − [small script l]₁ − ρ)e^−t/τ₁ for 0 < t < T₁ and x(t) = [small script l]₀ + ([small script l]₁ − [small script l]₀ + ρ)e^{−(t−T₁)/τ₀} for T₁ < t < T_c. Here, T_c is the oscillation period of the active dimer and T₁ is the duration of the cycle time when the dimer is in the ligand-bound state s = 1. If transition windows are narrow, i.e., the condition ρ ≪ ([small script l]₀ − [small script l]₁) is satisfied, we approximately haveandThe respective time-dependent force dipole is m(t) = k₁([small script l]₁ − x)x for 0 < t < T₁ and m(t) = k₀([small script l]₀ − x)x for T₁ < t < T_c. 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₁[small script l]₀([small script l]₀ − [small script l]₁) is taken at t = 0, i.e., in the state s = 1 just after substrate binding, and the maximum value m_max = k₀[small script l]₁([small script l]₀ − [small script l]₁) is reached at t = T₁, 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)²〉_det = C_det(0). In the limit ρ → 0, we approximately have

When k₁ ∼ k₀, this equation yields the scaling 〈m(t)²〉_det ∼ k₀².

C. The limit of w₀τ₀ ≪ 1 and w₁τ₁ ≫ 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₀[small script l]₀² ≫ k_BT and k₁[small script l]₁² ≫ k_BT are satisfied, i.e., 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₀(x).

In this case, the dependence x(t) consists of a sum of statistically independent rare pulses, each corresponding to one reaction cycle:

where z(t) = [small script l]₁ + ([small script l]₀ − [small script l]₁)e^−t/τ₁ for 0 < t < T₁ and z(t) = [small script l]₀ + ([small script l]₁ − [small script l]₀)e^{−(t−T₁)/τ₀} for t > T₁, with T₁ given by eqn (22). The pulses appear at random time moments t_j and the probability of their appearance per unit time is w₀.

Moreover, we also have

where ζ(t) = k₁([small script l]₁ − z(t))z(t) for 0 < t < T₁ and ζ(t) = k₀([small script l]₀ − z(t))z(t) for t > T₁.

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₁ ∼ k₀, the scaling 〈m²(t)〉 ∼ k₀^3/2 should hold.

D. The limit of w₀τ₀ ≫ 1 and w₁τ₁ ≪ 1. 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₀[small script l]₀² ≫ k_BT and k₁[small script l]₁² ≫ k_BT, we approximately have 〈m(t)〉 = 0 andIf we take into account eqn (12), it can be noticed that, when k₁ ∼ k₀, scaling 〈m²(t)〉 ∼ k^3/2₀ is again obtained.

2.3 Numerical simulations

Numerical simulations can yield statistical properties of force dipoles for selected parameter values in the regions where there are no approximate analytical results. Below in this section, we focus on the situation under substrate saturation, but 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_BT ≪ k₀[small script l]₀² is satisfied, so that thermal fluctuations are relatively weak.

Before proceeding to simulations, the model was non-dimensionalized. The dimensionless variables were t̃ = t/τ₀, x̃ = x/[small script l]₀, and m̃ = m/(k₀[small script l]₀²). The dimensionless transition rates were ṽ₀ = v₀τ₀ and ṽ₁ = v₁τ₀, while the dimensionless temperature was θ = k_BT/(k₀[small script l]₀²). Stochastic differential eqn (5) was numerically integrated, complemented by transitions between the ligand states.

In the simulations, we had [small script l]₁ = 0.55[small script l]₀, k₁ = 2k₀, and ρ = 0.01[small script l]₀. We have kept constant ṽ₁ = 2, but varied the parameter ṽ₀. A relatively low dimensionless temperature θ = 0.0018 was chosen to satisfy the condition k_BT ≪ k₀[small script l]₀². Under such choice, 〈m(t)²〉_det/〈Δm²〉_eq = 19.3 and w₁τ₁ = 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₀, 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.

Fig. 3 Time dependence of dimensionless force dipoles m̃ = m/(k₀[small script l]₀²) on time for (a) ṽ₀ = 0.03 and (b) ṽ₀ = 3. Dashed lines show the lower bound m̃_min = −0.9 for the deterministic oscillatory dimer and the absolute upper bound m̃_max = 0.25 for force dipoles.

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 ṽ₀ = v₀τ₀ = 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.

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

The dependence of the non-equilibrium part of the fluctuation intensity of force dipoles, eqn (19), on the substrate binding rate v₀, 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).

Fig. 5 Dependence of the non-equilibrium part 〈Δm²〉_A of the fluctuation intensity of force dipoles on the substrate binding rate v₀ (dots). The solid curve is a fit to the Michaelis–Menten function.

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) 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.

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

The correlation functions could be fitted (dashed curves in Fig. 6) to the dependence

Fig. 7 shows how the dimensionless relaxation time 1/(Γτ₀), the dimensionless oscillation period 2π/(Ωτ₀) and the phase shift α depend on the substrate binding rate. The oscillation period under saturation conditions is still larger than T_c/τ₀ = 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/(Γτ₀) = 2.

Fig. 7 The dependences of the relaxation time 1/Γ (circles), oscillation period 2π/Ω (triangles) and phase shift α (squares) on substrate binding rate v₀.

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 non-harmonic, as seen in Fig. 3.

There are two effects that make the dimer model stochastic, i.e., the thermal noise in the dynamical eqn (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 eqn (7) and (8).

2.4 Estimates for hydrodynamic force dipoles of enzymes

Above, statistical properties of force dipoles were analyzed in the framework of an idealized model of the active dimer. Now, the obtained results can be applied to approximately estimate the force dipoles for real enzymes and protein machines. To 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.²⁴ 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.²⁵ 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,²⁸ thus yielding continuous slow conformational dynamics described by a set of effective collective coordinates.

In a detailed study of adenylate kinase,²⁹ 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₁-ATPase³⁰ and HCV helicase.³¹

Typically, mechanochemical enzymes and molecular machines represent proteins with domain structure. Slow functional conformational dynamics in these proteins consists in relative motions of the domains that can be often characterized by a single coordinate, such as a hinge angle or a distance between the centers of mass of two protein domains. This leads to the reduced models for proteins, with just one or a few mechanical coordinates.²³ 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.⁷ 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 non-equilibrium 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 = [small script l]₁) 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.²⁹

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.³² 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₀ and k₁ and inter-domain distances [small script l]₀ and [small script l]₁ 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 [small script l]₀ and [small script l]₁ between the domains. Moreover, if the open and closed states are distinctly different, as, for example, in adenylate kinase or myosin, the change Δ[small script l] = [small script l]₀ − [small script l]₁ is comparable in magnitude to [small script l]₀ and [small script l]₁. As characteristic values for order-of-magnitude estimates, one can, for example, choose [small script l]₀ = 10 nm and [small script l]₁ = 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

which depends on the deviation of the hinge angle Θ from the equilibrium angle Θ₀. This can also be approximately written asso that the hinge is described as an elastic spring with x = LΘ, [small script l]₀ = LΘ₀ and the effective stiffness k = K/L², where L 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³³ to be about K = 5 k_BT rad⁻² both in the open and the closed states. On the other hand, the data of high-speed AFM observations by Ando with coworkers³⁴ corresponds to a higher value of K = 23 k_BT rad⁻². 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_BT rad⁻². Choosing L = 10 nm, this leads to k = 0.1 k_BT nm⁻².

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³ s⁻¹). The maximum intensity of force dipoles can be estimated by using eqn (25). If the parameter values k₀ = k₁ = 0.1 k_BT nm⁻², [small script l]₀ = 10 nm, [small script l]₁ = 5 nm and ρ = 1 nm are chosen, the non-equilibrium mean-square fluctuation intensity 〈Δm²〉_A of force dipoles in such enzymes is estimated as approximately 80 [pN nm]². This is similar to the previous estimate of 100 [pN nm]² 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₁ of product formation and release is approximately the same as the overall turnover rate, whereas τ₁ corresponds to the conformational transition time. Choosing the turnover rate of 15 s⁻¹, 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]². 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,³⁵ and enzymes that are very fast, such as catalase (17 μs) or urease (59 μs).³⁶ 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.³⁷ This was also observed in coarse-grained MD simulations for PGK.³⁸ Therefore, it is interesting to discuss under what general conditions stronger force dipoles can be expected in enzymes.

Eqn (10) relates the energy (generally, enthalpy) dissipated 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₀ = k₁ = 0.1 k_BT nm⁻² and [small script l]₀ − [small script l]₁ = 10 nm for myosin, we obtain ΔE = 10 k_BT, which is in reasonable agreement with the energy of about 20 k_BT 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₀ = k₁ = 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 Δ[small script l] = [small script l]₀ − [small script l]₁. 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

〈Δm²〉 = ζ₀k²[small script l]₀²([small script l]₀ − [small script l]₁)², (34)

where ζ₀ is a dimensionless factor of order unity that also includes the logarithmic term and we have taken k₀ = k₁ = 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 ζ₁ is another dimensionless factor of order unity. Moreover, by using eqn (21) and (33), we furthermore getif the condition k[small script l]₀² ≫ k_BT 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.³⁶ Hence, large energies of 42 k_BT or 25 k_BT 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.³⁹ It has been previously proposed³⁶ 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.⁷ These enzymes do not have a domain structure and therefore the results of our analysis based on the dimer model are not directly applicable to them. Nonetheless, they suggest that hydrodynamic force dipoles in them may be very strong.

3 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.

3.1 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 bywhere n is the concentration of active enzymes, μ is viscosity, and [small script l]_cut is a microscopic cut-off length of the order of a protein size. Moreover, we havewhere 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

σ(t) = exp(−t/τ_rot), (39)

where τ_rot is the orientational correlation time. As seen from eqn (37)–(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–42 For active dimers, the orientational correlation times can be estimated by using eqn (40) with R = [small script l]₀.

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 τ₀ and τ₁. In the discussion of conformational relaxation phenomena in proteins in Section 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²〉 in eqn (38), we approximately find

χ − χ_eq = τ_rot〈Δm²〉_A. (41)

Hence, the change in the diffusion coefficient can be estimated asSubstituting approximate analytical expressions for 〈Δm²〉_A in different limiting regimes, obtained in Section 2.2, or using the numerical data from Section 2.3, we can now estimate and analyze diffusion enhancement.

Particularly, it was found in Section 2.3 that 〈Δm²〉_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 Sections 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²〉_A is given by eqn (35) provided that the condition k[small script l]₀² ≫ k_BT holds. Substituting this expression into eqn (42), we obtain

Here, R₀ is the radius of a tracer particle, [small script l]₀ is the characteristic size of an enzyme (i.e., the dimer length), Δ[small script l] specifies the magnitude of the conformational change, ΔE is the free energy supplied with the substrate and dissipated within each cycle, and Δ = n^−1/3 is the mean distance between active enzymes. For the equilibrium diffusion constant D_T, the Stokes equationhas been used. Moreover, the microscopic cut-off length is chosen as [small script l]_cut = [small script l]₀ + R₀ and ν = 4πζ₁/5 is a numerical factor of order unity.

Eqn (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_BT and Δ[small script l] = 0.1[small script l]₀. As the enzyme concentration, we take n = 1 μM. This corresponds to a non-crowded solution where the mean distance between the enzymes is about ten times larger than their size ([small script l] ∼ 10[small script l]₀). Moreover, we consider passive particles with the sizes comparable to that of an enzyme (R₀ ∼ [small script l]₀). 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.

3.2 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.⁸ 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.¹²

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₀ located in the center of the disc is^12,14

where ζ_m = (1/32π)ln(R_m/[small script l]_cut), [small script l]_cut = R₀ + [small script l]₀, and χ is given by the integral in eqn (38) where, however, σ(t) is the planar orientational correlation function for proteins inside a membrane.

The viscosity μ_3D of lipid bilayers is about 10³ 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⁻¹⁰ cm² s⁻¹, i.e., about 10³ 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³ 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²〉_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 = ζ₁ζ_m and [small script l]_2D = n_2D^−1/2 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_BT and Δ[small script l] ∼ [small script l]₀. The orientational correlation time is taken to be τ_rot = 100 μs and the mean lateral distance between the proteins is [small script l]_2D = 10 nm. For such parameter values, the maximal possible diffusion enhancement under substrate saturation conditions is about D_A = 10⁻⁹ cm² s⁻¹. For comparison, Brownian diffusion constants for proteins in lipid bilayers are of the order of 10⁻¹⁰ cm² s⁻¹ and diffusion constants for lipids are about 10⁻⁸ cm² s⁻¹.

4 Discussion

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

4.1 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–46 With the exception of aldolase⁴³ (for which, however, the enhancement could not be independently confirmed⁴⁷), all these enzymes were exothermic and had high turnover rates of about 10⁴ s⁻¹. 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.⁴⁶

It does not seem plausible that such experimental data can be understood in the framework of the original theory¹² 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 perhaps explained by assuming that, for some reasons, the force dipoles of specific enzymes with high catalytic turnover rates were exceptionally strong. However, the absence of a dependence of the experimentally observed diffusion enhancement on the enzyme concentration clearly contradicts the theory¹² 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⁴⁸ 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.¹²

Experiments on optical tracking of particles in animal cells⁴⁹ and in bacteria or yeast⁵⁰ 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.⁵⁰ Strong reduction of diffusion under metabolism suppression was moreover found in various cytoplasm extracts.⁵²

It should be also noted that diffusion enhancement has been experimentally observed within chromatin in a living biological cell.⁵³ This was explained by active operation of molecular machines involved in transcription and translation of DNA.⁵⁴

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,⁵⁵ 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¹² 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.⁵⁰ Therefore, the active gel mechanism⁵⁶ 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.⁵⁷ Indeed, such behavior could be observed both in the cells⁵⁰ and in the extracts⁵² 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_BT per a protein molecule per a cycle, can lead to diffusion enhancement by one order of magnitude.⁵⁸ 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¹² for passive particles caused by collective catalytic activity of enzymes could not have been reliably confirmed.

4.2 Computational data

Large-scale computer simulations for colloids of active dimers have been performed by Dennison, Kapral and Stark.²⁰ In these simulations, the solvent was explicitly included and the multiparticle collision dynamics (MPCD) approximation⁶⁰ 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,²⁰ using the current notations employed by us. The natural lengths of the dimer in two ligand states were [small script l]₀ and [small script l]₁ = [small script l]₀/2, and the spring constants were k₀ and k₁ = 2k₀. The dimensionless spring constant k[small script l]₀²/(k_BT), characterizing stiffness of the dimer, was varied between 144 and 1440. The energy ΔE = (1/2)(k₀ + k₁)([small script l]₀ − [small script l]₁)², supplied to a dimer and dissipated by it as heat within a single cycle, was changing therefore between 121.5 k_BT and 1215 k_BT. The simulations were performed under substrate saturation conditions. Product formation and release were possible within a window of half-width ρ = 0.025[small script l]₀ near x = [small script l]₁. The rate v₁ of this transition could be varied in the simulations by a factor of 5.

The Langevin eqn (5) with viscous friction and thermal noise was not used. Instead, collisions between the two beads of the dimer and the solvent particles were explicitly taken into account in the framework of MPCD. For a single passive dimer, the equilibrium correlation function of force dipoles 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₀[small script l]₀²/(k_BT) ≫ 1, the relaxation time τ₀ = (γk₀)⁻¹ of the dimer should be close to this correlation time. Moreover, we have τ₁ = (γk₁)⁻¹ = τ₀/2. Using such estimates, it can be shown that w₁τ₁ varied between 0.001 and 0.1 in the simulations.²⁰ Because substrate saturation was assumed, conditions w₀τ₀ ≫ 1 and w₁τ₁ ≪ 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.²⁰ They showed damped oscillations and were similar to the correlation function for v₀τ₀ = 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 force-dipole intensity 〈Δm²〉 of active dimers was by about an order of magnitude larger than 〈Δm²〉_eq for the passive ones. Depending on the parameters, it scaled as k^α₀ 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.²⁰ 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.²⁰ In the simulations, the truncated potential

for r < 2^1/24(2r₀) and zero otherwise, with ε = 2.5 k_BT and r₀ = 1.075[small script l]₀, was used to describe steep repulsive interactions between the beads belonging to different dimers. The interaction radius r₀ was chosen as defining the radius of a bead.

Since distances [small script l]₀ and [small script l]₁ = 0.5[small script l]₀ in the open and closed dimer conformations were both smaller than 2r₀ = 2.15[small script l]₀, 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[small script l]₀.

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 [small script l]₀. 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.⁵⁷

When the dimers were active, diffusion of tracers was increasing instead with the dimer volume fraction ϕ. For the most stiff active dimers with k₀[small script l]₀²/(k_BT) = 1440 and the kinetic regime with w₁τ₁ about 0.1, relative diffusion enhancement of D_A/D_T = 0.3 could be observed²⁰ at the dimer volume fraction of ϕ = 0.266. For the least stiff dimers with k₀[small script l]₀²/(k_BT) = 144, diffusion enhancement by 5 percent was seen at ϕ = 0.133.

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²⁰ 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_BT per a turnover cycle). It would be therefore important to undertake such simulations also for the parameters closer approaching those of the real enzymes.

5 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,²⁰ 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¹² 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 all-atom or coarse-grained molecular dynamics simulations for specific 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 eqn (7) and (8) can be approximately sought in the form

p_s(x,t) = π_s(t)p(x,t|s), (49)

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₀ and w₁ 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₀ = ν₀ p(x = [small script l]₀|s = 0), w₁ = ν₁ p(x = [small script l]₁|s = 1). (56)

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²〉. In a steady state, its time derivative is zero. On the other hand, by using eqn (7) and (8) and integrating by parts, we findThus, we straightforwardly obtain that, for an active dimer in any statistically steady state, 〈m〉 = −k_BT.

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₀(x) and p₁(x).

C Force-dipole correlation function

Introducingwe can write the system of two master eqn (7) and (8) concisely aswhereandandand

L̂₀₁ = −u₁(x), L̂₁₀ = −u₀(x). (66)

The general solution of eqn (62) is

where λ_n and q⁽ⁿ⁾ are eigenvalues and eigenvectors of the linear operator L̂,

L̂q⁽ⁿ⁾ = λ_nq⁽ⁿ⁾, (68)

and decomposition coefficients A_n are determined by initial conditions.

Because the master equation must have a stable stationary solution, the operator L̂ should always possess a zero eigenvalue λ₀ = 0 and, furthermore, condition Reλ_n > 0 should hold for all other eigenvalues n.⁶¹ 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 λ₁ ≤ Re λ₂ ≤ Re λ₃ ≤ …). The stationary probability distribution p̄(x) coincides with the eigenvector q⁽⁰⁾(x).

The conditional probability G(x,s,t|x₀,s₀) gives the probability to find the dimer in various states (x,s) at time t provided that it was in the state (x₀,s₀) at time t = 0. It represents a special solution of the master eqn (62) given by

where a_n(x₀,s₀) are the coefficients of decomposition of this initial condition over eigenvectors q⁽ⁿ⁾.

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 stationary state, the correlation function of force dipoles iswhere the complex coefficients B_n are

If we retain in this decomposition only the first, most slowly decaying term, this yields

Therefore, the normalized correlation function iswhere Γ = Re λ₁, Ω = Im λ₁, and B₁ = C(0)e^iα/cos α.

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,⁶² 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.

D Dependence on orientational correlation time

Suppose that the force-dipole correlation function C(t) and the orientational correlation function σ(t) are given by eqn (74) and (39). By taking the integral in eqn (38), we findThis 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 = (Ω − Γ)⁻¹ and we havewhereis the limit of χ when τ_rot ≫ Γ⁻¹ and τ_rot ≫ Ω⁻¹.

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 non-monotonously on the orientational correlation time. It increases linearly with τ_rot at short times, then reaches a maximum at τ_rot = (Ω − Γ)⁻¹ and finally saturates at large orientational correlation times. For example, if we take the values Γ ≈ 1/(2τ₀) and Ω ≈ π/(3τ₀) corresponding to substrate saturation in Fig. 7, the maximum diffusion enhancement would be reached at τ_rot = 1.8τ₀ and, at the maximum, it will be larger by about 30 percent than in the limit τ_rot ≫ τ₀.

Acknowledgements

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

References

1. A. Cressman, Y. Togashi, A. S. Mikhailov and R. Kapral, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 050901.
2. C. Echeverria, Y. Togashi, A. S. Mikhailov and R. Kapral, Phys. Chem. Chem. Phys., 2011, 13, 10527.
3. R. Golestanian and A. Ajdari, Phys. Rev. Lett., 2008, 100, 038101.
4. R. Golestanian and A. Ajdari, J. Phys.: Condens. Matter, 2009, 21, 204104.
5. M. Iima and A. S. Mikhailov, EPL, 2009, 85, 44001.
6. T. Sakaue, R. Kapral and A. S. Mikhailov, Eur. Phys. J. B, 2010, 75, 381.
7. X. Bai and P. Wolynes, J. Chem. Phys., 2015, 143, 165101.
8. H. Diamant, J. Phys. Soc. Jpn., 2009, 78, 041002.
9. M.-J. Huang, A. S. Mikhailov, R. Kapral and H.-Y. Chen, J. Chem. Phys., 2012, 137, 055101.
10. M.-J. Huang, A. S. Mikhailov, R. Kapral and H.-Y. Chen, J. Chem. Phys., 2013, 138, 195101.
11. M.-J. Huang, H.-Y. Chen and A. S. Mikhailov, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 119.
12. A. S. Mikhailov and R. Kapral, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, E3639.
13. R. Kapral and A. S. Mikhailov, Phys. D, 2016, 318–319, 104.
14. A. S. Mikhailov, Y. Koyano and H. Kitahata, J. Phys. Soc. Jpn., 2017, 86, 101013.
15. Y. Koyano, H. Kitahata and A. S. Mikhailov, Phys. Rev. E, 2016, 94, 022416.
16. Y. Hosaka, K. Yasuda, R. Okamoto and S. Komura, Phys. Rev. E, 2017, 95, 052407.
17. K. Yasuda, R. Okamoto, S. Komura and A. S. Mikhailov, EPL, 2017, 117, 38001.
18. K. Yasuda, R. Okamoto and S. Komura, Phys. Rev. E, 2017, 95, 032417.
19. Y. Hosaka, S. Komura and D. Andelman, Phys. Rev. E, 2020, 101, 012610.
20. M. Dennison, R. Kapral and H. Stark, Soft Matter, 2017, 13, 3741.
21. F. Kogler, Interactions of artificial molecular machines, Diploma thesis, Technical Univ. of Berlin, Berlin, 2012.
22. P. Illien, T. Adeleke-Larodo and R. Golestanian, EPL, 2017, 119, 40002.
23. H. Flechsig and A. S. Mikhailov, J. R. Soc., Interface, 2019, 16, 20190244.
24. J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, C. Schütte and F. Noé, J. Chem. Phys., 2010, 134, 174105.
25. D. E. Shaw, P. Maragakis, K. Lindorff-Larsen, S. Piana, R. O. Dror, M. P. Eastwood, J. A. Bank, J. M. Jumper, J. K. Salmon, Y. Shan and W. Wriggers, Science, 2010, 330, 341.
26. H. P. Lu, L. Xun and X. S. Xie, Science, 1998, 282, 1887.
27. H.-P. Lerch, R. Rigler and A. S. Mikhailov, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 10807.
28. M. Gur, E. Zomot and I. Bahar, J. Chem. Phys., 2013, 139, 121912.
29. K. A. Henzler-Wildman, V. Thai, M. Lei, M. Ott, M. Wolf-Watz, T. Fenn, E. Pozharski, M. A. Wilson, G. A. Petsko, M. Karplus, C. G. Hübner and D. Kern, Nature, 2007, 450, 838.
30. Y. Togashi and A. S. Mikhailov, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 8697.
31. H. Flechsig and A. S. Mikhailov, Proc. Natl. Acad. Sci. U. S. A., 2013, 107, 20875.
32. M. Hinczewski, R. Tehver and D. Thirumalai, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, E4059.
33. K. Shiroguchi, H. F. Chin, D. E. Hannemann, E. Muneyuki, E. M. D. L. Cruz and K. Kinosita Jr., Proc. Natl. Acad. Sci. U. S. A., 2011, 9, e1001031.
34. N. Kodera, D. Yamamoto, R. Ishikawa and T. Ando, Nature, 2010, 468, 72.
35. J. M. Berg, L. Stryer, J. Tymoczko and G. Gatto, Biochemistry, WH Freeman, 9th edn, 2019.
36. C. Riedel, R. Gabizon, C. A. M. Wilson, K. Hamadani, K. Tsekouras, S. Marqusee, S. Pressé and C. Bustamante, Nature, 2015, 517, 227.
37. R. Inoue, R. Biehl, T. Rosenkranz, J. Fitter, M. Monkenbusch, A. Radulescu, B. Farago and D. Richter, Biophys. J., 2010, 99, 2309.
38. J. Schonfield, P. Inder and R. Kapral, J. Chem. Phys., 2012, 136, 205101.
39. P. Gouet, H.-M. Jouve, P. A. Williams, I. Andersson, P. Andreoletti, L. Nussaume and J. Hajdu, Nat. Struct. Biol., 1996, 3, 951.
40. S. von Bülow, M. Siggel, M. Linke and G. Hummer, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 9843.
41. Z. Bashardanesh, J. Elf, H. Zhang and D. van der Spoel, ACS Omega, 2019, 4, 20654.
42. G. Nawrocki, A. Karaboga, Y. Sugita and M. Feig, Phys. Chem. Chem. Phys., 2019, 21, 876.
43. P. Illien, X. Zhao, K. Dey, P. J. Butler, A. Sen and R. Golestanian, Nano Lett., 2017, 17, 4415.
44. K. K. Dey, Angew. Chem., Int. Ed., 2019, 58, 2208.
45. A.-Y. Jee, S. Dutta, Y.-K. Cho, T. Tlusty and S. Granick, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 14.
46. M. Xu, J. L. Ross, L. Valdez and A. Sen, Phys. Rev. Lett., 2019, 123, 128101.
47. Y. Zhang, M. J. Armstrong, N. M. B. Kazeruni and H. Hess, Nano Lett., 2018, 18, 8025.
48. H. Wang, M. Park, R. Dong, J. Kim, Y.-K. Cho, T. Tlusty and S. Granick, Science, 2020, 369, 537.
49. M. Guo, A. J. Ehrlicher, M. H. Jensen, M. Renz, J. R. Moore, R. D. Goldman, J. Lippincott-Schwartz, F. C. MacKintosh and D. A. Weitz, Cell, 2014, 158, 822.
50. B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O'Hem, E. R. Dufresne and C. Jacobs-Wagner, Cell, 2014, 156, 183.
51. E. Fodor, M. Guo, N. S. Gov, P. Visco, D. A. Weitz and F. van Wijland, EPL, 2015, 110, 48005.
52. K. Nishizawa, K. Fujiwara, N. Ikenaga, N. Nakajo and D. Mizuno, Sci. Rep., 2017, 7, 15143.
53. S. C. Weber, A. J. Spakowitz and J. A. Theriot, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 7338.
54. R. Bruinsma, A. Y. Grosberg, Y. Rabin and A. Zidovska, Biophys. J., 2014, 106, 1871.
55. A. Vendeville, D. Lariviere and E. Fourmentin, FEMS Microbiol. Rev., 2010, 35, 395.
56. F. C. MacKintosh and A. J. Levine, Phys. Rev. Lett., 2008, 100, 018104.
57. G. L. Hunter and E. R. Weeks, Rep. Prog. Phys., 2012, 75, 066501.
58. Y. Koyano, H. Kitahata and A. S. Mikhailov, EPL, 2019, 128, 40003.
59. N. Oyama, T. Kawasaki, H. Mizuno and A. Ikeda, Phys. Rev. Res., 2019, 1, 032038.
60. R. Kapral, Adv. Chem. Phys., 2008, 140, 89.
61. H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications, Springer, Berlin, 1989.
62. Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer, Berlin, 1984.

---