Coherence of Biochemical Oscillations is Bounded by Driving Force and Network Topology

TL;DR

This paper demonstrates that the maximum number of coherent biochemical oscillations is universally limited by the system's thermodynamic driving force and network topology, with broad implications for biological oscillators.

Contribution

It establishes a universal bound on biochemical oscillation coherence based on thermodynamics and network structure, applicable to arbitrary Markov processes.

Findings

Bound on coherence increases with thermodynamic force

Network topology constrains oscillation precision

Universal eigenvalue relation conjectured for stochastic matrices

Abstract

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 non-trivial eigenvalue of a stochastic matrix.

Equations

$$E_1\underset{k_2^{-}}{\overset{k_1^{+}}{\rightleftharpoons }}{E}_2\underset{k_3^{-}}{\overset{k_2^{+}}{\rightleftharpoons }}{E}_3\dots E_{N-1}\underset{k_N^{-}}{\overset{k_{N-1}^{+}}{\rightleftharpoons }}{E}_N\underset{k_1^{-}}{\overset{k_N^{+}}{\rightleftharpoons }}{E}_1,$$

$$\mathcal{A}\equiv \ln \prod _{i=1}^N{k}_i^{+}/k_i^{-},$$

$${\mathbf{L}}_{j,i}\equiv k_i^{+}{\delta }_{i,j-1}+k_i^{-}{\delta }_{i,j+1}-(k_i^{-}+k_i^{+}){\delta }_{i,j},$$

$$C_{1,1}(t)\equiv {\left[\exp (\mathbf{L}t)\mathbf{P}(0)\right]}_1,$$

$$\mathcal{R}\equiv X_I/X_R,$$

$$\mathcal{R}\le \cot (\pi /N)\tanh [\mathcal{A}/(2N)]\equiv f(\mathcal{A},N).$$

$$\mathcal{R}\le f({\mathcal{A}}^{\ast },N^{\ast })\le {\mathcal{A}}^{\ast }/(2\pi ).$$

$${\mathbf{L}}_{j,i}\equiv k_i^{+}{\delta }_{i,j-1}+k_i^{-}{\delta }_{i,j+1}-(k_i^{-}+k_i^{+}){\delta }_{i,j},$$

$$\mathcal{R}\equiv X_I/X_R=\sqrt{4C_1/C_2^2-1},$$

$$C_1\equiv$$

$$+k_1^{-}{k}_2^{+}+k_2^{-}{k}_3^{+}+k_3^{-}{k}_1^{+}$$

$$C_2\equiv k_1^{+}+k_2^{+}+k_3^{+}+k_1^{-}+k_2^{-}+k_3^{-}.$$

$$\mathcal{A}\equiv \ln \left(\prod _{i=1}^N{k}_i^{+}/k_i^{-}\right).$$

$$\Lambda (\{k_i^{+}\},\{k_i^{-}\},\alpha )=4C_1/C_2^2-\alpha (k_1^{+}{k}_2^{+}{k}_3^{+}-k_1^{-}{k}_2^{-}{k}_3^{-}{\text{e}}^{\mathcal{A}}),$$

$$\frac{d\Lambda }{d{k}_1^{+}}=\frac{4(k_2^{+}+k_3^{+}+k_3^{-})}{C_2^2}-8\frac{C_1}{C_2^3}-k_2^{+}{k}_3^{+}\alpha$$

$$\frac{d\Lambda }{d{k}_1^{-}}=\frac{4(k_2^{-}+k_2^{+}+k_3^{-})}{C_2^2}-8\frac{C_1}{C_2^3}+{\text{e}}^{\mathcal{A}}{k}_2^{-}{k}_3^{-}\alpha .$$

$$\alpha =\frac{4{\text{e}}^{-2\mathcal{A}/3}({\text{e}}^{\mathcal{A}/3}-1)}{9({\text{e}}^{\mathcal{A}/3}+1{)}^3(k^{-}{)}^3}.$$

$${\mathcal{A}}_1\equiv \ln \frac{k_{12}{k}_{24}{k}_{41}}{k_{21}{k}_{42}{k}_{14}},$$

$${\mathcal{A}}_2\equiv \ln \frac{k_{13}{k}_{34}{k}_{41}}{k_{31}{k}_{43}{k}_{14}}.$$

$${\mathcal{A}}_3={\mathcal{A}}_1-{\mathcal{A}}_2.$$

$${\mathcal{A}}_1=\ln \frac{k_{12}{k}_{23}{k}_{34}{k}_{41}}{k_{21}{k}_{32}{k}_{43}{k}_{14}}=\Delta \mu ,$$

$${\mathcal{A}}_2=\ln \frac{k_{12}{k}_{23}{k}_{31}}{k_{21}{k}_{32}{k}_{13}}=0,$$

$${\mathcal{A}}_3=\ln \frac{k_{12}{k}_{24}{k}_{41}}{k_{21}{k}_{42}{k}_{14}}=0,$$

$${\mathcal{A}}_4=\ln \frac{k_{13}{k}_{34}{k}_{41}}{k_{31}{k}_{43}{k}_{14}}=\Delta \mu ,$$

$${\mathcal{A}}_5=\ln \frac{k_{23}{k}_{34}{k}_{42}}{k_{32}{k}_{43}{k}_{24}}=\Delta \mu .$$

$${\mathcal{A}}_1=\ln \frac{k_{12}{k}_{23}{k}_{36}{k}_{65}{k}_{54}{k}_{41}}{k_{21}{k}_{32}{k}_{63}{k}_{56}{k}_{45}{k}_{14}}=\Delta \mu .$$

$${\mathcal{A}}_2=\ln \frac{k_{12}{k}_{25}{k}_{54}{k}_{41}}{k_{21}{k}_{52}{k}_{45}{k}_{14}}=\Delta \mu .$$

$${\mathcal{A}}_3=\ln \frac{k_{23}{k}_{36}{k}_{65}{k}_{52}}{k_{32}{k}_{63}{k}_{56}{k}_{25}}=0.$$

$${\mathcal{A}}_1=\ln \frac{k_{12}{k}_{23}{k}_{34}{k}_{45}{k}_{51}}{k_{21}{k}_{32}{k}_{43}{k}_{54}{k}_{15}}=\Delta \mu .$$

$${\mathcal{A}}_2=\ln \frac{k_{23}{k}_{34}{k}_{45}{k}_{52}}{k_{32}{k}_{43}{k}_{54}{k}_{25}}=\Delta \mu ,$$

$${\mathcal{A}}_3=\ln k_{12}{k}_{24}{k}_{45}{k}_{51}{k}_{21}{k}_{42}{k}_{54}{k}_{15}=0.$$

Full text

I Introduction

Circadian rhythms Goldbeter (1997), the cell cycle Ferrell et al. (2011) and gene expression in somitogenesis lewi03 constitute examples of biochemical oscillations that are of central importance for the functioning of living systems. While older observations of biochemical oscillations were made with glycosis Pye and Chance (1966), more recent advances include the observation of 24-h oscillations of the phosphorylation level of the Kai proteins that form the circadian clock of a cyanobacterium Nakajima et al. (2005); Dong and Golden (2008). Synthetically engineered genetic circuits can also display oscillatory behavior Potvin-Trottier et al. (2016). On the theoretical side, the basic conditions for biochemical oscillations to set in are well understood for deterministic rate equations that ignore fluctuations Novák and Tyson (2008). Such rate equations correspond to an effective description of an underlying set of chemical reactions that is fully described by a stochastic chemical master equation.

In principle, biochemical oscillations can happen in a small system with large fluctuations in the chemical species that oscillates, leading to variability in the period of oscillations. Hence, stochastic biochemical oscillations cannot be coherent for an indefinite time. The number of coherent oscillations, which quantifies the precision of the biochemical oscillator, is given by the time for which oscillations remain coherent divided by the period of oscillation more07 ; Cao et al. (2015). In such a context, a relevant question is as follows: Given a biochemical system with significant fluctuations, what is the number of coherent oscillations that can be sustained?

In a recent work related to this question, Cao et al. Cao et al. (2015) have investigated several stochastic models that display biochemical oscillations. They demonstrated that this number of coherent of oscillations increases with a larger rate of entropy production, which quantifies the free energy consumption of the biochemical system. Their work can be seen as part of the recent effort to understand the relation between a certain kind of precision and free energy consumption in biological systems, which include studies on kinetic proofreading Qian (2007), adaptation Lan et al. (2012), cellular sensing Mehta and Schwab (2012); Barato et al. (2013); De Palo and Endres (2013); Lang et al. (2014); Govern and ten Wolde (2014a, b); Hartich et al. (2015), information processing Sartori et al. (2014); Barato et al. (2014); Bo et al. (2015); Ito and Sagawa (2015); mcgr17 , and cost of precision in Brownian clocks Barato and Seifert (2016). In particular, we have recently shown a general relation that establishes the minimal energetic cost for a certain precision associated with a random variable like the output of a chemical reaction. This thermodynamic uncertainty relation Barato and Seifert (2015a); piet16 ; Gingrich et al. (2016) can be used to infer an unknown enzymatic scheme in single molecule experiments Barato and Seifert (2015b) and yields a bound on the efficiency of a molecular motor Pietzonka et al. (2016).

In this paper, we obtain a universal bound on the number of coherent oscillations that can be sustained in any biochemical system that can be modelled as a Markov process with discrete states. This universal bound depends on the thermodynamic forces that drive the system out of equilibrium and on the topology of the network of states. Our results are derived from a conjecture about the first non-trivial eigenvalue of a stochastic matrix that we support with thorough numerical evidence. Specifically, we obtain a bound on the ratio of the imaginary and real parts of this eigenvalue that quantifies the number of coherent oscillations. We illustrate our results with a model for a single KaiC protein van Zon et al. (2007) and with an activator-inhibitor model with several molecules Cao et al. (2015).

The paper is organized as follows. We consider the simple case of a unicyclic network in Sec. II. Our general bound for arbitrary multicyclic networks is formulated in Sec. III. In Sec. IV we apply our results to the two models. We conclude in Sec. V. In Appendix A we provide evidence for the bound for the case of unicyclic networks. The relation between the number of coherent oscillations and the Fano factor is discussed in Appendix B. Numerical evidence for our conjecture is presented in Appendix C. Appendix D is dedicated to the model for a single KaiC. The relation between number of coherent oscillation and the entropy production in analyzed in Appendix E. Finally, Appendix F is dedicated to the activator-inhibitor model.

II Unicycic network

As a simple model for a biochemical oscillation we start with a single enzyme $E$ with the unicyclic reaction scheme

[TABLE]

where $k_{i}^{\pm}$ are transition rates. A generic transition from state $E_{i}$ to $E_{i+1}$ can represent, for example, a conformational change, binding of substrate to the enzyme or the release of a product from the enzyme. The thermodynamic force driving this system out of equilibrium is given by the affinity Seifert (2012)

[TABLE]

where Boltzmann’s constant $k_{B}$ multiplied by the temperature $T$ is set to $k_{B}T=1$ throughout in this paper. For example, if one $ATP$ is consumed and $ADP+P_{i}$ generated in the cycle in Eq. (1), then the affinity is the chemical potential difference $\mathcal{A}=\mu_{ATP}-\mu_{ADP}-\mu_{P_{i}}$.

The model from Eq. (1) follows the master equation $d\mathbf{P}(t)/dt=\mathbf{L}\mathbf{P}(t)$, where $\mathbf{P}(t)=\{P_{1}(t),P_{2}(t),\ldots,P_{N}(t)\}^{T}$ is the vector of probabilities to be in a certain state. The stochastic matrix $\mathbf{L}$ is defined by

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta, $j-1=N$ for $j=1$, and $j+1=1$ for $j=N$. Let us assume that the enzyme is phosphorylated only in state $E_{1}$. The precision of oscillations in the phosphorylation level of an enzyme that is phosphorylated at time $t=0$ is characterized by the number of coherent oscillations in the correlation function $C_{1,1}(t)$ plotted in Fig. 1, which is the probability that the enzyme is in state $E_{1}$ at time $t$ given that the enzyme was in state $E_{1}$ at time [math], i.e.,

[TABLE]

where $\mathbf{P}(0)=\{1,0,0,\ldots,0\}$ and the subscript $1$ indicates the first component of the vector $\exp(\mathbf{L}t)\mathbf{P}(0)$. For large $t$, this correlation function tends to $P^{st}_{1}$, which is the stationary distribution for state $1$. This stationary distribution is the right eigenvector of the stochastic matrix $\mathbf{L}$ that is associated with the eigenvalue [math].

The first nontrivial eigenvalue of the stochastic matrix $\lambda=-X_{R}\pm X_{I}i$, gives the decay time $X_{R}^{-1}$ and the period of oscillations $2\pi/X_{I}$ in Fig. 1. We characterize the coherence of oscillations by the ratio Qian and Qian (2000)

[TABLE]

where the number of coherent oscillations more07 ; Cao et al. (2015) is $X_{I}/(2\pi X_{R})=\mathcal{R}/(2\pi)$ .

For general Markov processes that fulfill detailed balance, which corresponds to $\mathcal{A}=0$ for the unicyclic model, $X_{I}=0$ and there are no oscillations in correlation functions. Hence, a non-zero driving affinity $\mathcal{A}$ is a necessary condition for biochemical oscillations. In particular, for the case of uniform rates in Eq. (1) given by $k_{i}^{-}=k^{-}$ and $k_{i}^{+}=k^{-}\textrm{e}^{\mathcal{A}/N}$, we obtain $X_{R}=[1-\cos(2\pi/N)](k_{+}+k_{-})$ and $X_{I}=\sin(2\pi/N)(k_{+}-k_{-})$.

For the general unicyclic scheme in Eq. (1) with fixed affinity $\mathcal{A}$ and number of states $N$, the ratio $\mathcal{R}$ is maximized for uniform transition rates, which leads to our first main result

[TABLE]

Thus, the maximal number of coherent oscillations in a unicyclic network is bounded by the thermodynamic force $\mathcal{A}$ and by the network topology through the number of states $N$. The evidence for this bound is as follows. For $N=3$ we can show analytically that uniform rates correspond to a maximum of $\mathcal{R}$, whereas for larger $N$ we rely on extensive numerical evidence as shown in Appendix A. Specifically, we have confirmed this conjecture up to $N=8$ with both numerical maximization of $\mathcal{R}$ and evaluation of $\mathcal{R}$ at randomly chosen rates. Similar to the ratio $\mathcal{R}$, the Fano factor associated with the probability current is extremized for uniform rates Barato and Seifert (2015a, b). However, as discussed in Appendix B, the bound in Eq. (6) and this earlier bound on the Fano factor are different results, i.e., one does not imply the other.

III Multicyclic networks

Biochemical networks are typically more complicated than a single cycle. We now extend the bound from Eq. (6) to general multicyclic networks. As an example, we consider an enzyme $E$ that consumes a substrate $S$ and generates a product $P$. The enzyme has two binding sites, leading to the network of states shown in Fig. 2, which is a common model in enzyme kinetics Barato and Seifert (2015b). The affinity that drives the system out of equilibrium is the chemical potential difference between substrate and product $\varDelta\mu=\mu_{S}-\mu_{P}$.

This network of states has four types of cycles, as shown in Fig. 2. There are cycles with three states and affinity $\varDelta\mu$, like the cycle $E+S\to ES\to EP\to E+P$; cycles with four states and affinity [math], like the cycle $E+S+P\to ES+P\to ESP\to EP+S\to E+S+P$; cycles with five states and affinity $\varDelta\mu$, like the cycle $ES+S\to ESS\to ESP\to EPP\to EP+P\to ES+P$; and one cycle with six states and affinity $2\varDelta\mu$, which is the outer cycle in Fig. 2 that goes through all states. Among all these cycles, the last one with $\mathcal{A}=2\varDelta\mu$ and $N=6$ leads to the largest value of the function $f(\mathcal{A},N)$. We have verified numerically that indeed $f(2\varDelta\mu,6)$ bounds the ratio $\mathcal{R}$ with numerical maximization and numerical evaluation at randomly chosen rates, as shown in Fig. 2. The bound is saturated if the transition rates for the cycle with six states are uniform and much faster than the rates associated with the three links in the middle that are not part of the six-state cycle. In this way, the multicyclic network corresponds effectively to a unicyclic network with six states. In Appendix C, we perform similar numerical tests for several multicyclic networks that do not share any symmetry, and in all cases the ratio $\mathcal{R}$ follows a similar bound.

Based on this numerical evidence we conjecture the following universal bound on the ratio $\mathcal{R}$. Consider an arbitrary Markov process with a finite number of states $N$ on an arbitrary multicyclic network. The cycles are labeled by $\alpha$, with a number of states $N_{\alpha}\leq N$ and affinity $\mathcal{A}_{\alpha}$, where $\textrm{e}^{\mathcal{A}_{\alpha}}$ is the product of forward transition rates divided by backward transition rates over all links in the cycle (see Appendix C). The affinity and number of states of the cycle with the maximal value of $f(\mathcal{A}_{\alpha},N_{\alpha})$, defined in Eq. (6), are denoted by $\mathcal{A}^{*}$ and $N^{*}$, respectively, i.e., $f(\mathcal{A}^{*},N^{*})=\textrm{max}_{\alpha}f(\mathcal{A}_{\alpha},N_{\alpha})$. The ratio $\mathcal{R}$ is then bounded by

[TABLE]

The basic idea behind this bound is that the simple unicyclic network in Eq. (1) is a building block for a generic multicyclic network: any two point correlation function cannot have a larger number of coherent oscillations than the bound determined by its “best” cycle. Hence, the number of coherent oscillations is bounded by the thermodynamic force $\mathcal{A}^{*}$ and the topology of the network of states, as characterized by $N^{*}$. Our bound in Eq. (7) leads to two general necessary conditions for a large number of coherent oscillations, a large number of states and a large maximal affinity. For biochemical models with irreversible transitions, e.g., the models in more07 ; van Zon et al. (2007), the affinity $\mathcal{A}^{*}$ formally diverges, and the bound in Eq. (7) becomes $\mathcal{R}\leq\cot(\pi/N^{*})\leq\cot(\pi/N)$. The weaker second inequality involving the total number of states $N\geq N^{*}$ follows from a known result about the eigenvalues of a discrete time stochastic matrix dmit45 ; dmit46 ; swif72 . For the case of a complex network of states where identifying the large number of states in a cycle is not feasible, like in the activator-inhibitor model below, we can use the second inequality in Eq. (7), based on $\lim_{N\to\infty}f(\mathcal{A},N)=\mathcal{A}/(2\pi)$. We now proceed to illustrate this second main result in two models.

IV Case studies

IV.1 Model for a single KaiC

First, we consider a model for a single KaiC hexamer along the lines of the model proposed in van Zon et al. (2007). The assumptions entering the model, which is depicted in Fig. 3, are the following. A phosphate can bind to each one of the six monomers, hence the phosphorylation level of the hexamer varies from $i=0$, with no phosphate, to $i=6$, with all monomers phosphorylated. Each of the six monomers can be either active or inactive. However, either all monomers are active or all monomers are inactive, since the energetic cost of having two monomers with different conformations is high enough to avoid such configurations. There are a total of $14$ states, denoted by $C_{i}$ for $i$ phosphorylated active monomers and $\tilde{C}_{i}$ for $i$ phosphorylated inactive monomers. If the hexamer is active, phosphorylation reactions occur and if the hexamer is inactive only dephosphorylation reactions occur.

The transition rates of this model for a single KaiC protein are given in the caption of Fig. 3. The parameter $\varDelta\mu$ is the chemical potential difference of ATP hydrolysis. The parameter $E$ sets the energy of a state that depends on the hexamer activity and on the phosphorylation level. If the hexamer is active the energy of a state $C_{i}$ is $Ei/6$ and if it is inactive the energy of an state $\tilde{C}_{i}$ is $E(6-i)/6$. This parametrization implies that the transition rate from $C_{6}$ to $\tilde{C}_{6}$ and the transition rate from $\tilde{C}_{0}$ to $C_{0}$ are both larger than the rates for the respective reversed transitions. The parameters $k$ and $\gamma$ are related the time-scales of changes in the phosphorylation level and conformational changes between active and inactive, respectively.

The phosphorylation level of the KaiC protein oscillates with the number of coherent oscillations given by $\mathcal{R}/(2\pi)$, as shown in in Appendix D. The cycle with the largest value of the function $f(\mathcal{A},N)$ is the cycle that goes through all $N=14$ states with $\mathcal{A}=6\varDelta\mu$, which is marked with the red arrows in Fig. 3. Hence, for this model we have $\mathcal{R}\leq f(6\varDelta\mu,14)$, as shown in Fig. 3.

For fixed $E$ and $\gamma$ we obtain the red dashed curve in Fig. 3 for $\mathcal{R}$ as a function of $\varDelta\mu$. Interestingly, while the number of coherent oscillations has a maximum, after which it decreases to zero with increasing $\varDelta\mu$, the entropy production from stochastic thermodynamics Seifert (2012) is an increasing function of $\varDelta\mu$. Hence, the number of coherent oscillations can also decrease with an increase of the rate of free energy consumption, which provides a counter example to the relation between the number of coherent oscillations and energy dissipation inferred in Cao et al. (2015). We discuss the relation between $\mathcal{R}$ and the entropy production further in Appendix E.

The maximal number of coherent oscillations $\mathcal{R}/(2\pi)$ that can be achieved in this model is strictly speaking less than $1$. If a single molecule does not have a large number of states, several coherent oscillations can only be sustained in a system with many molecules as we discuss next in our second example.

IV.2 Activator-inhibitor model

We consider the activator-inhibitor model from Cao et al. (2015), see Fig. 4. The main components of this model are inhibitors $X$, activators $R$ and enzymes $M$ that can be in four different states. The enzyme goes through a phosphorylation cycle over these four states, hydrolyzing one ATP thus liberating a free energy $\varDelta\mu$. The enzyme $M$ in its phosphorylated form ($M_{p}$) activates $R$ and $X$, whereas $X$ inhibits $R$. Furthermore, the enzyme $M$ must bind an $R$ in order to phosphorylate. Hence, $R$ activates the production of $R$ and $X$, while $X$ inhibits $R$. This feedback loop leads to oscillations in, for example, the number of species $X$. Finally, there is a phosphatase $K$ that must bind to the enzyme $M$ for the dephosphorylation reaction. Further details of the model are given in Appendix F.

Two important aspects about the behavior of this model are the following. First, the number of oscillations increases with $\varDelta\mu$ and saturates for large enough $\varDelta\mu$. Second, in order for different enzymes $M$ to synchronize their cycles, they must compete for the smaller number of phosphatase $N_{K}<N_{M_{T}}$, where $N_{M_{T}}$ is the total number of enzymes, as explained in Appendix F. However, if $N_{K}$ is too small, the number of enzymes that synchronize is also too small to generate oscillations. Hence, there is an optimal value for $N_{K}$. These features are shown Fig. 4.

Due to the complex network of states of this model we bound the ratio $\mathcal{R}$ with the second inequality in Eq. (7). The largest affinity $\mathcal{A}^{*}$ is given by $\mathcal{A}^{*}=\varDelta\mu N_{M_{T}}$, which corresponds to a cycle where all enzymes $M$ go through their cycles in a synchronized way.

As shown in Fig. 4, the values of $\mathcal{R}$ obtained with numerical simulations are approximately one order of magnitude below the fundamental limit set by our bound of $N_{M_{T}}\varDelta\mu/(2\pi)$, which gives $\mathcal{R}\simeq 796$ for $\varDelta\mu=10$. This result is reasonable, as saturating the bound in a multicyclic network requires transition rates such that an optimal cycle dominates, which is not the case for the present model. The realization of this optimal cycle in a stochastic trajectory would require an unlikely sequence of events that all enzymes $M$ go through their own cycles in a synchronized way.

For the case of a close to optimal value of $N_{K}$ that maximizes $\mathcal{R}$, the number of enzymes $M$ that synchronize is roughly $N_{K}$. Hence, cycles with an affinity $\varDelta\mu N_{K}$ should be typical. Guided by our bound, it is then interesting to compare the value of $\mathcal{R}$ with the estimate $\varDelta\mu N_{K}/(2\pi)$. As shown in Fig. 4, indeed this number can give an approximate value of $\mathcal{R}$, typically overestimating $\mathcal{R}$. The estimate is best for $N_{K}=30$ and for $\varDelta\mu=10$. For large $\varDelta\mu$ the bound grows linearly, while $\mathcal{R}$ saturates. Our results then indicate that this estimate works well in the regime where $\varDelta\mu$ is close to its saturation value and $N_{K}$ is close to its optimal value. Even though this heuristic argument is restricted to this model, if competition for some scarce molecule that lead to synchronization is present in the biochemical system, e.g., in the model for the Kai system from van Zon et al. (2007), a similar reasoning that leads to an estimate of the number of coherent oscillations should be valid.

V Conclusion

In summary, we have conjectured a new bound on the number of coherent biochemical oscillations for systems with large fluctuations. Our universal result depends only on the thermodynamic forces that drive the system out equilibrium and on the network topology through the cycle with the largest function $f$ from Eq. (6). Knowledge of the chemical potential differences and the network of states of a biochemical system thus leads to a bound on the number of coherent oscillations. As illustrative examples, we obtained the largest number of oscillations that can be sustained by a single KaiC hexamer and analyzed the activator-inhibitor model, showing that our bound is also valid in models with a number of molecules that is large enough to make the network of states complicated but small enough to keep fluctuations relevant and, therefore, make a description in terms of deterministic rate equations inappropriate.

It remains to be seen whether and how our bound can be used as a guiding principle to understand how systems like circadian clocks have evolved and to engineer systems with precise oscillations in synthetic biology. Our results apply to autonomous biochemical oscillators. Analyzing the relation between precision and thermodynamics for biochemical oscillators that are coupled to an external periodic signal is an interesting direction for future work. Finally, a rigorous mathematical proof of our conjecture about the first excited eigenvalue of a stochastic matrix is an open problem for the theory of Markov processes.

Acknowledgements.

We thank F. Jülicher for helpful discussions.

Appendix A Evidence for unicyclic network

We discuss the evidence for the bound in Eq. (6) for the unicyclic scheme. The mathematical problem is to calculate the first non-trivial eigenvalue of the stochastic matrix

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta, $j-1=N$ for $j=1$, and $j+1=1$ for $j=N$. The absolute value of the imaginary (real) part of this eigenvalue is denoted $X_{I}$ ($X_{R}$).

For the case $N=3$ this eigenvalue can be exactly calculated with some algebra, leading to

[TABLE]

where

[TABLE]

and

[TABLE]

If $4C_{1}<C_{2}^{2}$, there are no oscillations in correlations functions. We want to find the transition rates that maximize $\mathcal{R}$ for a fixed affinity

[TABLE]

This maximum can be found with the Lagrange function

[TABLE]

where $\alpha$ is a Lagrange multiplier. The derivatives of $\Lambda$ with respect to the transition rates are given by

[TABLE]

and

[TABLE]

Due to symmetry, it is easy to deduce the derivatives with respect to $k_{2}^{\pm}$ and $k_{3}^{\pm}$ from the above expressions. If we substitute uniform rates $k_{i}^{-}=k^{-}$ and $k_{i}^{+}=\textrm{e}^{\mathcal{A}/3}k_{-}$ in the above expressions, we obtain that both derivatives become zero with a Lagrange multiplier

[TABLE]

Hence, we have proved that $\mathcal{R}$ is extremized for uniform rates. We can easily evaluate $\mathcal{R}$ for specific rates and check that uniform rates indeed correspond to a maximum.

For larger $N$, up to $N=8$, we have calculated this eigenvalue numerically. We have maximized the ratio $\mathcal{R}$ numerically and observed that it is maximized for uniform rates in all cases, providing convincing evidence for the bound in Eq. (6). As an independent check we have also evaluated $\mathcal{R}$ numerically for randomly chosen rates. As an example, we show a scatter plot obtained with this method in Fig. 5.

Appendix B Relation between $\mathcal{R}$ and the Fano factor

We now explain the difference between the bound in Eq. (6) and a bound on the Fano factor $F$ obtained in Barato and Seifert (2015a, b). This Fano factor is given by $F=2D/J$, where $J$ is the average probability current and $D$ the diffusion constant associated with the current Barato and Seifert (2015a). This bound on the Fano factor can be written as $F^{-1}\leq N\tanh[\mathcal{A}/(2N)]$, where the quantity $F^{-1}$ is also maximized for uniform rates. Furthermore, for uniform rates $2D=(k_{+}+k_{-})/N^{2}$ and $J=(k_{+}-k_{-})/N$, implying $X_{I}=\sin(2\pi/N)NJ$ and $X_{R}=[1-\cos(2\pi/N)]N^{2}2D$. Nevertheless, for arbitrary transition rates there is no such simple relations, with a prefactor that only depends on $N$, between $X_{I}$ ($X_{R}$) and $J$ ($D$).

For uniform rates $\mathcal{R}=F^{-1}N^{-1}\cot(\pi/N)$. Since $\mathcal{R}$ can be zero even out of equilibrium and $F^{-1}$ becomes zero only in equilibrium, we know that $F^{-1}N^{-1}\cot(\pi/N)$ can be larger than $\mathcal{R}$. Evaluating $\mathcal{R}$ and $F$ at different rates we find that $\mathcal{R}$ can also be larger than $F^{-1}N^{-1}\cot(\pi/N)$. Hence, the bound on the Fano factor from Barato and Seifert (2015a, b) does not imply our result in Eq. (6). Their main similarity is that both the Fano factor $F$ and the ratio $\mathcal{R}$ are extremized for uniform rates: it is common to find a function of several variables that is extremized at a symmetric point.

Appendix C Evidence for multicyclic networks

In this appendix we explain the numerical evidence for the bound in Eq. (7) for multicyclic networks. For all cases, we have confirmed our bound with numerical calculation of the first nontrivial eigenvalue of the stochastic matrix. We have confirmed the bound with both numerical maximization of $\mathcal{R}$ and by evaluating $\mathcal{R}$ at randomly chosen transition rates.

As a first example of a multicyclic network, we consider the network with four states shown in Fig. 6. The numbers represent states and the links between them represent transition rates that are not zero. A transition rate from state $i$ to $j$ is denoted $k_{ij}$. This network has three cycles: two cycles with three states $\mathcal{C}_{1}=(1,2,4,1)$ and $\mathcal{C}_{2}=(1,3,4,1)$, and one cycle with four states $\mathcal{C}_{3}=(1,2,3,4,1)$. The affinity of cycle $\mathcal{C}_{1}$ is

[TABLE]

and the affinity of cycle $\mathcal{C}_{2}$ is

[TABLE]

The affinity of $\mathcal{C}_{3}$, which is not independent of $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$, can be written as

[TABLE]

For the results shown in Fig. 6, which confirm our bound for this network, we have set the affinities of the cycles as $\mathcal{A}_{1}=3\varDelta\mu/2$ and $\mathcal{A}_{2}=\varDelta\mu$, which leads to $\mathcal{A}_{3}=\varDelta\mu/2$. As shown in Fig. 6, the cycle with the largest value for the function $f(\mathcal{A}_{\alpha},N_{\alpha})$, where $N_{\alpha}$ is the number of states of cycle $\alpha$, depends on the value of $\varDelta\mu$. For $\varDelta\mu<10.38$ the cycle with largest value for this function is $\mathcal{C}_{1}$ with $f(3\varDelta\mu/2,3)$. For $\varDelta\mu>10.38$ the cycle with largest value for this function is $\mathcal{C}_{3}$ with $f(\varDelta\mu/2,4)$.

The second multicyclic network has four states and is fully connected, as shown in Fig. 7. For this network we have one cycle four states and four cycles with three states. We fix the affinities of these cycles as

[TABLE]

[TABLE]

[TABLE]

where from now on we define the cycles through their affinities. These three affinities determine the values of the two remaining affinities as

[TABLE]

[TABLE]

The cycle leading to the maximal value of $f(\mathcal{A}_{\alpha},N_{\alpha})$ is the cycle with four states and affinity $\mathcal{A}_{1}=\varDelta\mu$. The numerical results illustrating the bound $\mathcal{R}\leq f(\varDelta\mu,4)$ for this network is shown in Fig. 7.

Our third example is the network with six states and three cycles shown in Fig. 8. The affinity of the cycle with six states is fixed as

[TABLE]

We also fix the affinity

[TABLE]

These two affinities determine the affinity of the third cycle as

[TABLE]

The dominant cycle is the one with affinity $\mathcal{A}_{1}=\varDelta\mu$ and six states. The numerical evidence for the bound $\mathcal{R}\leq f(\varDelta\mu,6)$ is shown if Fig. 8

The fourth and last example is the network with five states and five cycles shown in Fig. 9. The affinity of the five states cycle is set to

[TABLE]

There are two four states cycles with affinities

[TABLE]

and

[TABLE]

These three affinities determine the affinity of the two remaining three states cycles,which are

[TABLE]

and

[TABLE]

The dominant cycle has affinity $\mathcal{A}_{1}=\varDelta\mu$ and five states. The numerical evidence for the bound $\mathcal{R}\leq f(\Delta\mu,5)$ is shown in Fig. 9.

For the multicyclic network given in Fig. 2 we have performed a similar analysis. This network has a total of 11 cycles. We identify the states $E,ES,EP,ESS,ESP,EPP$ in Fig. 2 as $1,2,3,4,5,6$, respectively. We have generated two sets of points. The first set has $10^{9}$ points and we just accepted results fulfilling $\mathcal{R}\geq\tanh(\varDelta\mu/8)$. For this set we have chosen the rates $k_{12}$, $k_{31}$, $k_{45}$, and $k_{56}$ as $\textrm{e}^{\varDelta\mu/6}10^{(2+x)/2}$; the rates $k_{21}$, $k_{25}$, $k_{32}$, $k_{35}$, $k_{36}$, $k_{42}$, $k_{52}$, $k_{53}$, $k_{54}$, $k_{56}$, and $k_{65}$ as $10^{x}$; with $x$ uniformly distributed between $-2$ and $2$. The four remaining rates $k_{13}$, $k_{23}$, $k_{24}$, and $k_{63}$ were determined by the constraints set by the affinities. The second set has $10^{7}$ points and the $14$ independent rates were chosen as $10^{x}$, with $x$ uniformly distributed between $-2$ and $2$.

In summary, we have confirmed our bound numerically for four different networks aside from the one displayed in Fig. 2. Since these networks do not share any symmetry, our full numerics provides strong evidence for our bound conjectured in Eq. (7).

Appendix D Phosphorylation level in the model for a single KaiC

We show that the phosphorylation level of the KaiC protein displays oscillations, with the number of coherent oscillations characterized by the ratio $\mathcal{R}$.

For the single KaiC model we analyze in Sec. IV.1, the stochastic matrix reads

[TABLE]

where $r_{1}\equiv k+\gamma\textrm{e}^{\varDelta\mu/2}$, $r_{2}\equiv k+\gamma\textrm{e}^{E/6}+\gamma\textrm{e}^{\varDelta\mu/2}$, $r_{3}\equiv k\textrm{e}^{E/3}+\gamma\textrm{e}^{E/6}+\gamma\textrm{e}^{\varDelta\mu/2}$, $r_{4}\equiv k\textrm{e}^{2E/3}+\gamma\textrm{e}^{E/6}+\gamma\textrm{e}^{\varDelta\mu/2}$, and $r_{5}\equiv k\textrm{e}^{E}+\gamma\textrm{e}^{E/6}$. The first seven states are related to the inactive form of the protein, with state $1$ corresponding to $C_{0}$ and state $7$ corresponding to $C_{6}$. The last seven states are related to the active form of the protein, with state $8$ corresponding to $\tilde{C}_{6}$ and state $14$ corresponding to $\tilde{C}_{0}$.

The phosphorylation level of the protein is a state function. Its average is given by the expression

[TABLE]

$P_{C_{i}}(t)$ is the probability of configuration $C_{i}$ at time $t$. Initially the phosphorylation level is zero with the protein in state $C_{0}$, i.e., $P_{C_{0}}(0)=1$. The probabilities at time $t$ can be calculated with the expression

[TABLE]

where $\mathbf{P}(0)=(1,0,0,0,0,0,0,0,0,0,0,0,0,0)^{T}$ and $\mathbf{L}$ given in Eq. (33). We have calculated the phosphorylation level of the KaiC protein as a function of time from Eq. (34) and Eq. (35). The result is shown in Fig. 10. Clearly the exponential decay of the amplitude and the period of oscillation are determined by the first eigenvalue of the matrix in Eq. (33) (see caption of Fig. 10).

Appendix E Relation between $\mathcal{R}$ and entropy production

In this appendix, we discuss the relation between the ratio $\mathcal{R}$ and the entropy production from stochastic thermodynamics Seifert (2012), which we denote by $\sigma$.

For a unicyclic network with uniform rates $k^{+}=e^{\mathcal{A}/N}k^{-}$ this entropy production is given by

[TABLE]

For large $N$, we obtain $\mathcal{R}=\mathcal{A}/2\pi$, as in Eq. (6), and, from $X_{I}=(\textrm{e}^{\mathcal{A}/N}-1)k_{-}\sin(2\pi/N)$, we obtain $\varDelta Q\equiv 2\pi\sigma/X_{I}=\mathcal{A}$. If $\sigma$ is interpreted as the rate of heat dissipated to the environment Seifert (2012), $\varDelta Q$ is the dissipated heat per period of oscillation. Hence, for a unicyclic network with large number of states $N$, we find

[TABLE]

This expression is a particular case of the relation found in Cao et al. (2015), which states that after some critical value $\Delta Q_{c}$, for which oscillations set in, the inverse of the number of biochemical oscillations decay to some plateau with $(\Delta Q-\Delta Q_{c})^{-1}$. For our particular case both $\Delta Q_{c}$ and the plateau are zero. In Cao et al. (2015) this relation was demonstrated to be fulfilled for several different models.

While this relation is true for an unicyclic network with uniform rates that maximize $\mathcal{R}$, for arbitrary rates $\mathcal{R}$ can also decrease with an increase in $\Delta Q$. The entropy production for the generic unicylic model in Eq. (1) reads

[TABLE]

where $P_{i}$ is the stationary probability of state $i$. As an example, we consider the unicyclic model for $N=3$ with $k_{1}^{-}=k_{2}^{-}=k_{3}^{-}=1$, $k_{1}^{+}=k_{3}^{+}=\textrm{e}^{\mathcal{A}/4}$ and $k_{2}^{+}=\textrm{e}^{\mathcal{A}/2}$. In Fig. 11 we show that $\mathcal{R}$ as a function of $\varDelta Q=2\pi\sigma/X_{I}$ has a maximum, where we vary the affinity $\mathcal{A}$. Therefore, the number of coherent oscillations can also decrease with an increase in energy dissipation. A similar behavior has been observed in Qian and Qian (2000), with the main difference that instead of varying $\varDelta\mu$ the authors vary the temperature. The maximum of $\mathcal{R}$ as a function of temperature was identified as stochastic resonance.

The same maximum was also observed for the single KaiC model, as shown in Fig. 11. In this case, the entropy production can be written as

[TABLE]

where $P_{C_{i}}$ is the stationary probability of state $C_{i}$. In this expression we used the Schnakenberg cycle decomposition of the entropy production Seifert (2012).

Appendix F Activator-inhibitor model

In this appendix, we define the activator-inhibitor model from Cao et al. (2015). The model has four different chemical species: the activator $R$, the inhibitor $X$, the enzyme $M$ and the phosphatase $K$. An enzyme $M$ can be in four different states, which form the phosphorylation cycle

[TABLE]

The concentrations of $ATP$, $ADP$ and $P_{i}$ are assumed to be fixed. From the generalized detailed balance relation, the rates in Eq. (40) fulfill

[TABLE]

where $\varDelta\mu$ is the free energy liberated in one $ATP$ hydrolysis. The activator $R$ catalyzes the phosphorylation of the enzyme $M$ and the phosphatase $K$ catalyzes the dephosphorylation of $M$.

The enzyme in the phosphorylated state $M_{p}$ catalyzes the creation of both the activator $R$ with rate $l_{0}$ and the inhibitor $X$ with rate $l_{3}$. The activator $R$ can also be spontaneously created with a rate $l_{1}$. The inhibitor $X$ catalyzes the degradation of $R$ with rate $l_{2}$ and can be spontaneously degraded with a rate $l_{4}$. Hence, we have the following chemical reactions

[TABLE]

These are equilibrium reactions and a cycle that involves only them must have zero affinity. The only way to get a cycle with nonzero affinity in this model is to use the chemical reactions in Eq. (40). The reversed rates $\epsilon$ are assumed to be very small so that we can set them to zero. Formally, they must be nonzero for thermodynamic consistency, however, in a numerical simulation we can just set them to zero, instead of using a very small $\epsilon$ that will lead to the same results.

The total number of enzymes $N_{M_{T}}=N_{M}+N_{M_{p}K}+N_{M_{p}}+N_{MR}$, and that of phosphatases $N_{K}$ are conserved. The number of activators $R$ fulfills $N_{R}\geq 1$ and the number of inhibitors $X$ fulfills $N_{X}\geq 1$. A state of the system is then determined by the vector $\mathbf{N}=(N_{R},N_{X},N_{M},N_{M_{p}K},N_{M_{p}},N_{MR})$. The volume of the system is written as $V$ and the concentration of the chemical species $X$, for example, is denoted by $n_{X}\equiv N_{X}/V$. The master equation that defines this model reads

[TABLE]

where $P(\mathbf{N})$ is the probability to be in state $\mathbf{N}$ at time $t$ and $\delta\equiv 1/V$. Note that the rates $l_{2}$, $k_{1}^{\pm}$, $k_{3}^{\pm}$ have dimension $V^{-1}t^{-1}$, whereas the other rates have dimension $t^{-1}$. In the above equation, we set to zero the probability of configurations that violate the constraints $N_{R}\geq 1$, $N_{X}\geq 1$, $N_{M_{T}}=N_{M}+N_{M_{p}K}+N_{M_{p}}+N_{MR}$ and $N_{K}\geq N_{M_{p}K}$.

We have performed continuous time Monte Carlo simulations of this model and calculated the correlation function

[TABLE]

where the brackets denote an average over stochastic trajectories and $\langle N_{X}\rangle$ is the average number of $X$ in the stationary state. The initial condition $N_{X}(0)$ was sampled from the stationary distribution: in a simulation we let the system reach the stationary state before the time $t=0$.

The oscillating correlation function is shown in Fig. 12. The results presented in Fig. 4 were obtained in the following way. We have adjusted an exponential function to the peaks of the oscillation as shown in Fig. 12. The exponent gives $T/\tau$, where $T$ is the period of oscillation and $\tau$ the decay time. The ratio $\mathcal{R}$ was then estimated as $\mathcal{R}=2\pi\tau/T$. The parameters of the model are similar to the parameters used in Cao et al. (2015). They were set to $V=50$, $N_{M_{T}}=10V$, $l_{0}=l_{2}=l_{3}=1$, $l_{1}=0.4$, $l_{4}=0.5$, $k_{2}^{+}=k_{4}^{+}=k_{2}^{-}=k_{4}^{-}=15$, $k_{1}^{+}=k_{3}^{+}=100$, and $k_{1}^{-}=k_{3}^{-}=100\textrm{e}^{-\varDelta\mu/2}$. For the number of phosphatases we took $N_{K}=30,40,50$ and for the driving force $\varDelta\mu=10,12,14,16,18,20$.

An important aspect of this model from Sec. IV.2 is that the competition for a small number of phophatase $K$ synchronizes the cycles of different enzymes $M$. If $N_{K}$ is too small we have no oscillations in $M_{p}$. If $N_{K}$ is too large, the lack of competition for the phosphatase $K$ hinders oscillations in $M_{p}$. This feature is demonstrated in Fig. 13, where we show two time series of the four different states of the enzyme for $N_{K}=30$ and $N_{K}=150$. For $N_{K}=30$ we see clear oscillations, with the number $N_{M_{p}K}$ oscillating roughly between [math] and $30$. Whenever $N_{M_{p}K}(t)=30$, there is no phosphatase left and several enzymes get stuck in state $M_{p}$, synchronizing the phosphorylation cycles of these enzymes. For $N_{K}=150$, the number $N_{M_{p}K}$ stays below $150$. Hence, there is always free phosphatase in the system, resulting in no synchrony between different enzymes.