Spatiotemporal patterns enhanced by intra- and inter-molecular fluctuations in arrays of allosterically regulated enzymes

In our previous paper [6], we assumed a system without noise ($\sigma =0$). Once a cycle starts, a product molecule is always released after a constant delay, and the cycle is finished after a constant duration (cycle time); thus, the system is represented by reaction–diffusion equations with constant delays. We reported spatiotemporal patterns such as global oscillations, traveling waves, and spirals, in which the operation of enzymes was synchronized. When the uniform steady state solution of the equations becomes unstable, small perturbations around the steady state can grow and form the observed patterns; hence, the possibility that a particular pattern would occur could be evaluated by stability analysis of the solution. Such delayed reaction–diffusion schemes have been also applied to other kinds of biochemical reactions [17].

However, when we consider intramolecular fluctuations, i.e., cases with $\sigma \gt 0$, we cannot use the delayed reaction–diffusion equations because these delays are no longer constant. Therefore, to investigate the effects of intramolecular fluctuations, it is necessary to improve the method. Here, in order to perform stability analysis similar to our previous study, we constructed an analytical scheme to include distributed delays in reaction–diffusion equations.

As the analysis depends on the reaction–diffusion equations (i.e., partial differential equations of concentrations without noise) and continuous probability distributions of the delays, it corresponds to a system with a high density of enzymes where stochastic effects (the intermolecular fluctuations introduced later) are negligible.

Let ${{F}_{0}}(x,t)$, ${{F}_{p}}(x,t)$, and ${{F}_{1}}(x,t)$ represent the rates of the initiation of the cycle, product release, and finish of the cycle, respectively, at position x in the space and at time t. For free enzymes with and without regulatory molecules,

$$\begin{eqnarray}&&\frac{\partial {{e}_{1}}}{\partial t}=\beta p{{e}_{0}}-\kappa {{e}_{1}}-{{\alpha }_{1}}{{e}_{1}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\frac{\partial {{e}_{0}}}{\partial t}=-\beta p{{e}_{0}}+\kappa {{e}_{1}}-{{\alpha }_{0}}{{e}_{0}}+{{F}_{1}},\end{eqnarray} \tag{ 3 }$$

where ${{e}_{1}}(x,t)$ and ${{e}_{0}}(x,t)$ are the concentrations of these enzymes, respectively; and for free product concentration $p(x,t)$,

$$\begin{eqnarray}&&\frac{\partial p}{\partial t}=-\beta p{{e}_{0}}+\kappa {{e}_{1}}+{{\alpha }_{1}}{{e}_{1}}-\gamma p+{{F}_{p}}+D\frac{{{\partial }^{2}}p}{\partial {{x}^{2}}}.\end{eqnarray} \tag{ 4 }$$

The rate of the initiation of the cycle is

$$\begin{eqnarray}&&{{F}_{0}}={{\alpha }_{0}}{{e}_{0}}+{{\alpha }_{1}}{{e}_{1}}.\end{eqnarray} \tag{ 5 }$$

For enzymes during the reaction cycle, we consider the process represented by equation (1). F_p is the rate of product release, i.e., arrival at $\phi ={{\phi }_{p}}$ for the first time in the cycle, while F₁ is the rate of arrival at $\phi =1$. Because ϕ does not revert beyond ${{\phi }_{p}}$ after product release, and the cycle is finished once ϕ reaches 1, the cycle is divided into two independent first passage processes: from $\phi =0$ to ${{\phi }_{p}}$ and from $\phi ={{\phi }_{p}}$ to 1.

We can construct a Fokker–Planck type equation corresponding to equation (1) for the density of the enzymes at phase ϕ, and combine that with equations (2)–(4). However, with intramolecular fluctuations ($\sigma \gt 0$), the solution does not take a simple form because of the boundary conditions at $\phi =0$, ${{\phi }_{p}}$, and 1. Therefore, we consider the relationship between the input (initiation of the cycle) and the output (product release or finish of the cycle) instead.

Let $Q(\phi ^{\prime} ,t^{\prime} )$ be the probability distribution function of the first passage time t' for the process represented by equation (1) from $\phi =0$ to $\phi =\phi ^{\prime}$ with reflecting boundary at $\phi =0$; this condition applies to enzymes in a phase section, either $[0,{{\phi }_{p}})$ or $[{{\phi }_{p}},1)$ (with offset ${{\phi }_{p}}$ for the latter). Then,

$$\begin{eqnarray}&&{{F}_{p}}(t)=\int _{0}^{\infty }Q({{\phi }_{p}},t^{\prime} ){{F}_{0}}(t-t^{\prime} ){\rm d}t^{\prime} ,\end{eqnarray} \tag{ 6 }$$

where the rate $Q({{\phi }_{p}},t^{\prime} ){{F}_{0}}(t-t^{\prime} )$ corresponds to enzymes starting the cycle at $t-t^{\prime}$ and releasing a product at time t, i.e., after duration t'. Using Laplace transforms $\tilde{F}$ and $\tilde{Q}$ of F and Q,

$$\begin{eqnarray}&&\tilde{{{F}_{p}}}(s)=\tilde{Q}({{\phi }_{p}},s)\tilde{{{F}_{0}}}(s).\end{eqnarray} \tag{ 7 }$$

Similarly, for F₁ in equation (3), the rate of the finish of the cycle is

$$\begin{eqnarray}\begin{array}{rcl} {{F}_{1}}(t) & = & \int _{0}^{\infty }Q(1-{{\phi }_{p}},t^{\prime\prime} ){{F}_{p}}(t-t^{\prime\prime} ){\rm d}t^{\prime\prime} \\ {} & = & \int _{0}^{\infty }\int _{0}^{\infty }Q({{\phi }_{p}},t^{\prime} )Q(1-{{\phi }_{p}},t^{\prime\prime} ){{F}_{0}}(t-t^{\prime} -t^{\prime\prime} ){\rm d}t^{\prime} {\rm d}t^{\prime\prime} , \\ \end{array}\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{{\tilde{F}}_{1}}(s)=\tilde{Q}({{\phi }_{p}},s)\tilde{Q}(1-{{\phi }_{p}},s){{\tilde{F}}_{0}}(s).\end{eqnarray} \tag{ 9 }$$

Using the derivation in [18], the explicit form of $\tilde{Q}$ for $\sigma \gt 0$ is

$$\begin{eqnarray}&&\tilde{Q}(\phi ^{\prime} ,s)=\frac{R{\rm exp} (\frac{v\phi ^{\prime} }{2\sigma })}{R{\rm cosh} (R\phi ^{\prime} )+\frac{v}{2\sigma }{\rm sinh} (R\phi ^{\prime} )},\end{eqnarray} \tag{ 10 }$$

where $R=\sqrt{{{v}^{2}}+4\sigma s}/2\sigma$. For $\sigma =0$ (ideal delay),

$$\begin{eqnarray}&&\tilde{Q}(\phi ^{\prime} ,s)={\rm exp} (-\frac{\phi ^{\prime} }{v}s).\end{eqnarray} \tag{ 11 }$$

The mean first passage time ${{\tau }_{\phi }}$ is

$$\begin{eqnarray}&&{{\tau }_{\phi }}=\langle t^{\prime} \rangle =-\mathop{{\rm lim} }\limits_{s\to 0}\frac{\partial \tilde{Q}}{\partial s}=\frac{\phi ^{\prime} }{v}-\frac{\sigma }{{{v}^{2}}}\left[ 1-{\rm exp} \left( -\frac{v\phi ^{\prime} }{\sigma } \right) \right],\end{eqnarray} \tag{ 12 }$$

and the variance is

$$\begin{eqnarray}\begin{array}{rcl} \langle {{(t^{\prime} -{{\tau }_{\phi }})}^{2}}\rangle & = & \mathop{{\rm lim} }\limits_{s\to 0}\frac{{{\partial }^{2}}\tilde{Q}}{\partial {{s}^{2}}}-\tau _{\phi }^{2} \\ {} & = & \frac{2\sigma \phi ^{\prime} }{{{v}^{3}}}-\frac{5{{\sigma }^{2}}}{{{v}^{4}}}+(\frac{4\sigma \phi ^{\prime} }{{{v}^{3}}}+\frac{4{{\sigma }^{2}}}{{{v}^{4}}}){\rm exp} (-\frac{v\phi ^{\prime} }{\sigma })+\frac{{{\sigma }^{2}}}{{{v}^{4}}}{\rm exp} (-\frac{2v\phi ^{\prime} }{\sigma }). \\ \end{array}\end{eqnarray} \tag{ 13 }$$

The mean cycle time τ is the sum of ${{\tau }_{\phi }}$ for the two sections before and after ${{\phi }_{p}}$, i.e.,

$$\begin{eqnarray}&&\tau =\frac{1}{v}-\frac{\sigma }{{{v}^{2}}}\left[ 2-{\rm exp} \left( -\frac{v{{\phi }_{p}}}{\sigma } \right)-{\rm exp} \left( -\frac{v(1-{{\phi }_{p}})}{\sigma } \right) \right].\end{eqnarray} \tag{ 14 }$$

Note that c, the total concentration of enzyme, is constant at each position x because the enzymes are not mobile; thus, from equations (2), (3), and (5), $\dot{c}=\dot{{{e}_{0}}}+\dot{{{e}_{1}}}+{{F}_{0}}-{{F}_{1}}=0$ should be satisfied. In the present case, c is also spatially constant as the enzymes are uniformly distributed. The behavior of the system does not change if all of the concentrations are scaled with c, and if β is adjusted so that $\beta c$ is kept unchanged: all the other parameters remain unaffected (note that in this model, only the binding of the product to the regulatory site is second-order, whereas the other reactions are first-order). In the stability analysis below, $\beta c$ is adopted as a parameter.

We first consider concentrations ${{\bar{e}}_{0}}$, ${{\bar{e}}_{1}}$, and $\bar{p}$ at the uniform steady state, satisfying

$$\begin{eqnarray}&&\frac{\partial {{e}_{1}}}{\partial t}=\frac{\partial {{e}_{0}}}{\partial t}=\frac{\partial p}{\partial t}=0.\end{eqnarray} \tag{ 15 }$$

To attain this, the flux should be constant, i.e., from equation (5),

$$\begin{eqnarray}&&{{F}_{1}}={{F}_{p}}={{F}_{0}}={{\alpha }_{0}}{{e}_{0}}+{{\alpha }_{1}}{{e}_{1}}.\end{eqnarray} \tag{ 16 }$$

Using the mean cycle time τ, the total concentration of enzyme is represented by

$$\begin{eqnarray}&&c={{\bar{e}}_{0}}+{{\bar{e}}_{1}}+({{\alpha }_{0}}{{\bar{e}}_{0}}+{{\alpha }_{1}}{{\bar{e}}_{1}})\tau .\end{eqnarray} \tag{ 17 }$$

Substituting equations (15)–(17) to equations (2)–(4), we obtain

$$\begin{eqnarray}&&{{\bar{e}}_{1}}=\frac{c(X-Y+Z)}{V(X+Y+Z)}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{\bar{e}}_{0}}=\frac{2c\gamma W}{X+Y+Z}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\bar{p}=\frac{X-Y+Z}{2\beta \gamma V}\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl} V & = & 1+{{\alpha }_{1}}\tau \\ W & = & \kappa +{{\alpha }_{1}} \\ X & = & c{{\alpha }_{1}}\beta \\ Y & = & \gamma (1+{{\alpha }_{0}}\tau )W \\ Z & = & \sqrt{{{(X-Y)}^{2}}+4c{{\alpha }_{0}}\beta \gamma VW} \\ \end{array}\end{eqnarray*}$$

(although there are two solutions, the other solution (substituting Z with $-Z$ in equations (18)–(20)) does not satisfy $\bar{p}\geqslant 0$).

This uniform steady state always exists. However, if it is unstable, small perturbations may grow into spatiotemporal patterns. To analyze the stability, we consider the deviations $\delta {{e}_{1}}={{e}_{1}}-{{\bar{e}}_{1}}$, $\delta {{e}_{0}}={{e}_{0}}-{{\bar{e}}_{0}}$, and $\delta p=p-\bar{p}$ from the steady state. We require the solutions for equations (2)–(4) in the following forms:

$$\begin{eqnarray}\begin{array}{rcl} \delta {{e}_{1}}(x,t) & = & {{A}_{1}}{\rm exp} (\lambda t-2\pi {\rm i}qx), \\ \delta {{e}_{0}}(x,t) & = & {{A}_{0}}{\rm exp} (\lambda t-2\pi {\rm i}qx), \\ \delta p(x,t) & = & {{A}_{p}}{\rm exp} (\lambda t-2\pi {\rm i}qx), \\ \end{array}\end{eqnarray} \tag{ 21 }$$

where q is the spatial wavenumber, and complex number λ corresponds to the growth rate and temporal (oscillation) frequency. Then, from equations (5) and (7),

$$\begin{eqnarray}&&\tilde{{{F}_{p}}}(s)=\tilde{Q}({{\phi }_{p}},s)\left[ \frac{{{\alpha }_{0}}{{{\bar{e}}}_{0}}+{{\alpha }_{1}}{{{\bar{e}}}_{1}}}{s}+\frac{({{\alpha }_{0}}{{A}_{0}}+{{\alpha }_{1}}{{A}_{1}}){{{\rm e}}^{-2\pi {\rm i}qx}}}{s-\lambda } \right]\end{eqnarray} \tag{ 22 }$$

so that

$$\begin{eqnarray}&&{{F}_{p}}(t)={{\alpha }_{0}}{{\bar{e}}_{0}}+{{\alpha }_{1}}{{\bar{e}}_{1}}+\tilde{Q}({{\phi }_{p}},\lambda )({{\alpha }_{0}}\delta {{e}_{0}}(t)+{{\alpha }_{1}}\delta {{e}_{1}}(t)).\end{eqnarray} \tag{ 23 }$$

Similarly,

$$\begin{eqnarray}&&{{F}_{1}}(t)={{\alpha }_{0}}{{\bar{e}}_{0}}+{{\alpha }_{1}}{{\bar{e}}_{1}}+\tilde{Q}({{\phi }_{p}},\lambda )\tilde{Q}(1-{{\phi }_{p}},\lambda )({{\alpha }_{0}}\delta {{e}_{0}}(t)+{{\alpha }_{1}}\delta {{e}_{1}}(t)).\end{eqnarray} \tag{ 24 }$$

For small $\delta {{e}_{1}}$, $\delta {{e}_{0}}$, and $\delta p$ around the steady state, equations (2)–(4) are linearized as:

$$\begin{eqnarray}&&\frac{\partial }{\partial t}\left( \begin{array}{ccccccccccccccc} \delta {{e}_{1}} \\ \delta {{e}_{0}} \\ \delta p \\ \end{array} \right)={\boldsymbol{ \Lambda }} \left( \begin{array}{ccccccccccccccc} \delta {{e}_{1}} \\ \delta {{e}_{0}} \\ \delta p \\ \end{array} \right),\end{eqnarray} \tag{ 25 }$$

with the linearization matrix

$$\begin{eqnarray}{\boldsymbol{ \Lambda }} =\left( \begin{array}{ccccccccccccccc} -\kappa -{{\alpha }_{1}} & \beta \bar{p} & \beta {{{\bar{e}}}_{0}} \\ \kappa +{{\alpha }_{1}}r & -\beta \bar{p}-{{\alpha }_{0}}(1-r) & -\beta {{{\bar{e}}}_{0}} \\ \kappa +{{\alpha }_{1}}(1+{{r}_{p}}) & -\beta \bar{p}+{{\alpha }_{0}}{{r}_{p}} & -\gamma -\beta {{{\bar{e}}}_{0}}-4{{\pi }^{2}}D{{q}^{2}} \\ \end{array} \right),\end{eqnarray} \tag{ 26 }$$

where ${{r}_{p}}=\tilde{Q}({{\phi }_{p}},\lambda )$ and $r=\tilde{Q}({{\phi }_{p}},\lambda )\tilde{Q}(1-{{\phi }_{p}},\lambda )$. Substituting equations (10), (14), (19), and (20) to equation (26), we obtain the explicit form of ${\boldsymbol{ \Lambda }}$. Here, λ is one of the eigenvalues of ${\boldsymbol{ \Lambda }}$, i.e., satisfying $|{\boldsymbol{ \Lambda }} -\lambda {\bf I}|=0$ (see equations (21) and (25)), whereas ${\boldsymbol{ \Lambda }}$ contains $\tilde{Q}({{\phi }_{p}},\lambda )$ and $\tilde{Q}(1-{{\phi }_{p}},\lambda )$ and, hence, depends on λ. Thus, ${\boldsymbol{ \Lambda }}$ has an infinite number of eigenvalues ${{\lambda }_{i}}$ ($i=1,2,3,\cdots$) that cannot be solved analytically. ${{\mu }_{i}}=\operatorname{Re}({{\lambda }_{i}})$ and ${{f}_{i}}={\rm Im}({{\lambda }_{i}})/(2\pi )$ represent the growth rate (i.e., the instability of the steady state) and frequency of the oscillatory mode, respectively. If ${{\mu }_{i}}$ is positive, the mode can grow when perturbation is applied to the steady state. We sort the modes in ascending order of f_i. The 1st mode (f₁) corresponds to the turnover frequency of each enzyme, and the ith mode ($i\geqslant 2$) has a frequency approximately i times higher. Typically, ${{f}_{i}}\approx iv$ (or slightly lower because the turnover time includes the waiting time between two consecutive cycles and, thus, is longer than the cycle time).

These discrete frequency modes (overtones) arise from the following mechanism. Let us consider a scenario with two groups of enzymes, as shown in figure 2. Each group of enzymes starts the cycle almost simultaneously, and the timing is shifted by approximately half a cycle between the groups. If ${{\phi }_{p}}\approx 0.5$, when one group arrives at $\phi ={{\phi }_{p}}$ and releases product molecules, the other group is approaching, or has just arrived at, $\phi =1$; thus, it can be efficiently activated by the product and start the next cycle. Consequently, such synchronized groups can be maintained, and the oscillations of concentrations are sustained. Although the turnover time of each enzyme cannot be shorter than the cycle time, the apparent oscillation period of the concentrations is halved, and this is the basis of the mode for which frequency is twice as high. Similarly, the ith mode corresponds to a scenario with i groups of enzymes in the phase space. We call this the i-group mode.

Zoom In Zoom Out Reset image size

As predicted (and already reported for $\sigma =0$ [7, 13]), the i-group mode is likely to grow if ${{\phi }_{p}}$ is close to a multiple of $1/i$; e.g., the 3-group mode may appear for ${{\phi }_{p}}\approx 1/3$ or $2/3$. Synchronization is also possible for ${{\phi }_{p}}$ close to 0, in which case a fraction of enzymes start the cycle spontaneously and soon release products, which in turn activate the remaining enzymes: this is typical in the 1-group mode. For ${{\phi }_{p}}\approx 1$, activation by products released by the same group of enzymes is possible, which also enables synchronization. The condition for synchronization is also affected by the diffusion of the product if we consider $q\gt 0$ (spatial waves), as discussed in the section 3.

In the following sections, we solve the eigenvalue problem numerically and self-consistently by iteratively converging (Newton–Raphson method) from different initial values of λ, corresponding to different modes.

When the density of enzymes is low, stochastic effects are not negligible. To consider such cases, we also adopted stochastic reaction–diffusion simulations. For simplicity, we assumed 1-dimensional space; in 2-dimensional space, target patterns, spiral waves, and wave turbulence appear [6, 7], but the basic properties discussed later are not affected.

The space is divided into L_x bins, and the enzymes are randomly distributed over the space and fixed. We set the volume of each bin to 1, there are ${{e}_{1}}(x)$ and ${{e}_{0}}(x)$ free enzyme molecules with and without regulatory products, and $p(x)$ product molecules in the xth bin. The total number of enzyme molecules is $c{{L}_{x}}$. The product molecules diffuse, with the diffusion coefficient D [bin²/unit-time]. Additionally, a no-flux condition is applied. Binding of the product to the regulatory site of the enzyme in the same bin occurs at a rate of $\beta p{{e}_{0}}$ [/unit-time], and dissociation from the site occurs at $\kappa {{e}_{1}}$ [/unit-time]. We assume the substrate is abundant and its concentration is kept constant. Each free enzyme starts a cycle at rate ${{\alpha }_{1}}$ (active enzyme) or ${{\alpha }_{0}}$ (inactive enzyme) [/unit-time]. For each enzyme in the reaction cycle, the progress of the phase represented by equation (1) is simulated, with release of a product molecule at $\phi ={{\phi }_{p}}$ and return to the free state at $\phi =1$. ϕ does not revert beyond the boundaries at $\phi =0$ and $\phi ={{\phi }_{p}}$. Upon dissociation from the regulatory site and release from the catalytic site of the enzyme, the product molecule is placed in the same bin. Initially (t = 0), all the enzymes are free and no product exists.

Numerically, the progress of the enzyme phases and the diffusion of the product are processed at constant time intervals, whereas stochastic binding and dissociation at the regulatory site, initiation of a reaction cycle, and decay of the product are executed in each bin by the next reaction method [19].

First, we show the uniform (well-mixed) case, which is similar to a previous study [14]. Bifurcation diagrams for the uniform case (zero wavenumber) for different strengths σ of intramolecular fluctuations are shown in figure 3. The upper boundary of β is drastically decreased by the intramolecular fluctuations, whereas the lower boundary is minimally affected. This is because if β is large, a small amount of the product released earlier via intramolecular fluctuations is sufficient to prematurely activate a significant proportion of free enzymes.

Zoom In Zoom Out Reset image size

When we consider non-uniform cases, the growth rate of each mode also depends on the spatial wavenumber. In the current case, only the product can diffuse; therefore, the wavenumber q can be scaled by the diffusion coefficient D without loss of generality, resulting in the normalized wavenumber $q^{\prime} =\sqrt{D}q$ (see equation (26)). We calculated q' at which the growth rate takes the maximum value in $q^{\prime} \geqslant 0$ (dominant wavenumber), as shown in figure 4 (the maximum growth rate and the frequency at the dominant wavenumber are also shown).

Zoom In Zoom Out Reset image size

For some ${{\phi }_{p}}$ (e.g., 1-group mode with ${{\phi }_{p}}=0.3$), the dominant wavenumber is 0 and growth rate is positive, which means that global synchronization is strongest among the 1-group modes. For other cases, e.g., 1- and 3-group modes with ${{\phi }_{p}}=0.455$, the dominant wavenumber is not zero (also see figure 5). Traveling or standing waves are expected to have a non-zero wavenumber.

Zoom In Zoom Out Reset image size

As previously reported [6, 7], wave bifurcation is possible even outside of the Hopf bifurcation (uniform oscillation) boundary shown in figure 3. In such cases, the growth rate is negative for $q^{\prime} =0$ but positive for a range of q'. When $\sigma =0$, the 1-group uniform mode vanishes at ${{\phi }_{p}}=0.409$, whereas 1-group waves are possible up to ${{\phi }_{p}}=0.480$ with large wavenumbers ($q^{\prime} \approx 1.5$). For traveling waves, because the cycle time is almost constant, a larger wavenumber corresponds to slower propagation or to a longer delay from neighboring sites; thus, it is reasonable to suggest that such modes allow a communication delay through larger ${{\phi }_{p}}$ than uniform oscillations.

For larger σ, the maximum growth rate is lower, and synchronization is possible only in smaller parameter regions, as shown in figure 4(b); however, the dominant wavenumber and frequency do not change significantly (figures 4(a) and (c)).

Predictably, synchronization with many groups, i.e., at a high frequency, is sensitive to intramolecular fluctuations. If σ is large, only a small number of synchronized groups are allowed. However, 1-group synchronization is not always the most fluctuation-tolerant mode, especially in spatially non-uniform cases. For example, in the case shown in figure 5, the 2-group uniform (q = 0) oscillations are more robust against intramolecular fluctuations than the 1-group wave patterns with a finite wavenumber. Such an effect is indeed seen in numerical simulations; 2-group uniform oscillations are dominant in a certain range of σ, as shown in figure 6 (center column).

Zoom In Zoom Out Reset image size

We have considered fluctuations that originate from the internal states of the enzymes; however, there is another potential source of fluctuations: the stochastic binding, dissociation, and decay of molecules. In a biological cell, the concentration of an enzyme species is not always high. If molecular events such as collision, binding, and release involve a low concentration of enzymes, they are inevitably stochastic and show large fluctuations in their rate.

Because we can change the enzyme concentration c without affecting the behavior in the continuum limit if $\beta c$ is kept unchanged, we can investigate the effects of stochasticity. In figures 6 and 7, spatiotemporal patterns for two different enzyme concentrations are shown. At a low concentration, synchronization between the enzymes is more stochastic and often disturbed. Of course, at a concentration much lower than one enzyme per the diffusion length of the product (i.e., the distance the product diffuses before decaying), enzymes can no longer communicate and synchronization is lost. However, such stochasticity (variation in the waiting time) may even strengthen certain synchronization modes, in contrast to the intramolecular fluctuations (variation in the cycle time), which merely hinder synchronization.

Zoom In Zoom Out Reset image size

In some parameter regions, two or more modes (frequencies) are possible. With such parameters, stochastic fluctuations affect the selection of modes: an example is shown in figure 7. In this case, the parameters allow mixed modes. Where the concentration of enzymes is high (see left panel in figure 7), synchronization with 1 group occurs first, while it grows very slowly with 3 or 5 groups. In contrast, at a low enzyme concentration (right panel in figure 7), synchronization with 3 or 5 groups is induced much more quickly by the stochastic noise. This behavior is similar to noise-enhanced phase synchronization [20]; however, the frequency depends on the intrinsic dynamics of the enzymes.

It is also important to understand whether synchronization waves can be induced in a system where synchronization is not observed without fluctuations. Near the boundary of parameters for synchronization, once the enzymes are synchronized, oscillatory patterns may sustain for some cycles. If fluctuations are present, weakly oscillating patterns may appear. On some occasions, complete synchronization can almost be recovered even if β is so high that synchronization is impossible in the continuum limit. For this type of synchronization, the discreteness of the enzymes plays an important role.

In the continuum limit, the concentrations can be infinitesimally small but non-zero. Even in synchronized states, there is always a small fraction of free enzymes that can start a cycle. Although the decay of the product is fast, a small fraction of the regulatory product may remain after cycling and activate the enzymes asynchronously. If β is too high, this residual product can gradually desynchronize the enzymes. However, with discrete molecules, the numbers of the product and free enzyme molecules should be integers and reach 0 at some point, thereby hindering such step-by-step desynchronization and sustaining the synchronized state.

The situation in this stochastic synchronization is analogous to discreteness-induced transitions, previously reported by one of the authors [10, 21]. No matter how fast the positive feedback reaction progresses, it completely stalls when the relevant molecules vanish. In the previously reported case [10, 21], for the transitions to occur, decay of a certain chemical species had to be faster than the interval of molecular inflow of the chemical; however, in the present case, decay must be sufficiently faster than the turnover time, although the threshold is also affected by β.

Intramolecular fluctuations decrease the upper limit of β for synchronization in the continuum limit, as shown in figure 3. In contrast, intermolecular fluctuations may induce synchronization even if β is larger than the limit, i.e., they increase the effective upper limit of β. For synchronization with 2 or more groups, these effects are significant. In spatially distributed cases, the boundary of ${{\phi }_{p}}$ and β for synchronization also depends on the wavelength. Hence, the combination of intra- and intermolecular fluctuations may determine the synchronization mode (frequency and wavelength) that appears.

Figure 6 shows an example of such effects observed in stochastic simulations. Without any fluctuations, 1-group standing-wave patterns and 2-group global synchronization are possible (figure 5), although the parameter values are close to the bifurcation boundary. Mixed behavior is indeed observed with $\sigma =0$ and c = 1000 (figure 6, left column). The 1-group standing-wave patterns are enhanced under stronger intermolecular fluctuations (c = 10 with $\beta c$ fixed; right column), whereas the 2-group global synchronization dominates under moderate intramolecular fluctuations ($\sigma =5\times {{10}^{-6}}$; center column).