Coherence of biochemical oscillations is bounded by driving force and network topology

Biochemical oscillations are prevalent in living organisms. Systems with a small number of constituents cannot sustain coherent oscillations for an indefinite time because of fluctuations in the period of oscillation. We show that the number of coherent oscillations that quantifies the precision of the oscillator is universally bounded by the thermodynamic force that drives the system out of equilibrium and by the topology of the underlying biochemical network of states. Our results are valid for arbitrary Markov processes, which are commonly used to model biochemical reactions. We apply our results to a model for a single KaiC protein and to an activator-inhibitor model that consists of several molecules. From a mathematical perspective, based on strong numerical evidence, we conjecture a universal constraint relating the imaginary and real parts of the first nontrivial eigenvalue of a stochastic matrix.

Figure 1

Correlation function $C_{1,1}\left(t\right)$. The number of states is $N=100,\mathcal{A}=200$, and the transition rates are uniform with $k^{-}=1$. For this case $X_I\simeq 0.401$, which gives a period of 15.66, and $X_R\simeq 0.01655$, as indicated by the red solid line.

Figure 2

Multicyclic network. (a) Network of states, with the three-state cycle marked with a magenta circle, the four-state cycle with the dashed red links, and the five-state cycle with the solid blue links. (b) Numerically evaluated $\mathcal{R}$ for randomly chosen rates against the bound (solid red line). The points were generated according to the method explained in Appendix pp3.

Figure 3

Model and results for a single KaiC protein. (a) For the vertical arrows, the transition rates are $\gamma e^{\mathit{\Delta }\mu /2}$ for the larger arrow and $\gamma e^{E/6}$ for the small arrows. For the horizontal arrows, the transition rate from $C_i$ to ${\tilde{C}}_i$ is $k{e}^{\chi E(i-3)/3}$ and the transition rate from ${\tilde{C}}_i$ to $C_i$ is $k{e}^{\tilde{\chi }E(3-i)/3}$, where $\chi$ ($\tilde{\chi }$) is an indicator function that is zero (one) for $i=0,1,2,3$ and one (zero) for $i=4,5,6$. (b) Ratio $\mathcal{R}$ as a function of $\mathit{\Delta }\mu$. The dots were obtained by numerical maximization of $\mathcal{R}$ with respect to the parameters $\gamma$ and $E$, where $k=1$. The dotted red line represents $\mathcal{R}$ for $k=1,\gamma =e^5$, and $E=10$. The blue solid line is the bound $f(\mathcal{A},N)$ for $\mathcal{A}=6\mathit{\Delta }\mu$ and $N=14$.

Figure 4

The activator-inhibitor model. (a) Scheme representing the chemical reactions in the model, where a blue arrow represents activation and the red lines with a square at the end represent inhibition. The four forms of the enzyme are free enzyme $M$, the species $R$ bound to the enzyme $MR$, phosphorylated enzyme $M_p$, and phosphorylated enzyme bound to a phosphatase $M_{p}K$. (b) Numerical results for the ratio $\mathcal{R}$. The number of enzymes is $N_M=500$. The points were obtained from numerical simulation as explained in Appendix pp6 and the solid lines represent the estimate of the $\mathcal{R}$ given by $N_K\mathit{\Delta }\mu /\left(2\pi \right)$. The optimal number of oscillations takes place close to $N_K=30$, with $N_K=20$ there were practically no oscillations and with $N_K=50$ the number of oscillations gets considerably smaller.

Figure 5

Scatter plot for the unicyclic model with $N=4$. The bound in Eq. (6) is represented by the solid red line. The rate $k_1^{+}$ was set to $k_1^{+}=e^A{k}_1^{-}{k}_2^{-}{k}_3^{-}/\left(k_2^{+}{k}_3^{+}\right)$, the other seven rates were randomly chosen as ${10}^x$, with $x$ uniformly distributed between $-3$ and 3. For this figure, we have evaluated $\mathcal{R}$ for ${10}^7$ sets of rates.

Figure 6

First example of a multicyclic network. (a) Network of states. (b) Numerical evidence of the bound. The dashed red line is the function $f(3\mathrm{\Delta }\mu /2,3)$ and the solid red line the function $f(\mathrm{\Delta }\mu /2,4)$. The large blue dots were obtained with numerical maximization of $\mathcal{R}$ and the black dots were obtained by evaluating $R$ at randomly generated transition rates. We have generated two sets of ${10}^7$ points each. For one set we have chosen eight independent rates as $k_{24},k_{43}$, and $k_{31}$ were set to ${\mathrm{e}}^{\mathit{\Delta }\mu /8};k_{21},k_{34}$, and $k_{42}$ were set to ${10}^x;k_{14}$ was set to $e^{-5\mathit{\Delta }\mu /4}{10}^{-5}{10}^x$ and $k_{41}$ as set to ${10}^{-5}{10}^x$, where $x$ is uniformly distributed between $-3$ and 3. The remaining rates $k_{12}$ and $k_{13}$ were determined by the relations (C1) and (C2), respectively. For the other set, the eight independent rates were chosen as ${10}^x$.

Figure 7

Second example of a multicyclic network. (a) Network of states. (b) Numerical evidence of the bound. The solid red line is the function $f(\mathrm{\Delta }\mu ,4)$. The large blue dots were obtained with numerical maximization of $\mathcal{R}$ and the black dots were obtained by evaluating $R$ at randomly generated transition rates. We have generated ${10}^7$ points, choosing nine independent transition rates as ${10}^x$, with $x$ uniformly distributed between $-2$ and 2. The remaining transition rates $k_{12},k_{13}$, and $k_{42}$ were determined by the affinities in Eqs. (C4), (C5), and (C6).

Figure 8

Third example of a multicyclic network. (a) Network of states. (b) Numerical evidence of the bound. The solid red line is the function $f(\mathrm{\Delta }\mu ,6)$. The large blue dots were obtained with numerical maximization of $\mathcal{R}$ and the black dots were obtained by evaluating $R$ at randomly generated transition rates. We have generated ${10}^7$ points, choosing 12 independent transition rates as ${10}^x$, with $x$ uniformly distributed between $-2$ and 2. The two remaining transitions, $k_{25}$ and $k_{36}$, were determined by the affinities in Eqs. (C9) and (C10).

Figure 9

Fourth example of a multicyclic network. (a) Network of states. (b) Numerical evidence of the bound. The solid red line is the function $f(\mathrm{\Delta }\mu ,5)$. The large blue dots were obtained with numerical maximization of $\mathcal{R}$ and the black dots were obtained by evaluating $R$ at randomly generated transition rates. We have generated ${10}^7$ points, choosing 11 independent transition rates as ${10}^x$, with $x$ uniformly distributed between $-2$ and 2. The three remaining transition rates, $k_{51},k_{52}$, and $k_{42}$, were determined by the affinities in Eqs. (C12), (C13), and (C14).

Figure 10

Phosphorylation level of the KaiC protein. The parameters of the model are set to $k=e^{0.48097},\gamma =e^{10.3475},E=19.86$, and $\mathcal{A}=20$. Numerical calculation of the first eigenvalue gives $X_R\simeq 6.82\times {10}^7$ and $X_I\simeq 29.8\times {10}^7$, which gives $\mathcal{R}\simeq 4.37$. The red solid line represents an exponential function with decay exponent $X_R$.

Figure 11

Relation between energy dissipation and number coherent oscillations. (a) Unicyclic model with $N=3,k_1^{-}=k_2^{-}=k_3^{-}=1,k_1^{+}=k_3+=e^{\mathcal{A}/4}$, and $k_2^{+}=e^{\mathcal{A}/2}$. (b) Model for a single KaiC with $k=1,\gamma =e^5$, and $E=10$.

Figure 12

Numerical simulation of the activator-inhibitor model. Parameters are set to $\mathit{\Delta }\mu =12$ and $N_K=30$. The variable $n$ in (b) labels the the peak in (a). The solid red line in (b) is an exponential fit to the data that gives $\mathcal{R}\simeq 2\pi /0.12301\simeq 51.1$.

Figure 13

Stochastic trajectories of the activator-inhibitor model. The black trajectory shows the number of enzymes in state $M$, the red in state $MR$, the green in state $M_p$, and the blue in state $M_{p}K$. Parameters are set to $\mathit{\Delta }\mu =20$ and for (a) $N_K=30$ and (b) $N_K=150$.