Functioning and robustness of a bacterial circadian clock

Cyanobacteria are the simplest known cellular systems that regulate their biological activities in daily cycles. For the cyanobacterium Synechococcus elongatus, it has been shown by in vitro and in vivo experiments that the basic circadian timing process is based on rhythmic phosphorylation of KaiC hexamers. Despite the excellent experimental work, a full systems level understanding of the in vitro clock is still lacking. In this work, we provide a mathematical approach to scan different hypothetical mechanisms for the primary circadian oscillator, starting from experimentally established molecular properties of the clock proteins. Although optimised for highest performance, only one of the in silico-generated reaction networks was able to reproduce the experimentally found high amplitude and robustness against perturbations. In this reaction network, a negative feedback synchronises the phosphorylation level of the individual hexamers and has indeed been realised in S. elongatus by KaiA sequestration as confirmed by experiments.

C i . The subscript indicates the phosphorylation state.
The time evolution equation for the hexamers in the i-th phosphorylation state is given Equivalently we obtain for the complexes The transition rates between hexamers and KaiBC complexes is denoted by c − and c + (see Each link drawn in the reaction schemes (c.f. Fig. S1) gives rise to a new reaction constant to be optimised for. The costfunction, H, to be minimised is given by with α i the weight constants for the different contributions H i .
The first term arises from optimisation for maximal amplitude and is given by with the fraction of phosphorylated sites of hexamers and complexes given by The second term punishes solutions that have long transients or show chaotic behaviour where A p is the amplitude within the p-th period, N p the number of periods, A Np the value for the last amplitude and σ = 0.2 the allowed standard deviation of amplitudes.
The last term minimises differences of phase between the native system and a system with concerted 5-fold increase of protein concentrations where T 1 and T 5 are the average period lengths of the system under native concentrations and 5-fold elevated concentrations, respectively. Typical orders of magnitude after optimisation for the cost-functions are O(H 1 ) = 1, O(H 2 ) = 0.001, O(H 3 ) = 0.01. In our case we took α 1 = α 2 = α 3 = 1. This means that the maximisation of amplitude is the most important effect, followed by the minimisation of phase difference between the oscillations of the control system and the over-expressed system. The assurance that transient dynamics has declined is the weakest condition. The position of the links in the networks (see Table S1 and Table S2)   three reaction parameters is close to zero). The two-link topologies are shown in Table S2.
The reaction constants obtained from the optimisation procedure for maximal amplitude for the reaction networks Fig. 3A-C are listed in Table S3. Here, E is the amplitude of oscillation after transients have declined in percent of the maximal achievable amplitude.

C. Parameter Optimisation for potential in vitro circadian oscillators
The potential oscillatory candidates found by the reaction network scan have to be re-translated in the biochemical reality. In particular: (i) KaiA and KaiB have to be included in the reaction network. (ii) The rates for phosphorylation and dephosphorylation in absence of KaiB have to lie within the measured error bars. (iii) It has to be accounted for the fact that complexes can exchange monomers.
In this section we constrain our model to the phosphorylation and dephosphorylation data published in (Kitayama et al., 2003) and (Xu et al., 2003), respectively. From a least-meansquare fit we obtained k + i = 3.9 · 10 −2 min −1 , i ∈ {0, 1, 2, 3, 4, 5}, for the phosphorylation rate if we assume equal reaction constants for all phosphorylation steps (Fig. S2A, solid line). Equivalently one could argue that the phosphorylation rates scale with the number of unphosphorylated residues, k + i = k + (N − i). The latter case gives k + = 1.0 · 10 −2 min −1 after least-square fitting and is shown by Fig. S2A, dashed line. From the data fits it cannot be decided which of two assumptions is the more realistic one but as we eventually optimise each of the reaction rates individually we stick to assumption that all reaction constants are equal for the rewiring approach. By a similar data fit we find k − i = 9.7 · 10 −3 min −1 for the dephosphorylation rate, assuming that 100% of the hexamers were fully phosphorylated at time t = 0 (Fig. S2B, solid line). This may not represent the physiological initial condition since the maximal achievable phosphorylation level via KaiA is close to 90%. A better fit is obtained by starting from the initial condition of 80% phosphorylated KaiC, where we found k − i = 7.8 · 10 −3 min −1 for the dephosphorylation rate (Fig. S2B, dashed line). Introducing these quantities in our model leads to oscillations with a period varying between 16 h and 30 h, depending on the setting of the parameters in the optimisation routine. This means that the principal time-scale of the Clock is fixed by the phosphorylation and dephosphorylation of the hexamers. Other processes (e.g. complex formation) can be expected to occur on a much more faster time-scale.
Since all parameters appear in linear form in the model equations, we can always rescale the time scale, t, of the system following t ′ = αt and all rates, k i , take then the form k ′ i = k i /α. The relevant information from the data fit is therefore that the rate k + i /k − i should be set between 3.9 · 10 −2 /9.7 · 10 −3 ≃ 4 and 3.9 · 10 −2 /7.8 · 10 −3 ≃ 5 in the model.
For the simulations in this work we have set k + i /k − i = 4. Leaving one parameter (the phosphorylation rate k + i ) free allows us to optimise the behaviour of the system without being constrained to a specific time-scale of oscillations.
with k + i ,k − i , k + i and k − i the kinetic rates for phosphorylation and dephosphorylation, respectively. The previous equations are valid for i = 1 . . . 5 and the boundary terms (i = 0, 6) of these equations read The complexes that appear in the experiments can be expressed in terms of our state variables as follows: were the stable complexes consist of one hexamer and six monomers of KaiB.
Assuming the unique stoichiometry for the KaiBC complexes of one hexamer and six KaiB dimers, the conservation equations for species KaiA and KaiB dimers are given by and lead to the concentrations of free KaiA and KaiB dimers (26) The model assumes that KaiA binding and dissociation is much faster than the phosphorylation steps, which is reasonable as the phosphorylation occurs on the basis of hours. Also the model assumes that the fraction of free hexamers is significantly larger than the fraction of hexamers existing in unstable complex form with KaiA dimers and stable complex form with 6 KaiB monomers. The model can only show significant oscillating behaviour if c + < c − .

E. Shuffling Process: Synchronisation by Monomers Exchange between Hexamers
We introduce the shuffling process through the master equation where P i (t) = C i (t)/ i C i (t) is the probability that a complex C i is in the i − th phosphorylation state and W(j → i) is the transition matrix, whose elements are the rates for a complex to change its phosphorylation level from j to i phosphorylated residues. The sum rule j W j→i = 1 implies i P i (t) = 1.
In our model we assume that two complexes C i and C j can exchange only one monomer when they get into contact. The rate equation for KaiBC phosphorylation then reads with boundary terms Here, i/N is the probability of a complex to loose one of its i phosphorylated monomers and (N − j)/N is the probability that the exchange partner donates one of its N − j nonphosphorylated monomers. The shuffling rate is here represented by the parameter k s i . For the optimisation for highest performance we set for simplicity all shuffling constants equal.
In Fig. S3 we suppressed the negative feedback loop due to KaiA sequestration ( K A i = ∞) and looked at the effect of shuffling alone as a synchronisation mechanism. We obtained through optimisation an amplitude of oscillations around a few percent. The corresponding reaction constants from the optimisation for highest amplitude and robustness including monomer shuffling (Fig. 4B)

II. A SIMPLIFIED ANALYTICAL MODEL OF THE CORE CLOCK
Assuming constant (KaiB independent) transition rates between hexamers and complexes at the highest and lowest phosphorylation state and equal phosphorylation rates we arrive at the following simplified equations: with k + (t) = k A K A A(t) and A(t) given by For the complexes we obtain In the following we assume that the transition rate from the lower to upper reaction chain is small, e.g. c − , c + ≪ k + (t), k − . The solution of Eq. (31) can then be approximated by its corresponding stationary state solution with N (t) a normalisation function determined from conservation of the total amount of hexamers, H T = N n=0 H n (t) + N n=0 C n (t), The time evolution of the complexes follows a Poisson process for 1 ≤ n ≤ N , whose explicit solution can be calculated by Laplace Transform of Eqs. (33) with G n (t − t ′ ) given by The solutions Eqs. (38 and 39) assume the initial conditions After optimisation the kinetic constants of the full problem Eqs. (12) and (13) for maximum amplitude of sustained oscillations, it turns out that the transition rates between the upper and lower reaction chains can be forced to follow c + ≪ c − without significant loss in amplitude. As a consequence the amount of hexamers (but not necessarily the contribution to the total phosphorylation level) is significantly smaller in the lower reaction chain than in the upper reaction chain. We use this fact to simplify the equations above and allow N i=0 C i (t) to be very small but with a very high binding affinity to KaiA. By allowing only the states n = 0 in the lower reaction chain to bind KaiA dimers the Eqs. (35-39) can be simplified to Here the time evolution for C n (t) is now represented by an integral equation, Eq. (42).
The three equations above can now be attributed the following functions: in a higher value for C n (t) that in turn decreases k + (t) ∼ 1/C n (t) that in turn lowers the level of H N (t) The items 1.-3. are essential ingredients for a negative feedback oscillator.