Synchronization and enhanced catalysis of mechanically coupled enzymes

We examine the stochastic dynamics of two enzymes that are mechanically coupled to each other, e.g., through an elastic substrate or a fluid medium. The enzymes undergo conformational changes during their catalytic cycle, which itself is driven by stochastic steps along a biased chemical free energy landscape. We find conditions under which the enzymes can synchronize their catalytic steps, and discover that the coupling can lead to a significant enhancement in their overall catalytic rate. Both effects can be understood as arising from a global bifurcation in the underlying dynamical system at sufficiently strong coupling. Our findings suggest that, despite their molecular scale, enzymes can be cooperative and improve their performance in metabolic clusters.

Introduction.-Since the observation of "the sympathy of two clocks" by Christiaan Huygens in 1665, synchronization phenomena have been observed in a variety of systems, at different time-and length-scales [1][2][3]. Generic theoretical frameworks, in particular the Kuramoto model, have been widely used to predict the conditions under which synchronization can occur [4]. However, the individual oscillators are usually coupled to each other through a physical medium, and it is often necessary to include the microscopic details of the mechanical coupling between them in order to obtain a useful description. One well-studied example is the synchronization of periodically beating flagella and cilia [5,6], for which hydrodynamic interactions have been shown to play a crucial role [7][8][9][10][11][12]. Another example is synchronization mediated by elastic stresses in a solid substrate, known to be important for algal flagella [13,14] and cardiac cells [15,16]. In these micro-scale examples, the cyclic motion of each individual oscillator is driven by a non-vanishing, deterministic driving force.
At the even smaller nanoscale, however, molecular motors and enzymes convert chemical energy into mechanical work following repeated thermodynamic cycles [17][18][19][20][21]. These processes take place in a noise-activated regime, where motion only occurs stochastically, in response to barrier-crossing events along a chemical free energy landscape. Recently, the relationship between the conformational changes of enzymes during their catalytic cycle [22,23] and their translational dynamics has been the subject of many experimental and theoretical studies [24][25][26][27][28][29][30][31][32][33]. While these studies have considered the effect of enzymatic activity on mechanical motion, how mechanical interactions feed back into enzymatic activity and whether synchronization of the catalytic cycles across enzymes is possible are questions that remain open. These questions are of high biological relevance considering that enzymes are frequently assembled into clusters [34][35][36][37][38][39].
In this Letter, we study the dynamics of two enzymes which undergo conformational changes during their cat-alytic cycle and interact with each other mechanically. We show that this coupling is sufficient to synchronize their stochastic catalytic steps, and moreover leads to a significant enhancement of their catalytic rate.
Mechanochemical coupling.-We use a minimal model where each enzyme α is considered to have a single mechanical degree of freedom L α , which might represent for example its elongation, see Fig. 1(a-c); and a single chemical degree of freedom or phase φ α , which is a reaction coordinate describing the state of the chemical reaction happening inside the enzyme. Both L α and φ α evolve together according to the potential Here, the first term describes conformational changes of the enzyme during the catalytic cycle, with L(φ α ) being the rest length of the enzyme as a function of the reaction coordinate and k the stiffness of the enzyme, while the second term represents The catalytic cycle of each enzyme is represented by a phase φα evolving in a biased free energy landscape V (φα) (solid black). The enzyme elongation Lα tries to adapt to a phasedependent rest length L(φα) (dotted red). the free energy of the reaction, described by a biased potential V (φ α ) that drives the phase forward [17,18]; see Fig. 1(d). Assuming that the overdamped medium surrounding the enzymes couples forces to velocities linearly, the dynamics of the elongation of enzyme α will be governed byL α = µ f α +hf β , where µ is the mobility associated to the elongation, and h is the (dimensionless) mechanical coupling between the two enzymes α and β = α.
As an example, for enzymes directly coupled into a complex as in Fig. 1(a) the coupling constant h can be easily calculated and is found to be negative with 0 > h > −1 [40]. The internal forces (force-dipoles) f α and f β generated by the corresponding enzymes can be calculated as In turn, the phase dynamics are given byφ where µ φ is the mobility along the chemical reaction coordinate. Phase equations.-The coupled dynamics of length and phase can be simplified further by projecting the dynamics of the lengths, assumed fast, onto the slow manifold of the configuration space defined by the phases φ α [40].

Staying to lowest order in
, which corresponds to the assumption that phase changes rather than conformational changes constitute the bottleneck in the dynamics of the enzyme [41], the deterministic dynamics for the phases readṡ with ω(φ α ) ≡ −µ φ V (φ α ). This shows that the reaction dynamics of two enzymes that undergo conformational changes during catalysis are coupled through the mechanical interaction. We note that, in contrast to usual models of synchronization, the interaction term is proportional to the driving force ω(φ α ). In fact, whereas in other models the coupling term can be understood as coming from the gradient of a potential (e.g. a term proportional to − cos(φ α −φ β ) in the Kuramoto model), such a description is not possible here: the velocity field corresponding to the right hand side of Eq. (1) has non-zero curl when h = 0. Moreover, while in typical descriptions of synchronization the driving force ω(φ α ) is either a constant or a positive-definite function of the phase that can be mapped onto a constant using a gauge transformation [4], in our system ω(φ α ) vanishes and changes sign twice through a complete catalytic cycle, because of the energy barrier in V (φ). This highlights that catalysis is an activated process that can only occur in the presence of noise. Stochastic dynamics.-Thermal noise can be systematically added to the deterministic dynamics in Eq. (1) [40], resulting in the stochastic dynamicṡ where the mobility tensor is defined as M 11 = M 22 = 1, and M 12 = M 21 = hA 1 A 2 /(1 − h 2 ), the square-root of the mobility tensor Σ is defined via M αβ = Σ αν Σ βν [40], and Einstein summation convention for repeated indices is used. Eq. (2) is to be interpreted in the Stratonovich sense. The first term is identical to the deterministic dynamics in Eq. (1), the last term is a (multiplicative) noise where ξ satisfies ξ α (t)ξ β (t ) = δ αβ δ(t − t ), and the second term is the spurious drift associated with this multiplicative noise. These dynamics are constructed such that they correspond to the Fokker-Planck independently of the value of h, whenever the system allows equilibration, e.g. when the potential V (φ) is unbiased or the range of φ 1 , φ 2 is bounded. Eq. (2) highlights that the forces on the phases are actually conservative, with associated potential V (φ 1 ) + V (φ 2 ), but are connected to the phase velocities via a non-diagonal mobility matrix, with the coupling constant h setting the magnitude of the off-diagonal terms, which become important when the system is out of equilibrium. As a side-note, we believe that these physically-suggestive properties make Eq.
(2) stand out in the context of synchronization (as an unexplored class of dynamical systems), not just in the particular form of the off-diagonal terms that we have derived here from mechanical coupling, but in its full generality [42].
Brownian dynamics simulations.-In what follows, we model the reaction free energy using a washboard potential V (φ) = −F φ − v cos [φ + arcsin(F/v)], which has minima at φ = 2πn for all n ∈ Z when F/v < 1.
Here, F and v determine the height of the free energy barrier E ba and the free energy difference of the chemical reaction E * [see Fig. 1 The washboard potential thus represents an unending sequence of substrate-to-product transformation reactions, assuming that substrate is abundant and that the binding and unbinding of substrate and product in between each reaction happen very fast compared to all other processes. The conformational changes of the enzyme mimic the washboard potential, with the rest length given by L(φ) = L 0 + cos [φ + arcsin(F/v)], so that the extrema of V (φ) coincide with those of L (φ). Here, represents the amplitude of the conformational changes and may be positive or negative depending on whether the enzyme expands or contracts during catalysis. The synchronization dynamics, however, is independent of the sign of , as the equations of motion are invariant under changes of this sign. Defining dimensionless time as τ ≡ tµ φ v, the system depends only on three dimensionless parameters: the rescaled coupling constant scape as determined by E ba /E * , and the noise strength Results of Brownian dynamics simulations of Eq. (2) are shown in Fig. 2 and Fig. S1 in the Supplemental Material [40]. As expected, the enzymes undergo stochastic, quasi-discrete steps. While for low or positive couplingh the two enzymes mostly do individual catalytic steps, for sufficiently negativeh the two enzymes tend to step in synchrony; see Fig. 2(a,b). Remarkably, for sufficiently low noise (k B T /E ba 1), we observe long synchronized runs, in which the two enzymes undergo a large number of joint catalytic steps after a thermal fluctuation kicks them out of a local free energy minimum (and before falling back into a minimum). An example of a 5-step run can be seen between the two black arrows for the red trajectory in Fig. 2(a).
As a quantitative measure of synchronization, we use the phase-difference diffusion coefficient D ∆ , calculated from the relation (φ 1 − φ 2 ) 2 ∼ 2D ∆ t. This can be compared to the single-phase diffusion coefficient D φ calculated from (φ α − φ α ) 2 ∼ 2D φ t. If φ 1 and φ 2 were independent variables, we would expect D ∆ /D φ = 2, as we observe forh = 0. However, negative values ofh lead to D ∆ /D φ < 2, implying that synchronous steps are favored; see Fig. 2(c,e). Synchronization is most pro-nounced for strong bias (E ba /E * 1) and low noise (k B T /E ba 1). For positiveh we find predominantly D ∆ /D φ 2, implying in this case that synchronous steps are inhibited. We then consider the total catalytic activity, i.e. the number of catalytic steps per enzyme per unit time over the whole time of the simulation Ω ≡ [φ 1 (τ tot ) + φ 2 (τ tot )]/(4πτ tot ); see Fig. 2(d,f). We find that mechanical coupling can enhance catalytic activity, particularly for strongly synchronized cases (with strong bias, low noise, and negativeh), in which case enhancements as large as 200% are seen, a remarkable observation for just two coupled enzymes.
Phase portrait.-The fact that synchronization and enhanced catalysis are strongest at low noise suggests that their emergence may be understood from the phase portraits of the underlying deterministic dynamical system, Eq. (1). Indeed, the dynamical system, which is defined on the torus, undergoes a global bifurcation for sufficiently strong negative couplingh; see Fig. 3(a,b). Forh >h * , the phase space is divided into four basins of attraction, separated by four heteroclinic orbits. At a critical valueh =h * , the heteroclinic orbits change topology, and forh <h * two of the heteroclinic orbits become two homoclinic orbits, between which a running band of periodic orbits emerges. This topological bifurcation has important consequences. Whereas forh >h * thermal fluctuations will typically kick a system which is initially at (0, 0) into either of the basins of attraction (2π, 0) or (0, 2π) corresponding to a single-enzyme step, forh <h * the system is instead kicked into either the basin of attraction (2π, 2π) corresponding to a synchronized two-enzyme step, or into the running band. In the latter case, the system will perform multiple synchronized steps until noise kicks it out of the band once again. The critical valueh * decreases with increasing E ba /E * , as seen in Fig. 3(c), and appears to be well described byh * −2(E ba /E * ) 1/4 . This value is plotted as the dashed line in Fig. 2(c-f), and correctly predicts the regions with enhanced catalysis and synchronization.
Discussion.-Using a minimal model, we have shown that enzymes that undergo conformational changes during their catalytic cycle can synchronize with each other through mechanical interactions, which moreover can significantly enhance their overall catalytic rate. These effects are favored for negative mechanical coupling h < 0, which implies that contraction of one enzyme favors expansion of the other, and vice versa. A negative coupling is guaranteed for complexed enzymes as in Fig. 1(a) [40], and should be expected in similar configurations such as the one in Fig. 1(c). Here, synchronization arises as an entrainment of the inherently stochastic, noise-activated catalytic steps of the two enzymes. While synchronization in excitable systems has been described before [43], particularly in the context of FitzHugh-Nagumo oscillators [44][45][46], the coupling in these systems was by means of an added coupling force (diffusive or Kuramoto-like) and resulted in a synchronization transition mediated by standard Hopf or saddle-node bifurcations. In contrast, in our system we find a novel form of coupling which arises from off-diagonal terms in the mobility matrix that connects forces to velocities, and thus leaves the equilibrium probability distribution of the system intact while introducing non-trivial effects in an out-of-equilibrium setting. The resulting transition is mediated by a global bifurcation which gives not only synchronization but also enhanced catalysis. We note that the mechanism for enhanced catalysis that we observe is also distinct from recent proposals for activated barrier crossing [47][48][49] which rely on colored noise.
Biological enzymes often form dense assemblies where mechanical interactions should be expected. These range from large-scale, three-dimensional clusters such as "metabolons" [34] and bacterial micro-compartments [35], to ordered filaments such as the cytoophidium [36], to oligomeric complexes of just a few enzymes [37]. While the functional benefit of these structures is not yet clear, it has been proposed that proximity favors channeling of reaction intermediates among different enzymes in the same catalytic pathway [34]. Our result of enhanced catalysis provides a possible additional advantage to close proximity between enzymes. In fact, many enzymes that are functional in their monomeric form but also assemble into homooligomeric forms [see Fig. 1(a)] are more catalytically active in the oligomeric form [37], a behavior which could be explained by mechanical coupling as proposed here. Moreover, synchronization effects could be particularly relevant in the context of rapid and robust signaling by membrane ion channels, which also operate in clusters [38,39]. Future experiments may also test our predictions in a controlled in vitro setting, by creating enzyme assemblies with designed geometry [50][51][52], or using single-molecule techniques that allow for the measurement of individual catalytic events and conformational changes [53][54][55][56].
This work has received support from the Max Planck School Matter to Life and the MaxSynBio Consortium, which are jointly funded by the Federal Ministry of Education and Research (BMBF) of Germany, and the Max Planck Society. T.A.-L. acknowledges the support of an EPSRC Studentship.
shorter than that of changes in the phase velocity. The dynamics of the deviation of the length from its preferred value δL α ≡ L α − L(φ α ) is governed by (S12) Assuming fast relaxation, we can set δL α = 0 and solve for δL α to obtain This shows that there is a difference between the elastic rest length of the enzyme and its dynamical rest length. The non-vanishing deformation leads to a finite force (-dipole) that is present during the enzyme activity. Introducing Eq. (S13) into Eq. (S11), and solving forφ α , we obtaiṅ the effective stochastic dynamics. To this end, starting from Eq. (S14), we build the stochastic dynamics of the system in form of the corresponding Fokker-Planck equation for the probability distribution P(φ 1 , φ 2 ; t), which ensures equilibration to the Boltzmann distribution P(φ 1 , φ 2 ) ∝ exp{−[V (φ 1 ) + V (φ 2 )]/k B T } when equilibration is possible. Equation (S23) belongs to the general class of stochastic systems with multiplicative noise, which implies that building a Langevin dynamics for the system requires special considerations. We first introduce the square-root of the mobility tensor Σ via M αβ = Σ αν Σ βν , and then construct the corresponding Langevin equations [associated with Eq. (S23)] asφ where ξ satisfies ξ α (t)ξ β (t ) = δ αβ δ(t − t ). In Eq. (S24), which coincides with Eq. (2) in the main text, the Stratonovich convention for the multiplicative noise is used and a corresponding spurious drift term is introduced. This completes the stochastic formulation of the system of two mechanically coupled enzymes. The square-root Σ of the mobility tensor M can be calculated explicitly in the two limiting cases for M described in the previous section. When the phase constitutes the bottleneck of the dynamics and the mobility matrix is given by (S18-S19), we find Σ 11 = Σ 22 = ν + + ν − and Σ 12 = Σ 21 = ν + − ν − , where we have defined ν ± ≡ (1/2) 1 ± hA 1 A 2 /(1 − h 2 ). When h is small and the mobility matrix is given by (S20-S22) we find, again to first order in h, Σ 11 = 1/ 1 +  S2. The same bifurcation as described in the main text occurs when neither the conformational changes nor the phase changes constitute a bottleneck. Here, we have set µ φ 2 /µ = 1 and used the full mobility matrix as given by (S15-S17). As the coupling is increased from (a) h = −0.6 to (b) h = −0.8, the bifurcation is crossed and a running band of periodic orbits emerges. The bias of the free energy landscape is set to E ba /E * = 10 −3 . The same bifurcation as described in the main text occurs in a simpler model, given by Eq. (S14) but with a constant mobility matrix M11 = M22 = 1 and M12 = M21 = g. As the coupling is increased from (a) g = 0.2 to (b) g = 0.4, the bifurcation is crossed and a running band of periodic orbits emerges. The bias of the free energy landscape is set to E ba /E * = 0.02.