Translational control of gene expression via interacting feedback loops

TL;DR

This paper presents a mathematical model showing how feedback loops regulate translation, revealing complex dynamics like oscillations and bistability, which expand understanding of cellular control mechanisms.

Contribution

It introduces a novel multi-component framework demonstrating the role of feedback interactions in translation regulation.

Findings

Feedback loops induce oscillations in protein synthesis

Bistability emerges from feedback interactions

Complex dynamics suggest richer cellular control mechanisms

Abstract

Translation is a key step in the synthesis of proteins. Accordingly, cells have evolved an intricate array of control mechanisms to regulate this process. By constructing a multi-component mathematical framework for translation we uncover how translation may be controlled via interacting feedback loops. Our results reveal that this interplay gives rise to a remarkable range of protein synthesis dynamics, including oscillations, step-change and bistability. This suggests that cells may have recourse to a much richer set of control mechanisms than was previously understood.

Equations

$$\underset{n}{\underbrace{\text{P}+\text{P}+\dots \text{P}}}\text{autorightleftharpoons}{r}_1{r}_{2}n\text{P},\text{P}\text{autorightarrow}r\varnothing$$

$$\frac{d[n\text{P}]}{dt}=r_1{\left([\text{P}]\right)}^n-r_2[n\text{P}]\approx 0⟹[n\text{P}]=\frac{r_1}{r_2}{\left([\text{P}]\right)}^n=\frac{r_1}{r_2}{\left(\frac{N}{V}\right)}^n,$$

$$n\text{P}+\text{SS}\text{autorightleftharpoons}{k}_1{k}_{2}n\text{P}\cdot \text{SS}.$$

$$\frac{d[n\text{P}\cdot \text{SS}]}{dt}$$

$$i.e.$$

$$\frac{[\text{SS}]}{[\text{SS}]+[n\text{P}\cdot \text{SS}]}=\frac{[\text{SS}]}{[\text{SS}]+\frac{k_1}{k_2}[n\text{P}][\text{SS}]}=\frac{k_2}{k_2+k_1[n\text{P}]},$$

$$f(N):=\frac{1}{1+{\left(N/\theta \right)}^n},$$

$${\rho }_1=1-\frac{1}{4{\alpha }_{eff}},{\rho }_L=\frac{1}{4{\beta }_{eff}},\rho =\frac{1}{2},J=\frac{1}{4}.$$

$${\rho }_1={\alpha }_{eff},{\rho }_L=\frac{{\alpha }_{eff}(1-{\alpha }_{eff})}{{\beta }_{eff}},\rho ={\alpha }_{eff},J={\alpha }_{eff}(1-{\alpha }_{eff}).$$

$${\rho }_1=1-\frac{{\beta }_{eff}(1-{\beta }_{eff})}{{\alpha }_{eff}},{\rho }_L=1-{\beta }_{eff},\rho =1-{\beta }_{eff},J={\beta }_{eff}(1-{\beta }_{eff}).$$

$$\frac{\alpha }{2\alpha +k\left(1+I^n\right)}⩽\{\begin{array}{rlrl}& \beta & ⩽\frac{\alpha k}{1+I^n-2\alpha }& \text{if}\alpha <\frac{1}{2}\left(1+I^n\right),\\ & \beta & & \text{otherwise},\end{array}$$

$${\alpha }_{eff}=\frac{\alpha }{1+I^n}+\frac{k}{2\left(\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+I^n\right)\right)}\right)}\text{mbox}and{\beta }_{eff}=\frac{1}{2}\left(\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+I^n\right)\right)}\right).$$

$$(\alpha ,\beta )=\left(\frac{1-k}{2}\left(1+I^n\right),\frac{1-k}{2}\right).$$

$$\alpha <\{\begin{array}{rlrl}& \beta \left(1+{\left(\frac{4I(\beta +k)(1-\beta )}{(1+k{)}^2}\right)}^n\right)& & \text{if}2\beta +k<1,\\ & \frac{\beta }{2\beta +k}\left(1+I^n\right)& & \text{otherwise}.\end{array}$$

$$\alpha >\{\begin{array}{rlrl}& \beta \left(1+{\left(\frac{4I(\beta +k)(1-\beta )}{(1+k{)}^2}\right)}^n\right)& & \text{if}2\beta +k<1,\\ & \frac{k\beta }{1-2\beta }\left(1+I^n\right)& & \text{otherwise},\end{array}\beta <\frac{1}{2}.$$

$${\beta }_{eff}$$

$$4^{n+1}\alpha (\theta r{)}^n{\beta }_{eff}^2-4^{n+1}\alpha \beta (\theta r{)}^n{\beta }_{eff}-k\beta \left((4\theta r{)}^n+1\right)=0.$$

$${\beta }_{eff}=$$

$${\beta }_{eff}=\frac{1}{2}\left(\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)\right)}\right).$$

$${\alpha }_{eff}$$

$${\beta }_{eff}=\frac{\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)\right)}}{2}⩾\frac{1}{2},$$

$$\beta ⩾\frac{\alpha }{2\alpha +k\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)},$$

$${\alpha }_{eff}=\frac{\alpha }{1+{\left(\frac{1}{4\theta r}\right)}^n}+\frac{k}{2\left(\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)\right)}\right)}⩾\frac{1}{2},$$

$$\frac{k}{2\left(\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)\right)}\right)}⩾\frac{1}{2}-\frac{\alpha }{1+{\left(\frac{1}{4\theta r}\right)}^n}.$$

$$k⩾2\left(\frac{1}{2}-\frac{\alpha }{1+{\left(\frac{1}{4\theta r}\right)}^n}\right)\left(\beta +\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)\right)}\right),$$

$$k-\left(1-\frac{2\alpha }{1+{\left(\frac{1}{4\theta r}\right)}^n}\right)\beta ⩾\left(1-\frac{2\alpha }{1+{\left(\frac{1}{4\theta r}\right)}^n}\right)\sqrt{\beta \left(\beta +\frac{k}{\alpha }\left(1+{\left(\frac{1}{4\theta r}\right)}^n\right)\right)}>0.$$

$$\beta ⩽\frac{\alpha k}{1+{\left(\frac{1}{4\theta r}\right)}^n-2\alpha }.$$

$${\alpha }_{eff}=\frac{\alpha }{1+{\left(\frac{{\alpha }_{eff}(1-{\alpha }_{eff})}{\theta r}\right)}^n}+k\frac{{\alpha }_{eff}(1-{\alpha }_{eff})}{\beta +k(1-{\alpha }_{eff})}.$$

$$P({\alpha }_{eff}):=\beta {\alpha }_{eff}^{n+1}(1-{\alpha }_{eff}{)}^n+(\theta r{)}^n(\beta +\alpha k){\alpha }_{eff}-\alpha (\theta r{)}^n(\beta +k)=0.$$

$${\alpha }_{eff}<{\beta }_{eff}=\beta +k(1-{\alpha }_{eff}),⟹(1+k){\alpha }_{eff}<\beta +k,$$

$${\alpha }_{eff}<\frac{\beta +k}{1+k}.$$

Full text

Translational control of gene expression via interacting feedback loops. Supplemental Material.

Liang Wang

University of Dundee, Dundee, United Kingdom

  

M. Carmen Romano

University of Aberdeen, Aberdeen, United Kingdom

  

Fordyce A. Davidson*

University of Dundee, Dundee, United Kingdom

Corresponding Author

1 Protein Dynamics and Effective Vacancy of the Start Site

Protein complex formation and loss, respectively, can be represented as follows:

[TABLE]

where P is protein, $n\geq 1$ proteins form a protein complex $n\text{P}$ and $r_{1},\,r_{2}$ and $r$ are rate constants, where $r_{1}$ has dimensions of $(\text{volume})^{n-1}/\text{time}$, and $r_{2}$ and $r$ have dimensions of 1/time. Applying a quasi-steady-state approximation for the complex concentration:

[TABLE]

where $[\,\cdot\,]$ represents concentration, $N$ denotes the number of molecules of P, and $V$ denotes the cell volume. The protein complex is assumed to bind and unbind from some location at the start site (SS) at rates $k_{1},k_{2}$, respectively:

[TABLE]

Writing (S2) as a differential equation and applying the quasi-steady state approximation as above yields

[TABLE]

This quasi-steady state approximation leads to an effective vacancy of the SS from the perspective of ribosome binding events. The probability that the SS is vacant is equivalent to

[TABLE]

where we have used (S3). This probability can be written in terms the number of proteins, $N$, as follows. Using (S1) and (S4) it follows that the probability that the SS is free is

[TABLE]

where $\theta=V\sqrt[n]{r_{2}k_{2}/r_{1}k_{1}}$.

2 Phases and Phase Boundaries

Maximal Current Phase: $\alpha_{eff},\,\beta_{eff}\geqslant 1/2$

[TABLE]

Low Density Phase: $\alpha_{eff}<\beta_{eff},\,\alpha_{eff}<1/2$

[TABLE]

High Density Phase: $\alpha_{eff}>\beta_{eff},\,\alpha_{eff}>1/2$

[TABLE]

In each case the corresponding steady state protein level is given by $N^{*}:=J/r$. We now define the phase boundaries in a series of lemmas. Proofs are lengthy and hence grouped in the next section for easier exposition.

Lemma 1

The Maximal Current Phase is defined by

[TABLE]

*where $I:=\frac{1}{4r\theta}$. Within this region, *

[TABLE]

Remark 1

Note that the stacking of the left and right inequalities (S9) in case $\alpha<\dfrac{1}{2}\left(1+I^{n}\right)$ is necessary and therefore generates a lower bound for $\alpha$ and $\beta$ marked by the intersection of the two functions of $\alpha$. It is straightforward to compute this intersection:

[TABLE]

Lemma 2

The Low Density Phase is defined by

[TABLE]

Within this region, there exists a unique, positive expression for $\alpha_{eff}$ that yields unique, positive expressions for $\beta_{eff}$ and $N^{*}$.

Lemma 3

The High Density Phase is defined by

[TABLE]

Within this region, there exist either one or three positive solutions, $\beta_{eff}$, that yield positive expressions for $\alpha_{eff}$ and $N^{*}$.

3 Proofs

Proof of Lemma 1 Substituting the effective rates into (S6) yields:

[TABLE]

Therefore, we obtain a quadratic equation for $\beta_{eff}$:

[TABLE]

Subsequently,

[TABLE]

Since we require $\beta_{eff}$ to be positive and $\sqrt{\beta^{2}+\dfrac{k\beta\left((4\theta r)^{n}+1\right)}{4^{n}\alpha(\theta r)^{n}}}>\beta$, it follows that:

[TABLE]

Then $\alpha_{eff}$ is easy to obtain:

[TABLE]

The maximal current phase is determined by the conditions $\alpha_{eff},\,\beta_{eff}\geqslant\dfrac{1}{2}$ i.e.:

[TABLE]

so $\sqrt{\beta\left(\beta+\dfrac{k}{\alpha}\left(1+\left(\dfrac{1}{4\theta r}\right)^{n}\right)\right)}\geqslant 1-\beta$. Since $1-\beta>0$ for $\beta<1$, we have

[TABLE]

and

[TABLE]

namely,

[TABLE]

If $\dfrac{1}{2}-\dfrac{\alpha}{1+\left(\dfrac{1}{4\theta r}\right)^{n}}\leqslant 0$, i.e., $\alpha\geqslant\dfrac{1}{2}\left(1+\left(\dfrac{1}{4\theta r}\right)^{n}\right)$, then (S17) is always true.

If $\alpha<\dfrac{1}{2}\left(1+\left(\dfrac{1}{4\theta r}\right)^{n}\right)$, we have:

[TABLE]

subsequently,

[TABLE]

From (S18), we have a weak condition

[TABLE]

The result follows by grouping the inequalities determined above and on writing the final expression in terms of the intensity factor.

$\Box$

Proof of Lemma 2 Substituting the effective rates into (S7) yields:

[TABLE]

Again after some manipulation we obtain:

[TABLE]

If there exist positive solutions, $\alpha_{eff}$, to (S20), then from (S7), $\beta_{eff}=\beta+k\left(1-\alpha_{eff}\right)$ and $N^{*}=\alpha_{eff}(1-\alpha_{eff})/r$. Next we consider the conditions for the existence of $\alpha_{eff}$ and thus define the phase boundary for the LD phase. LD phase is by defined by $\alpha_{eff}<\beta_{eff},\,\alpha_{eff}<\dfrac{1}{2}$. From $\alpha_{eff}<\beta_{eff}$, we have:

[TABLE]

therefore

[TABLE]

On the other hand, we require $\alpha_{eff}<\dfrac{1}{2}$. Therefore $\alpha_{eff}<\min\left(\dfrac{1}{2},\,\dfrac{\beta+k}{1+k}\right)$.

If $\dfrac{\beta+k}{1+k}\geqslant\dfrac{1}{2}$, i.e., $2\beta+2k\geqslant 1+k$, $2\beta+k\geqslant 1$, then $\alpha_{eff}<\dfrac{1}{2}$. The derivative of $P$ with respect to $\alpha_{eff}$ is

[TABLE]

Since both $1-\alpha_{eff}$ and $1-2\alpha_{eff}$ are positive in $\left(0,\,\dfrac{1}{2}\right)$, so $P^{\prime}>0$ in $\left(0,\,\dfrac{1}{2}\right)$. Moreover, $P(0)=-\alpha(\theta r)^{n}(\beta+k)<0$. Therefore, a real, positive solution of (S20) exists if and only if $f\left(\dfrac{1}{2}\right)>0$ (in this case, the solution is unique). Now,

[TABLE]

so $P\left(\dfrac{1}{2}\right)>0$ requires $\beta\left(1+\left(\dfrac{1}{4\theta r}\right)^{n}\right)-\alpha(2\beta+k)>0$, namely,

[TABLE]

If $\dfrac{\beta+k}{1+k}<\dfrac{1}{2}$, i.e., $2\beta+k<1$, then $\alpha_{eff}<\dfrac{\beta+k}{1+k}$. Since $P$ is increasing in $\left(0,\,\dfrac{1}{2}\right)$ and $P(0)<0$, so a solution of (S20) exists if and only if $P\left(\dfrac{\beta+k}{1+k}\right)>0$ (and again the solution is unique).

[TABLE]

provided

[TABLE]

The proof is complete by writing the inequalities in terms of the Intensity Factor.

$\Box$

Proof of Lemma 3 Substituting the effective rates into (S8) yields:

[TABLE]

Further algebra yields the equation

[TABLE]

With $\beta_{eff}$ a solution of (S24), then from (S8), $\rho=1-\beta_{eff}$, $J=\beta_{eff}(1-\beta_{eff})$ and $N^{*}=\beta_{eff}(1-\beta_{eff})/r$. Hence, the effective initiation rate is given by

[TABLE]

One of the HD conditions is $\beta_{eff}<\dfrac{1}{2}$, i.e. $\dfrac{1}{2}>\beta_{eff}=\beta+k(1-\rho_{1})>\beta$, namely, $\beta<\dfrac{1}{2}$. Moreover, from the second condition for HD $\alpha_{eff}>\beta_{eff}$, we have

[TABLE]

Therefore, we require $\beta_{eff}<\min\left\{\dfrac{1}{2},\,\dfrac{\beta+k}{1+k}\right\}$.

If $\dfrac{\beta+k}{1+k}\geqslant\dfrac{1}{2}$, i.e., $2\beta+k\geqslant 1$ then we need to check the existence of solutions to (S24) in $\left(0,\,\dfrac{1}{2}\right)$. To this end, note that

[TABLE]

Thus, $Q^{\prime\prime}>0$ in $\left(0,\,\beta_{eff}^{-}\right)$ and $Q^{\prime\prime}<0$ in $\left[\beta_{eff}^{-},\,\dfrac{1}{2}\right)$ with the maximum of $Q^{\prime}$ obtained at $\beta_{eff}^{-}$. Moreover,

[TABLE]

If $Q^{\prime}\left(\beta_{eff}^{-}\right)<0$, then $Q^{\prime}<0$ ($Q$ is decreasing) in $\left(0,\,\dfrac{1}{2}\right)$. Besides, $Q(0)=\beta(\theta r)^{n}(\alpha+k)>0$. Thus, the solution of (S24) $Q=0$ exists if and only if $Q\left(\dfrac{1}{2}\right)<0$, and the solution is unique.

If $Q^{\prime}\left(\beta_{eff}^{-}\right)>0$, then there are two solutions to $Q^{\prime}=0$. Let $\beta_{eff}^{1}$ and $\beta_{eff}^{2}$ denote the solutions and $0<\beta_{eff}^{1}<\beta_{eff}^{2}<\dfrac{1}{2}$. Thus, $Q^{\prime}<0$ ($Q$ is decreasing) in $\left(0,\,\beta_{eff}^{1}\right)\cup\left(\beta_{eff}^{2},\,\dfrac{1}{2}\right)$, $Q^{\prime}>0$ ($Q$ is increasing) in $\left(\beta_{eff}^{1},\,\beta_{eff}^{2}\right)$. Moreover, $Q(0)=\beta(\theta r)^{n}(\alpha+k)>0$, and $Q\left(\beta_{eff}^{1}\right)$ is the local minimum. We have following cases.

(I) If $Q\left(\beta_{eff}^{1}\right)<0$ and $Q\left(\beta_{eff}^{2}\right)<0$, then $Q<0$ in $\left(\beta_{eff}^{1},\,\dfrac{1}{2}\right)$, there exists a unique solution to (S24) in $\left(0,\,\beta_{eff}^{1}\right)$. In this case, $Q\left(\dfrac{1}{2}\right)<0$.

(III) If $Q\left(\beta_{eff}^{1}\right)<0$ and $Q\left(\beta_{eff}^{2}\right)>0$, then (S24) has three solutions if $Q\left(\dfrac{1}{2}\right)<0$ or two solutions if $Q\left(\dfrac{1}{2}\right)>0$. If $Q\left(\dfrac{1}{2}\right)>0$, we have:

[TABLE]

which gives:

[TABLE]

However, from (S9), (S25) sets the system in the MC phase and hence is not a feasible condition. Therefore, in this case, $Q\left(\dfrac{1}{2}\right)<0$ and there are three solutions to (S24).

Overall, for $2\beta+k\geqslant 1$, the existence of the solution of $\beta_{eff}$ requires $Q\left(\dfrac{1}{2}\right)<0$:

[TABLE]

which gives:

[TABLE]

If $\dfrac{\beta+k}{1+k}<\dfrac{1}{2}$, i.e., $2\beta+k<1$, then we require $\beta_{eff}<\dfrac{\beta+k}{1+k}$. Similar to the case for $2\beta+k\geqslant 1$, the existence of the solution of $\beta_{eff}$ requires $Q\left(\dfrac{\beta+k}{1+k}\right)<0$ and one or three solutions can be obtained. Note,

[TABLE]

provided

[TABLE]

The result follows directly as above.

$\Box$

4 Simulation

In every Monte-Carlo step (MCS), $L+2$ random numbers were chosen, corresponding to all possible reactions. Protein degradation was modelled using $n_{r}$ (the current number of proteins) attempts, each with corresponding degradation rate, $r$. The initiation and termination probabilities were governed by $\alpha_{eff}$ and $\beta_{eff}$. The processes of complex formation and subsequent interaction with the Start Site were implemented using the quasi-steady state approximations detailed above. [Test cases where all processes were explicitly accounted, such as protein complex formation and protein complex binding, etc., displayed little qualitative or quantitative difference to the hybrid Monte-Carlo used here, but required significantly greater computational time.]

Time average was calculated from a simulation consisting of $10^{8}$ Monte-Carlo steps (MCS), where the first $10^{6}$ MCS were disregarded to ensure steady-state. Then, the integration time was divided in windows of $10^{4}$ MCS, so that the standard deviation could be calculated for each value of $\alpha$.

Typical simulations are shown in Fig.S1A. The lag before the first protein is made is approximately 500 time units corresponding to the length of the lattice. After a transient rise in protein number, the model predicts steady (average) protein number. The average long term behaviour of the Monte-Carlo simulations matched closely the corresponding analytical steady states $N^{*}$ derived above.

In Fig. S1B, the phase diagram for the average particle density is shown. Three separate phases are clearly identified by relatively sharp changes in average density. This phase separation is well-approximated by the analytical phase boundaries for the modified TASEP derived above. Within the phases the qualitative match of the phase density with that predicted by the steady state theory is generally in good agreement - this is easiest identified in the maximal current phase where the steady state theory predicts the value $\rho=1/2$. Agreement is generally good in the other phases.

5 Dynamics in the Modified TASEP

5.1 Onset of oscillations is characterised by a Hopf Bifurcation

The mean time for a particle to transit the lattice can be estimated as

[TABLE]

in the LD regime. After substitution and some rearranging, the effective ribosome binding rate in LD is

[TABLE]

Appealing to the mean field approximation and setting $\langle N\rangle=N^{*}$, the mean transit time is then defined as

[TABLE]

The delay differential equation for protein copy number in the LD regime is:

[TABLE]

where we have written $T$ for $\langle T\rangle$ for ease of exposition here and below. Steady states, $N^{\ast}$, of (S30) are exactly the solutions generated by substituting solutions of (S20) into the definition $N^{*}=\alpha_{eff}(1-\alpha_{eff})/r$.

To check the stability of $N^{\ast}$, we consider perturbations, $\delta(t)$, with $|\delta|\ll 1$. Substituting $N(t)=N^{\ast}+\delta(t)$ into (S30), using Taylor expansion and dropping the higher order term yields

[TABLE]

where

[TABLE]

as $\alpha_{eff}<1/2$ in the LD phase by definition. Setting $\delta(t)=Ae^{\lambda t}$ and on substitution into (S31), we get (hereon setting $B(N^{*})=B$ for ease of notation)

[TABLE]

Setting $\lambda=\mu+i\omega$, standard arguments yield a critical wave number at which stability is lost, namely, $\omega^{2}=B^{2}-r^{2}$ and hence a necessary condition is $B>r$. Some algebra eliminates $\omega$ to defines the Hopf locus in terms of system parameters thus:

[TABLE]

Hence, the existence of periodic solutions to (S30) is guaranteed, at least for values of $(\alpha,\beta)$ sufficiently close to the locus. The Hopf locus in the $\alpha-\beta$ plane is shown in Figure SS4 along with the effects of varying $I$ and $k$ (Figs. S3 A and B). Finally, note that the necessary condition $B-r>0$, yields:

[TABLE]

where

[TABLE]

Hence, $B-r>0\implies g\left(\left(\dfrac{N^{\ast}}{\theta}\right)^{n}\right)>0$. Notice that if $n=1$, $g\left(\left(\dfrac{N^{\ast}}{\theta}\right)^{n}\right)=\left(-(\alpha+1)\beta-2\alpha k\right)\dfrac{N^{\ast}}{\theta}-(\beta+\alpha k)(1-\alpha)<0$. Hence the Hopf bifurcation cannot occur in this case. Therefore, co-operativity is a necessary condition for Hopf bifurcation to occur.

For $n>1$, the discriminant of $g\left(\left(\dfrac{N^{\ast}}{\theta}\right)^{n}\right)$ is given by:

[TABLE]

and $\Delta>0$ certainly when $\alpha<1$. Therefore, in this case $g\left(\left(\dfrac{N^{\ast}}{\theta}\right)^{n}\right)=0$ has two distinct real roots:

[TABLE]

where $\Delta$ is given in (S36) and $\left(\dfrac{N^{\ast}}{\theta}\right)^{n}_{-}$ is negative. Hence, it is easy to conclude that $g\left(\left(\dfrac{N^{\ast}}{\theta}\right)^{n}\right)>0$ requires:

[TABLE]

The inequality (S37) and setting the right hand side of (S30) to zero, together imply

[TABLE]

Rewriting this in terms of the Intensity Factor $I$ yields the condition $I>F(\bf{\Lambda})$ where $F(\bf{\Lambda})$ is a multiple of the reciprocal of the right hand side of (S38).

5.2 Numerical Simulations of the Characteristic DDE

The DDE (S30) was solved numerically using the code Matlab dde. Typical behaviour for a range of values of $\alpha$ is shown in the Figures below.

The amplitude and period of the oscillations generated by (S30) are in good agreement with those generated by Monte-Carlo simulations (cf. Fig. S3D and Fig.2A in the main text).

5.3 Periodic oscillations are restricted to the LD Phase

In the MC phase, $J\equiv 1/4$ and $\alpha_{eff}$ is given by (S15) and thus the loading rate is independent of time. Therefore $N\equiv 1/4r$ is the unique asymptotic solution related to any non-zero initial data. Similarly, in the HD phase, $J=\beta_{eff}(1-\beta_{eff})$, with $\beta_{eff}$ given by the solutions of (S24). In this case, $\beta_{eff}=\beta_{eff}(N(t))$ and thus $J=J(N(t))$ and the characteristic equation $dN/dt=J-rN$ is a scalar, ordinary differential equation with constant coefficients. Standard theory precludes periodic solutions. The behaviour of the characteristic equations therefore supports the numerical simulations to suggest that (almost) periodic solutions are only realised in the LD phase.

We computed an autocorrelation function based measure $ACF(\tau):=E[(N(t)-\mu)(N(t+\tau)-\mu)]/\sigma^{2}$, where $\mu$ is the average and $\sigma$ the standard deviation of the time series $N(t)$ of the number of proteins. $N(t)$ means the time series at time $t$, and $N(t+\tau)$ the time series at time $t+\tau$. Results show that the Hopf locus associated with the DDE provides a reasonable approximation for the parameter region in the $\alpha-\beta$ plane for which oscillations are observed (Fig. S4).

6 Bistability