Robust circadian clocks from coupled protein modification and transcription-translation cycles
aa r X i v : . [ q - b i o . M N ] A p r Robust circadian clocks from coupled protein modification andtranscription-translation cycles
David Zwicker, David K. Lubensky, and Pieter Rein ten Wolde FOM Institute for Atomic and Molecular Physics (AMOLF),Science Park 104, 1098 XG Amsterdam, The Netherlands Department of Physics, University of Michigan, Ann Arbor, MI 48109-1040
The cyanobacterium
Synechococcus elongatus uses both a protein phosphorylation cycle and atranscription-translation cycle to generate circadian rhythms that are highly robust against bio-chemical noise. We use stochastic simulations to analyze how these cycles interact to generatestable rhythms in growing, dividing cells. We find that a protein phosphorylation cycle by itself isrobust when protein turnover is low. For high decay or dilution rates (and co mpensating synthe-sis rate), however, the phosphorylation-based oscillator loses its integrity. Circadian rhythms thuscannot be generated with a phosphorylation cycle alone when the growth rate, and consequentlythe rate of protein dilution, is high enough; in practice, a purely post-translational clock ceases tofunction well when the cell doubling time drops below the 24 hour clock period. At higher growthrates, a transcription-translation cycle becomes essential for generating robust circadian rhythms.Interestingly, while a transcription-translation cycle is necessary to sustain a phosphorylation cy-cle at high growth rates, a phosphorylation cycle can dramatically enhance the robustness of atranscription-translation cycle at lower protein decay or dilution rates. Our analysis thus predictsthat both cycles are required to generate robust circadian rhythms over the full range of growthconditions.
I. INTRODUCTION
Many organisms use circadian clocks to anticipatechanges between day and night [1]. It had longbeen believed that these clocks are driven primarily bytranscription-translation cycles (TTCs) built on nega-tive feedback. However, while some circadian clocks canmaintain robust rhythms for years in the absence of anydaily cue [1], recent experiments have vividly demon-strated that gene expression is often highly stochastic[2]. This raises the question of how these clocks can berobust against biochemical noise. In multicellular organ-isms, the robustness might be explained by intercellularinteractions [3, 4], but it is now known that even unicel-lular organisms can have very stable circadian rhythms.The clock of the cyanbacterium
Synechococcus elongatus ,for example, has a correlation time of several months [9],even though the clocks of the different cells in a popu-lation hardly interact with one another [6, 9]. Interest-ingly, it has recently been discovered that the
S. elon-gatus clock also includes a protein phosphorylation cycle(PPC) that can run independently of the TTC [7, 8]. Ithas been suggested that this protein modification oscil-lator [8], possibly in combination with the transcription-translation oscillator [15], is responsible for the clock’sstriking noise resistance. Here, we use mathematicalmodeling to study how these two clocks interact in grow-ing, dividing cells. We find that the PPC alone is verystable when cell growth is negligible, but it begins to failas the growth rate increases, and for doubling times typ-ical of
S. elongatus it cannot function without help froma TTC. We then use stochastic simulations to show thata PPC can, however, dramatically enhance the robust-ness of a TTC, especially when cells are growing slowly.The two mechanisms thus perform best in complemen- tary situations, and it is likely that both oscillators arenecessary to maintain even minimal clock stability acrossthe full range of conditions and growth rates that bacte-ria encounter in the wild. Importantly, the TTC and thePPC in
S. elongatus are much more tightly intertwinedthan conventional coupled phase oscillators; as a result,the combination of the two far outperforms not just eachof its two components individually, but also a hypotheti-cal system in which the two parts are coupled in normaltextbook fashion [10].In
S. elongatus , the central components of the clockare the three genes kaiA , kaiB , and kaiC [11]. Undercontinuous light conditions, the levels of mRNA from the kaiBC operon and of the protein KaiC oscillate in a cir-cadian fashion with a 6 hour lag between their peaks [12].Overexpression of phosphorylated KaiC abolishes kaiBC expression [13, 14], while a transient increase in KaiC re-sets the phase of the oscillator [12]. These observationsled to the proposal that the Kai system is a transcription-translation oscillator, with KaiC negatively regulating itsown transcription. In 2005, however, Kondo and cowork-ers showed that KaiC is phosphorylated in a cyclical man-ner with a period of 24 hours, even when kaiBC transcrip-tion is inhibited [7]. Still more remarkably, the rhythmicphosphorylation of KaiC could be reconstituted in thetest tube in the presence of only KaiA, KaiB, and ATP[8]. This raised the possibility that the principal pace-maker of the clock is not a transcription-translation cycle,but a protein phosphorylation cycle [8]. Yet, in 2008, thesame group showed that circadian oscillations of gene ex-pression persist even when KaiC is always held in a highlyphosphorylated state [15]. They thus concluded that theclock is driven by both a TTC and a PPC, and suggestedthat the interactions between the two oscillators may en-hance the robustness of the clock [15].The PPC has been characterized experimentally inconsiderable detail. It is known, for example, that KaiCforms a homo-hexamer [15], KaiA a dimer [15], and KaiBa dimer [5, 15] or a tetramer [17]. KaiC has two phospho-rylation sites per protein monomer, which are phospho-rylated and dephosphorylated in a definite sequence as aresult of KaiC’s autokinase and autophosphatase activ-ity [7, 8]. KaiA stimulates KaiC phosphorylation [20, 21],while KaiB negates the effect of KaiA [20–23]. Thanksto the wealth of available experimental data, the PPChas proven an unusually fruitful system for mathemati-cal modeling [1, 7, 16, 24–26, 28–31].The TTC appears to encompass several different regu-latory pathways and is much less well understood. Multi-ple lines of evidence indicate that circadian promoter ac-tivity can be achieved without any specific cis -regulatoryelement [33–35], suggesting that the clock may regu-late transcription at least in part through changes inchromosome compaction and superhelicity [34, 36, 37].Several proteins important for transcriptional regula-tion of the kaiBC operon have also been identified. Inparticular, LabA and the sensor histidine kinase CikAhave been implicated in negative feedback of phospho-rylated KaiC on its own transcription during the sub-jective night [11, 13, 14, 38, 39], while the SasA histi-dine kinase, acting through the RpaA response regulator,seems to mediate positive feedback during the subjec-tive day [40, 41]. Epistasis analysis suggests that LabAmay likewise act through RpaA [38]. The picture thatemerges is thus that during the phosphorylation phase ofthe PPC, KaiC activates kaiBC expression, while in thedephosphorylation phase, it represses kaiBC expression[39] (Fig. 1).In this study, we perform stochastic simulations to in-vestigate the roles of the PPC and the TTC. We firstexamine a situation in which the total number of eachKai protein is constant —they are neither produced nordestroyed—and only the PPC is operative. We show thatin this case the PPC is highly robust against noise aris-ing from the intrinsic stochasticity of chemical reactions.Even for reaction volumes smaller than the typical vol-ume of a cyanobacterium, the correlation time is longerthan that observed experimentally [9]. Living cells, how-ever, constantly grow and divide, and proteins must thusbe synthesized to balance dilution. In fact, dilution canbe thought of as introducing an effective protein degra-dation rate set by the cell doubling time. We thereforenext study a PPC in which the Kai proteins are producedand degraded with rates that are constant in time. Thesimulations reveal that protein synthesis and decay dra-matically reduce the viability of the PPC; we predict thatfor a cell doubling time of 24 hours and a bacterial volumeof 1 µ m , the PPC dephases in roughly 10 days, muchfaster than real S. elongatus [9]. The same loss of oscil-lations is seen even in the limit of infinite volume, indi-cating it is not due solely to noise. Instead, the constantsynthesis of proteins, which we assume are all initiallycreated in the same phosphorylation state, necessarily in-
FIG. 1: A cartoon of how the TTC might interact with thePPC [39]. Active RpA activates kaiBC expression. Duringthe phosphorylation phase, KaiC activates RpaA, while dur-ing the dephosphorylation phase, KaiC represses RpaA. jects KaiC hexamers with the “wrong” phosphorylationlevels at certain phases of the cycle [1, 42]; if these appearfast enough, they can destroy the oscillation. One roleof the TTC is thus to introduce proteins only when thephosphorylation state of the freshly made KaiC matchesthat of the PPC. Our simulations of a combined PPC andTTC reveal that a TTC can indeed greatly enhance therobustness of the PPC, yielding correlation times consis-tent with those measured experimentally [9]. Finally, weconsider whether the PPC is needed at all, or whether onecould build an equally good circadian clock using only aTTC. We find that it is possible to construct a TTC witha period of 24 hours and the observed correlation timeof a few months [9]. However, this comes at the expenseof very high protein synthesis and decay rates, which im-pose an extra energetic burden on the cell. Our resultsthus suggest that a PPC allows for a more robust oscilla-tor at a lower cost. Although our models are simplified,we argue in the Discussion section that our qualitativeresults are unavoidable consequences of the interactionbetween a circadian clock and cell growth and so shouldhold far more generally.
II. RESULTSA. A protein phosphorylation cycle with constantprotein concentrations is highly robust
In recent years, the PPC has been mathematicallymodelled with considerable success [1, 7, 16, 24–26, 28–31] (for a review, see [43]). All models endow KaiC withan intrinsic ability to cyclically phosphorylate and de-phosphorylate itself, but they differ in how the cycles ofthe individual KaiC proteins are synchronized. In thismanuscript, we adopt the model developed by us [1],which is one of several based on a mechanism that wecalled “differential affinity” [1, 7, 24, 28]: KaiA stim-ulates KaiC phosphorylation, but the limited supply ofKaiA dimers binds preferentially to those KaiC moleculesthat are falling behind in the cycle, allowing them tocatch up. Specifically, in our model each KaiC hexamercan switch between an active conformational state C i ,where the number i of phosphorylated monomers tendsto increase, and an inactive state e C i , where i tends todecrease (Fig. 1); KaiA stimulates phosphorylation ofactive KaiC, but is sequestered by complexes containingKaiB and inactive KaiC. KaiC in the inactive state canthus delay the progress of fully dephosphorylated hexam-ers that have already switched back to the active stateand are ready to be phosphorylated again. With A and Bdenoting, respectively, a KaiA dimer and a KaiB dimer,the model becomes:C i f i ⇄ b i e C i , C i + A k Af i ⇄ k Ab i AC i k pf → C i +1 + A (1) e C i + B k Bf i ⇄ ˜ k Bb i B e C i , B e C i + B ˜ k Bf i ⇄ k Bb i B e C i (2)B x e C i + A x ˜ k Af i ⇄ ˜ k Ab i AB x e C i , AB e C i + A ˜ k Af i ⇄ k Ab i A B e C i (3)C i k ps ⇄ k dps C i +1 , e C i ˜ k ps ⇄ ˜ k dps e C i +1 (4)B x e C i ˜ k ps ⇄ ˜ k dps B x e C i +1 , A y B x e C i ˜ k ps ⇄ ˜ k dps A y B x e C i +1 . (5)This model reproduces the phosphorylation behaviourof KaiC in vitro not only when all Kai proteins arepresent, but also when KaiA and/or KaiB are absent [1].It moreover correctly predicted the experimentally ob-served disappearance of oscillations when the KaiA con-centration is raised [6, 7], a success that strongly supportsthe idea that KaiA sequestration is the primary driver ofsynchronization. Our model does not feature monomerexchange between KaiC hexamers, an alternative meansof synchronisation [25] that has been observed in experi-ments [5, 29]; we and others find that monomer exchangeis not critical for stable oscillations [1, 7, 28]. In the Sup-porting Information ( SI ), we show that similar results areobtained with a model that focuses on the phosphoryla-tion cycle of individual KaiC monomers [7, 8] rather thanof KaiC hexamers and that includes a positive feedbackloop.We quantify our model’s robustness to chemical noiseby performing kinetic Monte Carlo simulations of thechemical master equation [4]. In our simulations, wevary the reaction volume, but keep the concentrations ofthe Kai proteins constant at levels comparable to thoseused in the in vitro experiments [5, 6]. Fig. 2A showsa time trace of the fraction p ( t ) of phosphorylated KaiCmonomers at a volume of order that of a cyanobacterium,while Fig. 2B gives, as a function of volume, the corre-lation number of cycles n / , defined as the number of P ho s pho r y l a t i on R a t i o Time [h]V=1.0 µ m V=0.1 µ m -1 C o rr e l a t i on N u m be r o f C yc l e s Volume V [ µ m ] FIG. 2: The KaiC phosphorylation oscillator is highly stablewhen proteins are not synthesized or degraded. (A) Timetrace of the phosphorylation ratio p ( t ), i.e. the fraction ofmonomers that are phosphorylated, for volumes V = 1 µ m and V = 0 . µ m . (B) Correlation number of cycles, n / , as afunction of the volume. The protein concentrations are thoseused in the in vitro experiments [5, 6]: [A] T = 0 . µM ; [B] T =1 . µ M; [C] T = 0 . µ M. For other parameters, see Table S1of SI . cycles after which the standard deviation in the phase ofthe oscillation is half a day [16]; in the SI we discuss ourerror estimates, which suggest that our computed n / values are typically accurate to within 15–20%. One is-sue that arises in comparing the simulation results tothe measured in vivo clock robustness is that the Kaiproteins appear to be present in living cells in a ratioat which the in vitro system would not oscillate. Ki-tayama et al. found that there are approximately 10,000KaiC monomers but only 250–500 KaiA monomers per S.elongatus cell [22], corresponding to a roughly 6:1 KaiChexamer to KaiA dimer ratio instead of the standard1:1 ratio in vitro [5, 6]. It has been suggested that thisdiscrepancy may indicate that the clock reactions are infact confined to a subdomain of the whole cell from whichsome of the KaiB and KaiC molecules are excluded [22];here, we adopt this hypothesis and assume that the Kaiproteins are found in the physiologically relevant reac-tion volume in proportions comparable to those used inthe in vitro experiments. If we take this volume to be ∼ µ m , comparable to the size of the entire cell, thenFig. 2B shows that n / ≈ ±
100 [9]. Even if we consider a volume ofonly 0 . µ m —small enough that the measured numberof KaiA molecules is adequate to give the in vitro KaiAdimer concentration of 0 . µ M—we find that n / ≈ SI ).Our model thus predicts that the PPC is very resistantto noise arising from the intrinsically stochastic nature ofchemical reactions. B. A phosphorylation cycle with constant proteinsynthesis and degradation rates is not stable
Fig. 2 shows that the phosphorylation cycle is highlyrobust when the total concentrations of the Kai proteinsare strictly constant. But in vivo proteins are constantlybeing synthesized and degraded. To study how this af- P ho s pho r y l a t i on R a t i o Time [h] µ = 0.01 1/h µ = 0.05 1/h µ = 0.1 1/h FIG. 3: A system in which the Kai proteins are producedand degraded with rates that are constant in time cannotsustain a stable phosphorylation cycle. (A) Time trace ofthe phosphorylation ratio p ( t ) for three different degradationrates µ and V = 1 µ m . (B) Correlation number of cycles, n / , as a function of the volume and degradation rate. Thewhite line denotes the bifurcation line µ = 0 . − , wherethe system undergoes a supercritical Hopf bifurcation in thedeterministic limit; the three colored dots correspond to thethree curves in panel A. All proteins are degraded with thesame rate; KaiA, KaiB and KaiC are produced with ratessuch that the average total concentrations equal those usedin the in vitro experiments (see Fig. 2). fects the PPC, we consider a model in which the Kai pro-teins are produced and degraded in a stochastic (memo-ryless) fashion with rates that are constant in time, withthe effects of active degradation [2] and of passive dilu-tion lumped into a single first-order decay rate µ (see SI ).Fig. 3 shows the performance of the phosphorylationcycle as a function of degradation rate and cell volume;in all cases, the synthesis rates are adjusted so that themean concentrations are constant and equal to thoseused in the previous section. The oscillator’s robust-ness clearly decreases dramatically with increasing pro-tein synthesis and decay rate. For a volume comparableto that of a cyanobacterium and a degradation rate of0 . − , the correlation time is less than 20 days, muchlower than that observed in vivo [9]. This degradationrate is precisely the effective rate arising from protein di-lution with a cell doubling time of 24 hours. It is known,however, that KaiC is also degraded actively at a rate ashigh as 0 . − [2], leading to still worse stability.The disappearance of the oscillations for higher proteinsynthesis and decay rates can be understood by notingthat fresh KaiC hexamers are made in a fixed phosphory-lation state, which then has to catch up with those of theproteins that are already in the cycle [1, 42]. When thedegradation rate is high, the new proteins are likely tobe degraded before the PPC can synchronise their phos-phorylation levels; indeed, in the limit that the proteinsynthesis and decay rates go to infinity, p ( t ) becomes con-stant in time and equal to the phosphorylation level offreshly made KaiC proteins. This is not a purely stochas-tic effect; the dashed bifurcation line of Fig. 3B showsthat even in a deterministic model the oscillations disap-pear when the synthesis and decay rates become too big(see SI ). C. A protein phosphorylation cycle with atranscription-translation cycle is very stable
To sustain a phosphorylation cycle, KaiC has to bemade in an oscillatory fashion: newly synthesized KaiCproteins have to be injected into the phosphorylationcycle at the moment that the phosphorylation states ofthe other proteins that are already in the cycle matchthe modification state of the fresh KaiC proteins [1].This is the principal role of the transcription-translationcycle. Here, we present a model for how such a cyclemight interact with the protein phosphorylation cycle(see Fig. 1). The model is inspired by that of Kondo etal. [39] and contains the following key ingredients:
1. RpaA activates kaiBC expression
Deletion of rpaA reduces the expression of clock-controlled genes, including that of kaiBC [38, 39, 41].Since neither promoters nor transcription orchromosome-compaction factors have been identifiedthat interact with RpaA [41], we make the phenomeno-logical assumption that RpaA directly activates kaiBC expression [39].
2. RpaA is activated by KaiC when KaiC is in the activestate
RpaA is activated via phosphorylation by the histidinekinase SasA, whose activity is in turn stimulated by KaiC[36, 41]; inactivation of SasA reduces kaiBC expression[38, 39, 41]. Moreover, SasA-RpaA phosphorylation oc-curs 4-8 hours before the peak of KaiC phosphorylation[41]. This suggests that partially phosphorylated KaiCthat is on the active branch activates RpaA throughSasA [39]. Since SasA phosphorylation, occurring ontime scales of minutes [36], is much faster than KaiCphosphorylation, occurring on time scales of hours, weassume that the SasA dynamics can be integrated out.
3. RpaA is inactivated by KaiC when KaiC is in theinactive state
Inactivation of LabA [38] or CikA [39] increases kaiBC expression, with inactivation of both having a stillstronger effect [39]. SasA inactivation can compen-sate for both LabA [38] and CikA inactivation [39],but RpaA inactvation cannot compensate for LabAinactivation [38]. Taken together, these results suggestthat SasA, LabA, and CikA control kaiBC expressionthrough different pathways, with at least the SasA andLabA pathways converging on RpaA [39]. Becausephosphorylation of KaiC is critical for negative feedbackon kaiBC expression [14], LabA and CikA appear toact downstream of phosphorylated KaiC [39]. Since themechanisms by which LabA and/or CikA repress RpaAactivation are unknown, we make the phenomenologicalassumption that KaiC on the inactive branch deactivatesRpaA.
4. KaiC is injected into the system as fully phosphory-lated hexamers
Imai et al. reported that newly synthesized KaiC isphosphorylated in vivo within 30 min [2], much fasterthan phosphorylation of KaiC hexamers in vitro , whichtakes about 6 hours [7]; we thus assume that newlysynthesized KaiC is injected into the system as fullyphosphorylated KaiC hexamers.
5. KaiA and RpaA are synthesized at constant rates
The mRNA levels of kaiA and rpaA exhibit muchweaker oscillations than those of kaiB and kaiC [33]. Wetherefore assume that kaiA and rpaA are expressed atconstant rates.
6. The phosphorylation cycle in vivo is similar to thatin vitro .Fig. 1 shows a cartoon of this model, which isdescribed by the reactions of Eqs. 1–5 for the PPCtogether with the following reactions for the TTC andthe coupling between them: e R + X k a → R + X , R + e X k i → e R + e X (6) ∅ β c [R] / ( K +[R] ) = ⇒ τ ; σ τ C + 3B (7) ∅ β a → A , ∅ β r → e R , (8)R , e R , A , B , AC i , B x e C i , A y B x e C i µ → ∅ (9)Eq. 6 models activation of inactive RpaA, e R, by X ∈{ AC , . . . AC } and inactivation of active RpaA, R, by e X ∈ { A y B x e C . . . A y B x e C } . The model is robust to theprecise choice of phosphoforms that activate and repressRpaA (see SI ). Because KaiA enhances kaiBC expres-sion [20], we assume that active KaiC bound to KaiAactivates RpaA; a model in which both C i and AC i acti-vate RpaA also generates oscillations, albeit less robustly(see SI ). Eq. 7 models the activation of kaiBC expres-sion by RpaA, using a Hill function with coefficient 4. Weassume a normally distributed delay, denoted by the dou-ble arrow, with mean τ = 5 hr and standard deviation σ τ = 0 . kaiBC transcrip-tion and the appearance of KaiB and KaiC protein. Thelength of the delay is dictated by the requirement thatnew KaiC protein be produced when its phosphorylationstate matches that of the PPC. Eq. 8 models constitutiveexpression of kaiA and rpaA . Eq. 9 describes degradationof all species with the same rate constant µ ; we ignorerhythmic KaiC degradation [2], which is not essential toproduce a robust clock. In the SI we discuss a more de-tailed model that includes cell growth and binomial par-titioning upon cell division, which appreciably increasesnoise without changing qualitative trends. Although onecan think of our model in loose terms as consisting of cou-pled transcriptional and post-translational oscillators, itcannot be mathematically decomposed into two separateoscillatory systems, each with its own variables, nor canthe strength of the coupling between the two cycles beindependently tuned. It is thus formally quite differentfrom textbook models of coupled oscillators [10].Fig. 4B shows the robustness of this PPC-TTC modelas a function of cell volume and protein degradation rate µ ; when µ is varied, all three protein synthesis rates β α C on c en t r a t i on [ µ M ] P ho s pho r y l a t i on R a t i o Time [h]Total KaiCPhosphorylation Ratio -1 Volume V [ µ m ]10 -3 -2 -1 D eg r ada t i on R a t e µ [ / h ] C o rr e l a t i on N u m be r o f C yc l e s FIG. 4: The combination of a TTC and a PPC can gener-ate stable circadian rhythms. (A) Time trace of the phos-phorylation ratio p ( t ) and total KaiC concentration [C] T for V = 1 µ m and µ = 0 . − (solid lines) and 0 . − (dashedlines). (B) Correlation number of cycles, n / , for p ( t ) as afunction of the volume and degradation rate µ . The proteinsynthesis rates are varied with the degradation rates such thatthe average protein concentrations equal those used in vitro (see Fig. 2). are adjusted so that the average protein concentrationsare unchanged and comparable to those of the modelsconsidered above. As expected, n / decreases with de-creasing cell volume. Its dependence on the degradationrate, however, is markedly different from that seen withconstant KaiC synthesis (Fig. 3B): A PPC sustained bya TTC becomes more robust with increasing decay rate.If we assume that proteins are lost only through dilu-tion, then for a bacterial volume of 1 µ m and a celldoubling time of 24 hours (corresponding to a decayrate µ = 0 . − ), the correlation time is about 200days, consistent with the value measured experimentally[9]. If proteins are also degraded actively, increasing µ ,this excellent behavior improves still further; even for V = 0 . µ m , n / ≈
120 for µ = 0 . − . This is instark contrast to the stability of a PPC without a TTC,which we have seen is hardly stable for a protein lifetimeof about 10 hours (Fig. 3B). D. A protein phosphorylation cycle dramaticallyenhances the robustness of atranscription-translation cycle
Figs. 3B and 4B show that a TTC can greatly improvethe robustness of a PPC. One might thus ask whetherthe PPC is needed at all, or whether an adequate clockcan be built with only a TTC. To address this question,we modify the PPC-TTC model (Eqs. 6–9) so that itconsists only of a TTC: ∅ β c [R] / ( K +[R] ) = ⇒ τ ; σ τ C , ∅ β r → R , C , R µ → ∅ (10) e R k ta → R , C + R k ti → C + e R (11)In the simulations, we adjust the delay τ and the synthe-sis rate β for each choice of the degradation rate µ suchthat the oscillation period is 24 hours and the averageKaiC concentration is comparable to that of the othermodels considered so far (see SI ). We also investigatedthe simplest possible TTC model, namely one in whichKaiC directly represses its own synthesis; this gave verysimilar results (see SI ).Fig. 5B shows for a bacterial volume of approximately1 µ m the robustness of this TTC-only model as a func-tion of the degradation rate, together with the resultsfor the PPC-TTC model (Fig. 4B) and for the PPC-only model with constant synthesis rates (Fig. 3B). TheTTC-only model’s behavior is the opposite of that ofthe PPC-only model; the TTC-only model is most stablefor high degradation rates, and its robustness falls dra-matically when µ drops below 0 . − . The combinedTTC-PPC model, however, is robust for all degradationrates—its correlation time interpolates between those ofthe TTC-only model for large µ and of the PPC-onlymodel for small µ . Importantly, the relative advantageof the combined model is greatest when 1 /µ is of order theoscillator period, and this is precisely the regime wherephysiological degradation and dilution rates for S. elon-gatus fall [2]. Further, the combined oscillator does muchbetter than either oscillator alone; in contrast, when twoconventional phase oscillators with comparable noise lev-els are coupled, one expects only about a factor of twogain in n / [10].Dilution puts a lower bound on the degradation rate,which means that a stable oscillator cannot be basedon a PPC only, especially when the growth rate of thebacterium is large. The degradation rate can, on theother hand, be increased by active degradation. For highdegradation rates, i.e. µ > . − , an oscillator basedonly on transcriptional feedback can be sufficiently ro-bust (Fig. 5B). However, to balance these high decayrates, the protein synthesis rates have to be correspond-ingly large, which is energetically costly [47]. Combin-ing a PPC with a TTC makes it possible to dramaticallyimprove the robustness while keeping the protein synthe-sis rates the same. While a phosphorylation cycle doesnot come entirely for free [48], this suggests that a PPCcombined with a TTC gives the best performance-to-costratio.The question remains why a PPC so effectively en-hances the robustness of a TTC. Fig. 5A shows, for thefull PPC-TTC model, time traces of the concentration ofRpaA, the sum of the concentrations of the KaiC phos-phoforms that activate RpaA, and the sum of those thatrepress RpaA; for comparison, we also show the KaiCconcentration in the TTC-only model. In the TTC-onlymodel, the KaiC concentration varies slowly, and the con-centration threshold for kaiBC repression is crossed at ashallow angle. In contrast, in the PPC-TTC model theconcentration of active RpaA switches rapidly betweena value that is well above the activation threshold for kaiBC expression and one that is well below it. Suchsharp threshold crossings minimize the effect of fluctua-tions in the concentration of the regulatory protein on thetiming of the switch between kaiBC activation and kaiBC repression. RpaA’s abrupt switching is a direct conse-quence of the sharp rise and fall in the concentrations ofthe KaiC phosphoforms that regulate RpaA. This in turn C on c en t r a t i on [ µ M ] Time [h]Active RpaARepressorActivatorTTC -3 -2 -1 C o rr e l a t i on N u m be r o f C yc l e s Degradation Rate µ [1/h]PPCTTCFull FIG. 5: A PPC in combination with a TTC generates ro-bust rhythms over a wide range of degradation rates. (A)Time traces of the concentrations of active RpaA and of theKaiC phosphoforms that activate (“activator”) and repress(“repressor”) RpaA for the PPC-TTC model ( V = 1 µ m and µ = 0 . − ), and KaiC for the TTC-only model ( V = 1 µ m and µ = 0 . − ; for µ = 0 . − the TTC does not os-cillate); (B) Correlation number of cycles as a function ofdegradation rate µ for the PPC, TTC, and combined TTC-PPC model ( V = 1 µm ); for the TTC, n / was computedfor [C]( t ) rather than for p ( t ). The protein synthesis ratesare varied with the degradation rates such that the averageconcentrations equal those used in the in vitro experiments(see Fig. 2). stems not so much from the oscillation in the total con-centration of KaiC, but from the phosphorylation cycleat the multiple phosphorylation sites of KaiC: The con-centrations of the various phosphoforms approach zero incertain parts of the cycle (Fig. 5A), even though the totalKaiC concentration remains relatively large throughout(Fig. 4A). Indeed, the stability of the circadian rhythm inthe PPC-TTC model is mostly determined by the KaiCphosphorylation cycle, which, as Fig. 2B shows, is highlyrobust. It thus seems that the phosphorylation cycle isthe principal pacemaker, even though a transcription-translation cycle is needed to sustain it and can, in fact,raise its robustness above the limiting value at zero pro-tein synthesis and decay rate (Fig. 5B). III. DISCUSSION
The evidence is accumulating that circadian rhythmsare driven by both transcription-translation and protein-modification cycles not just in cyanobacteria but even inhigher organisms [1, 49]. Our analysis suggests that bothcycles are required to generate stable circadian rhythmsin growing cells over a wide range of conditions. On theone hand, it is intrinsically difficult for a TTC to generateoscillations when the average protein decay time is longcompared to the oscillation period. Indeed, the speed ofprotein degradation (and/or dilution) sets a simple up-per bound on the size of any oscillation; if only a smallfraction of the proteins can be degraded in one oscillationcycle, then any modulation of the total protein concen-tration must necessarily have a low amplitude and so bevery susceptible to noise. In cells whose doubling time isof order tens of hours (comparable to the 24 hour clockperiod) or longer, a TTC alone can thus function robustlyonly at the expense of high rates of protein synthesis andactive protein degradation. On the other hand, a pro-tein modification cycle must fail when protein turnoverbecomes faster than the clock period: Proteins are thentypically degraded and replaced by new molecules in thewrong modification state before they can pass throughthe full cycle, destroying the oscillation. While a TTC isthus required to sustain a protein modification cycle athigher growth rates, a modification cycle can also greatlyenhance the performance of a TTC, especially at lowerprotein synthesis and decay rates, where the modifica-tion cycle can be exploited to set the rhythm of geneexpression. The protein phosphorylation cycle of KaiCis intrinsically very robust, because KaiC is present inrelatively large copy numbers and because a KaiC hex-amer must go through many individual phosphorylationand dephosphorylation steps in each oscillation period.These two characteristics together are responsible for thereliable and pronounced oscillations in the concentrationsof the various phosphoforms of KaiC (Fig. 5A), which,in turn, drive the robust rhythm of kaiBC expression.Our combined PPC-TTC model is consistent witha number of experimental observations. It not onlymatches the average oscillations of the phosphorylationlevel and the total KaiC concentration in wildtype cells(Fig. 4A), but also reproduces the observation that evenin the presence of an excess of KaiA, the total KaiC con-centration still oscillates with a circadian period [15] (see SI ). Yet, because quite a few elements of the TTC havenot been characterized experimentally, our model of theTTC is necessarily rather simplified and phenomenolog-ical; not surprisingly, some observations thus cannot bereproduced. For example, our model predicts that thephase of the oscillation in KaiC abundance lags behindthat of the phosphorylation level by a few hours, while ex-periments seem to show that these oscillations are morein phase [7, 15]. This may be due to our simplifyingassumption that KaiC is produced as fully phosphory-lated hexamers. While phosphorylation of fresh KaiC hasbeen reported to occur within 30 minutes [2], fragmen-tary evidence suggests that hexamerisation is slow [50];one would then expect to detect KaiC monomers in ex-periments several hours before our simulations report thepresence of fully functional hexamers. We believe, how-ever, that while our model is simplified, it does capturethe essence of a coupled TTC and PPC: that is, it isbuilt around a protein that not only regulates its ownsynthesis, but also undergoes a protein modification cy-cle, where the different modification states either activateor repress protein synthesis. Moreover, the basic trendswe observe can largely be explained by the argumentsof the preceding paragraph, which depend only on thecomparison of a protein decay timescale with the oscilla-tion period; we thus expect that our qualitative resultsshould apply to any system that exploits both a protein modification cycle and a transcription-translation cycleto generate circadian rhythms.In S. elongatus , the period of the circadian rhythm isinsensitive to changes in the growth rate [51, 52]. Sinceprotein synthesis and decay rates, on the other hand,tend to vary with the growth rate [53], the question ariseswhether the period of the oscillator as predicted by ourmodel is robust to such variations. In the SI we showthat a ten-fold increase in the synthesis and degradationrates of all proteins only decreases the period by 10-20%.The reason is that the period of the oscillator is mostlydetermined by the intrinsic period of the phosphoryla-tion cycle (Fig. 5A). The latter does not depend on theabsolute protein synthesis and decay rates, although itdoes depend on the ratio of the concentrations of the Kaiproteins [1, 5, 6]. Clearly, it would be of interest to in-vestigate how the ratio of the concentrations of the Kaiproteins varies with the growh rate.Finally, how could our predictions be tested experi-mentally? The most important testable prediction ofour analysis is that the PPC requires a TTC when thegrowth rate is high. Kondo and coworkers have demon-strated that a PPC can function in the absence of a TTC[7]. These experiments were performed under constantdark (DD) conditions, which means that the growth ratewas (vanishingly) small, protein synthesis was halted,and even protein decay was probably negligible, sincethe KaiC level was rather constant [7]. These experi-ments are consistent with our predictions, but the criti-cal test would be to increase the growth rate while avoid-ing any circadian regulation of kaiBC expression. Thismeans that the cyanobacteria must be grown under con-stant light (LL) conditions. Moreover, one would need tobring kaiBC expression under the control of a promoterthat is constitutively active. However, most promotersof S. elongatus [33, 34], and even many heterologous pro-moters from
E. coli [35], are influenced by the circadianclock. Nevertheless, a number of promoters have beenreported that exhibit arhythmic activity [33], and thesemight be possible candidates. A still more challengingexperiment would be to express the phosphorylation cy-cle in the bacterium
E. coli [54]. Our analysis predictsthat in growth-arrested cells, the phosphorylation cycleshould be functional, while in normal growing
E. coli cells, with cell doubling times of roughly 1 hour, the os-cillations should cease to exist.
A. Acknowledgments
We thank Tom Shimizu for critically reading themanuscript. This work was supported in part byFOM/NWO (D.Z. and P.R.t.W.) and by NSF grantPHY05-51164 to the KITP (D.K.L.). [1] Johnson CH, Egli M, Stewart PL (2008) Structural in-sights into a circadian oscillator.
Science
Science
Cell
Science
Na-ture
Proc Natl Acad SciUSA
Science
Science
Genes Dev
Synchro-nization, a universal concept in nonlinear sciences (Cam-bridge University Press).[11] Ishiura M, et al. (1998) Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria.
Science
EMBO J
ProcNatl Acad Sci USA
Synechococcus elongatus
PCC 7942.
Proc Natl Acad Sci USA
J Biol Chem
Mol Cell
J Biol Chem
EMBO J
Science
Proc NatlAcad Sci USA kaiBC promoter in regulating KaiC.
EMBO J
EMBO J
Synechococcus elongatus : a potential clock inputmechanism.
Proc Natl Acad Sci USA
J Biol Rhythms
Phys Rev Lett
PLoS Comput Biol
Proc Natl Acad Sci USA
Mol Syst Biol
PLoS Biol
J Biol Rhythms
PLoS One
Biophys J
Synechococcus elongatus . Proc Natl Acad SciUSA
Proc Natl Acad Sci USA
Annu Rev Genet
Synechococcus elongatus . ProcNatl Acad Sci USA
Proc Natl Acad Sci USA labA : a novel gene required fornegative feedback regulation of the cyanobacterial circa-dian clock protein KaiC.
Genes Dev based oscillator cooperate to generate robust circadian kaiBC expression in cyanobacteria.
Proc Natl Acad SciUSA .[40] Iwasaki H, et al. (2000) A KaiC-interacting sensory his-tidine kinase, SasA, necessary to sustain robust circadianoscillation in cyanobacteria.
Cell
Proc Natl Acad SciUSA
Nat Struct MolBiol
FEBS Lett
FEBS Lett
J. Phys. Chem.
J Biol Chem
Nature
Proc Natl Acad Sci USA
Curr Biol
ProcNatl Acad Sci USA
Proc Natl Acad SciUSA
Science
Cell
Int JBioinform Res Appl
Supporting Information
Robust circadian clocks from coupled protein modification and transcription-translation cycles
David Zwicker, David K. Lubensky and Pieter Rein ten Wolde
I. THE MODELA. Protein phosphorylation cycle
The model for the KaiC protein phosphorylation cycle (PPC) is described by reactions (1)-(5) of the main text,which we repeat here:C i f i ⇄ b i e C i , C i + A k Af i ⇄ k Ab i AC i k pf → C i +1 + A , (S1) e C i + B k Bf i ⇄ ˜ k Bb i B e C i , B e C i + B ˜ k Bf i ⇄ k Bb i B e C i , (S2)B x e C i + A x ˜ k Af i ⇄ ˜ k Ab i AB x e C i , AB e C i + A ˜ k Af i ⇄ k Ab i A B e C i , (S3)C i k ps ⇄ k dps C i +1 , e C i ˜ k ps ⇄ ˜ k dps e C i +1 , (S4)B x e C i ˜ k ps ⇄ ˜ k dps B x e C i +1 , A y B x e C i ˜ k ps ⇄ ˜ k dps A y B x e C i +1 . (S5)Here, C i denotes a KaiC hexamer in the active conformational state, in which the number i of phosphorylatedmonomers tends to increase, and e C i denotes a KaiC hexamer in the inactive conformational state in which i tendsto decrease; A denotes a KaiA dimer, and B denotes a KaiB dimer. The reactions C i ⇄ e C i in Eq. S1 modelthe conformational transitions between active and inactive KaiC; the second set of reactions in Eq. S1 describephosphorylation of active KaiC that is stimulated by KaiA; the reactions in Eqs. S2 and S3 model the sequestrationof KaiA by inactive KaiC that is bound to KaiB; note that an inactive KaiC hexamer can bind up to two KaiAdimers; the reactions in Eqs. S4 and S5 model spontaneous phosphorylation and dephosphorylation of active andinactive KaiC. For a more detailed discussion of the model, we refer to [1]. B. A PPC with constant protein synthesis and degradation rates
In the main text, we also discuss the performance of a model that includes not only a PPC, but also synthesis anddegradation of the Kai proteins. These reactions are given by ∅ β c → C + 3B , ∅ β a → A , (S6)A , B , AC i , B x e C i , A y B x e C i µ → ∅ (S7)As explained in the main text, we assume that fresh KaiC is injected into the system as fully phosphorylated hexamersbecause phosphorylation of fresh KaiC proteins has been reported to be fast, i.e. occurring within 30 minutes [2].However, the precise choice for the phosphoform of fresh KaiC is not so important in this model; it does not affectthe robustness of this model. Because the KaiB and KaiC proteins are both products of the kaiBC operon, wechoose to model the production of both proteins as a single reaction. We note that while in the model with thetranscription-translation cycle, discussed in section I D, the delay in the synthesis reactions is critical, in the abovemodel, where the Kai proteins are produced with rates that are constant in time, a delay would have no effect; thesynthesis reactions are therefore modeled as simple, Poissonian birth reactions. C. Deterministic PPC model with synthesis and degradation
To verify that the disappearance of oscillations as the degradation rate is increased is not a purely stochastic effect,we consider the model of Eqs. (S1)–(S7) in the deterministic limit of infinite volume and protein number. In this limit,the concentrations of the different proteins evolve according to deterministic rate equations. We make two furthersimplifying assumptions: First, we replace the two-step binding of KaiB to e C i with a trimolecular reaction that turns e C i directly into B e C i , and making a similar change for binding of KaiA to the inactive branch. Second, we assumethat binding and unbinding reactions are fast enough that they are effectively in steady state and thus explicitly keeptrack only of the concentrations of the various KaiC species; the concentrations of free KaiA and KaiB can then beinferred from conservation laws. The dynamical equations are then essentially the same as those given in Eqs. 44 – 47of our previous publication [1], with the addition of a linear decay term with rate µ for each species and of synthesisof C with rate β c : d [C i ] T dt = σ ps i − [C i − ] T + σ dps i +1 [C i +1 ] T − ( σ ps i + σ dps i )[C i ] T − σ Ff i [C i ] T + σ Fb i [ e C i ] (S8)+ β c δ i, − µ [C i ] T (S9) d [ e C i ] dt = ˜ k ps [ e C i − ]+˜ k dps [ e C i +1 ] − (˜ k ps +˜ k dps )[ e C i ]+ σ Ff i [C i ] T − σ Fb i [ e C i ] − ˜ κ Bf i ([B] T − P i [B e C i ] T ) [ e C i ] + ˜ κ Bb i e K i [B e C i ] T e K i + [A] − µ [ e C i ] (S10) d [B e C i ] T dt = ˜ k ps [B e C i − ] T +˜ k dps [B e C +1 ] T − (˜ k ps +˜ k dps )[B e C i ] T +˜ κ Bf i ([B] T − P i [B e C i ] T ) [ e C i ] − ˜ κ Bb i e K i [B e C i ] T e K i + [A] − µ [B e C i ] T , (S11)where the concentration of free KaiA, [A], is given by[A] + N X i =0 [A][C i ] T K i + [A] + 2 N X i =0 [A] [B e C i ] T e K i + [A] − [A] T = 0 . (S12)Here [C i ] T is the total concentration of KaiC hexamers with i phosphorylated monomers in the active state, whether ornot complexed with KaiA, i.e [C i ] T = [C i ] + [AC i ]; [B e C i ] T is defined similarly. The effective rate constants appearingin these equations depend on the concentration [A] of free KaiA and are defined in terms of the more microscopic rateconstants as follows: The effective (de)phosphorylation rates on the active branch are σ ps i = ( k ps K i + k pf [A]) / (K i +[A])and σ dps i = K i k dps / (K i + [A]). The effective flipping rates are given by σ F fi = f i K i / ( K i + [ A ]) and σ F bi = b i , where f i and b i are the forward and backward flipping rates. The parameters ˜ κ Bf i and ˜ κ Bb i differ from ˜ k Bf i and ˜ k Bb i , respectively,in that the κ ’s are rate constants for trimolecular reactions which are broken down into two successive bimolecularreactions in the stochastic simulations. The dissociation constants K i satisfy K i = k Ab i /k Af i ; the e K i could be definedsimilarly in terms of forward and backward rates for KaiA binding to the inactive branch, but (just as with ˜ κ Bf i and ˜ κ Bb i ) these rates would differ from those used in the stochastic simulations, so we choose instead to quote thedissociation constants directly. Following [1], we choose values for the new parameters associated with trimolecularinteractions such that time dependence of the phosphorylation ratio p ( t ) matches the average behavior of the stochasticmodel.To determine where oscillations disappear as µ is increased, we analyzed these equations using the XPPAUTimplementation of the AUTO continuation package [3]. We found that, for the parameter values given in Table S1,the system undergoes a supercritical Hopf bifurcation at µ = 0 . / h , as noted in the main text. D. The PPC and TTC combined
The full model of the main text consists of a PPC, a transcription-translation cycle (TTC), and a pathway thatcouples these two cycles. This model is described by the reactions of Eqs. S1-S5 for the PPC, together with the
Constant ValuePPC (Eqs. 1-5 and S1-S5): k ps , ˜ k ps / h k dps , ˜ k dps / h k pf / h f i { − , − , − , − , − , − , } / h b i / h k Af i . ·
10 1 / M · h k Ab i { , , , , , , } / h ˜ k Bf i . · × { . , . , , , , , } / M · h ˜ k Bb i { , , , , , , } / h ˜ k Af i . · × { . , , , , , . , . } / M · h ˜ k Ab i { , , , , , , } / h Deterministic PPC (Eqs. S8-S12):˜ κ Bf i . · × { . , , , , , , } / M · h ˜ κ Bb i × { , , , , , , } / h K i . · − × {∞ , , , , , ∞ , ∞} M RpaA activation (Eqs. 6 + S13): k a . · / M · h k i . · / M · h TTC (Eqs. 7-9, S14,S15): K . µ M β a µ − × . µ M β r µ − × . µ M β c from optimization for h [C] i = 0 . µ M τ σ τ . k ta / M · h k ti / M · h TABLE S1: The parameters used for the models of the main text. The degradation rate µ is a free parameter that we vary toexplore different growth conditions. The numbers between the curly brackets in the right column correspond to the differentKaiC phosphorylation states i in ascending order; values of ∞ for K i indicate that a particular binding reaction is not allowed. following reactions for the TTC and the coupling between them: e R + X k a → R + X , R + e X k i → e R + e X , (S13) ∅ β c [R] / ( K +[R] ) = ⇒ τ ± σ τ C + 3B , ∅ β a → A , ∅ β r → e R , (S14)R , e R , A , B , AC i , B x e C i , A y B x e C i µ → ∅ . (S15)Here, R and e R denote the RpaA protein in its active and its inactive form, respectively. The X and e X in Eq. (S13)denote any of the phosphoforms of KaiC that mediate the activation and repression of RpaA, respectively; in thenext section, we discuss the dependence of the results on precisely which phosphoforms are chosen. The double arrowindicates a reaction with a Gaussian distributed delay with mean τ and variance σ τ . We thus assume that kaiBC expression is activated by RpaA, where the activity of RpaA is modulated by the PPC. In contrast, the expression ofKaiA, KaiB, and RpaA is assumed to occur constitutively. E. Parameters
Table S1 gives the values of the kinetic parameters used in the stochastic simulations based on the kinetic MonteCarlo algorithm developed by Gillespie [4]. Unless otherwise noted, we choose the total concentrations of the Kaiproteins to match common conditions for the in vitro reaction system [5, 6]: [A] T = 0 . µ M, [B] T = 1 . µ M, and[C] T = 0 . µ M. Model Activator Repressor Threshold
K n / a AC , AC , AC , AC A y B x e C , . . . , A y B x e C . µ M 195b AC , AC A y B x e C , . . . , A y B x e C . µ M 118c AC , AC A y B x e C , . . . , A y B x e C . µ M 180d C , AC , C , AC A y B x e C , . . . , A y B x e C . µ M 39e C x , AC x , x ∈ { , , , } A y B x e C , . . . , A y B x e C . µ M 48TABLE S2: Models with different output pathways from the PPC to the TTC. These models differ in the choice of phosphoformsthat activate and repress RpaA, respectively. The maximal production rate β c has been modified such that the average concen-tration of KaiC is 0 . µ M, as used in the in vitro experiments [5, 6]. In each case, n / is given for a volume V = 1 µ m and adecay rate µ = 0 .
03 h − . Model “a” is the full model from the main text. To make the simulations tractable, we neglected repres-sion of RpaA activation by KaiC phosphoforms that occur in negligible concentrations; consequently, the full list of phosphoformsthat has the potential to repress RpaA is { B e C , B e C , AB e C , A B e C , B e C , AB e C , A B e C , B e C , AB e C , A B e C , B e C , AB e C , A B e C } . Other parameters are given in Table S1. II. RESULTS ARE INDEPENDENT OF DETAILS OF THE OUTPUT PATHWAY
In this section, we show that the precise choice of the KaiC phosphoforms that activate and repress RpaA is notcritically important for the existence of oscillations. Table S2 shows the different models that we have considered,and Fig. S1 shows their time traces. It is seen that the time traces are very similar to those of the model in the maintext, which is model “a” in Table S2. The most significant difference can be observed for the time trace of RpaA inmodels “d” and “e”. In these models, not only C x A activates RpaA, but also C x , thus KaiC that is not bound toKaiA. The concentration of C x A reaches zero during the dephosphorylation phase, and, as a result, the concentrationof RpaA becomes zero during this phase in models “a”-”c”. However, the concentration of C x does not reach zeroduring the dephosphorylation phase, and consequently, there is some residual activation of RpaA during this phasein models “d” and “e”. Nevertheless, RpaA activation during the dephosphorylation phase in these models does notmanifest itself in the time traces of KaiC, because the concentration of active RpaA is still below the threshold for kaiBC expression. The oscillations of the phosphorylation level and total KaiC concentration are thus fairly similarin all models, although models “d” and “e” are less robust. (A) P ho s pho r y l a t i on R a t i o Time [h] Model aModel bModel cModel dModel e (B)
700 800 900 1000 1100 1200 1300 1400 0 12 24 36 48 60 72 84 96 T o t a l K a i C C on c en t r a t i on Time [h] Model aModel bModel cModel dModel e (C) C on c en t r a t i on o f A c t i v e R ap A Time [h] Model aModel bModel cModel dModel e
FIG. S1: Time traces of the stochastic simulations of models with different output pathways from the PPC to the TTC (seeTable S2). The phosphorylation ratio (panel A) and the total concentration of KaiC (panel B) are fairly similar for all modelsstudied. -1 Volume V [ µ m ]10 -3 -2 -1 D eg r ada t i on R a t e µ [ / h ] C o rr e l a t i on N u m be r o f C yc l e s FIG. S2: Contour plot of the correlation number of cycles of the PPC model of Rust et al. [7] with production and degradationof Kai proteins with rates that are constant in time, as a function of the volume and the degradation rate.
III. AN ALTERNATIVE PPC MODEL
In the main text, we argue that the synergy between a transcription-translation cycle and a protein modificationcycle is a generic feature of clocks that exploit both cycles. To support this claim, we have studied a model in whichour model of the PPC is replaced by that of Rust et al. [7]. This model describes a phosphorylation cycle at thelevel of KaiC monomers, rather than KaiC hexamers as in our model. In the Rust model, each KaiC monomercycles between an unphosphorylated state “U”, a singly phosphorylated state “T” where KaiC is phosphorylated atthreonine 432, a doubly phosphorylated state “ST” where KaiC is phosphorylated at threonine 432 and serine 431,and a singly phosphorylated state “S” where KaiC is phosphorylated at serine 431 [7, 8]. This cycle is described bythe reactions U ↔ T , T ↔ ST , ST ↔ S , S ↔ U (S16)with reaction rates given by Eq. 5 of the supplementary material of Rust et al. [7]. These rates depend on theconcentration of free KaiA, which is sequestered by KaiC in the S-state. We model KaiA sequestration explicitly:S + A ↔ AS , AS + A ↔ A S . (S17)We picture sequestration to be fast and we picked a forward rate of 0 . / µ M · s and a backward rate of / s for bothequations above. Dephosphorylation of KaiC in the S-state might occur even when KaiA is bound, in which case theKaiA protein is released from the complex. We define the output signal as p ( t ) = [T] + [ST] + [S] + [AS] + [A S][U] + [T] + [ST] + [S] + [AS] + [A S] , (S18)which resembles the phosphorylation ratio in the case where we cannot distinguish between singly and doubly phos-phorylated KaiC. The denominator in the above expression is also the total KaiC monomer concentration. We use thesame concentrations as Rust et al. , [KaiA] = 1 . µ M (active KaiA monomers) and [KaiC] = 3 . µ M (KaiC monomers),and simulate this model using the Gillespie algorithm [4].When we simulate this model for a volume V = 1 µ m , we find a period of L = 21 . τ d = 128 h. The corresponding correlation number of cycles is n / = 30, which is lowerthan that observed experimentally, n / = 166 ±
100 [9], and lower than that of the PPC developed by us [1], for which n / ≈ et al. features a phosphorylation cycle at the level of KaiC monomersrather than KaiC hexamers as in our model. The concomitant reduction in the total number of phosphorylation stepsin the cycle reduces the robustness in the model of Rust et al. [7]. A. The Rust model with constant protein synthesis and degradation
To study the behavior of the model of Rust et al. [7] under conditions in which cells grow and divide, we have toinclude protein degradation and make up for this by protein synthesis. As in the main text, when we vary the protein P r o t e i n C on c en t r a t i on [ µ M ] P ho s pho r y l a t i on R a t i o Time [h]KaiC MonomersPhosphorylation Ratio10x RapA
FIG. S3: Time trace of the Rust model [7] extended to include a TTC for V = 1 µ m . Both the total concentration of KaiCand the phosphorylation level of KaiC show stable oscillations. The messenger protein RpaA concentration is shown tenfoldand crosses the threshold of K = 0 . µ M ( ≈
100 molecules) reliably. degradation rates, we adjust the protein synthesis rates such that the average protein concentrations are unchangedand similar to those used in the in vitro experiments [5, 6].Figure S2 shows the results for this model. They are qualitatively the same as those of Figure 3 of the main text:the robustness decreases with decreasing volume and increasing degradation rate. Hence, not only in our model butalso in that of Rust et al. [7], protein degradation can cause the oscillations to disappear. This supports our claimthat a protein modification oscillator cannot function on its own when the cell’s growth rate is high enough.
B. The Rust model with a TTC
We will now show that a TTC can resurrect the PPC of Rust et al. [7]. We model the TTC as: ∅ β c [R] / ( K +[R] ) = ⇒ τ,σ τ ST , ∅ β a → A , ∅ β r → e R , (S19) e R + T k a → R + T , R + A x S k i → e R + A x S , x ∈ { , , } , (S20)A , R , e R , U , S , A x S , T , ST µ → ∅ . (S21)where the first line describes the production of proteins to counteract their degradation and the second line summarizesthe RpaA signaling pathway, where KaiC that is phosphorylated at the T site activates RpaA and KaiC that isphosphorylated at the S site represses RpaA activation. The parameters in this model are β c = 1 . µ M / h , K =0 . µ M, β a = 0 . µ M / s , β r = 0 . µ M / s , k a = k i = 0 . / µ M · s , and the delay is τ = (3 ± .
3) h. Figure S3 shows theresults of this model at a volume V = 1 µ m . Analysing the autocorrelation function, we find a period L = 22 . τ = 1832 h leading to a correlation number of cycles of n / = 402. Clearly, a TTC canalso resurrect the PPC of Rust et al. [7], supporting our claim that the qualitative results of the main text shouldapply to any biological system that exploits both a protein modification cycle and a protein synthesis cycle to generatecircadian rhythms. IV. THE EFFECT OF BURSTS
In the model of the main text, described in section I of this SI , we have concatenated transcription and translationinto one gene-expression step. Moreover, we have ignored promoter-state fluctuations. Allowing for the explicitformation and translation of mRNA [10], as well as for slow promoter-state fluctuations [11, 12], could lead to burstsin protein synthesis, which are expected to lower the robustness of the TTC. This could potentially lower the stabilityof the clock. To address this, we have performed simulations of a model in which KaiB and KaiC are produced inbursts. We assume that 5 KaiC hexamers and 15 KaiB dimers are formed in each gene expression reaction (rather -1 Volume V [ µ m ]10 -3 -2 -1 D eg r ada t i on R a t e µ [ / h ] C o rr e l a t i on N u m be r o f C yc l e s FIG. S4: The correlation number of cycles as a function of the degradation rate µ and the volume V for a TTC-PPC modelthat exhibits bursts in gene expression; upon a gene expression event, 5 KaiC molecules are produced and 15 KaiB dimers. than the 1 and 3 of Eq. S14); this corresponds to typical burst sizes observed experimentally [10] in E. coli . Eq. S14is thus replaced by: ∅ β c [R] / ( K +[R] ) = ⇒ τ ± σ τ + 15B . (S22)Fig. S4 shows the resulting phase diagram. As expected, it is similar to Fig. 4 of the main text, which shows theresults of the PPC-TTC model without bursts in gene expression. The robustness of the model with bursts is lower,but not very much so: n / = 150 for the model with bursts versus n / = 195 for the model without bursts shownin the main text ( µ = 0 . − and V = 1 µ m in both cases). We believe that this relatively small reduction in theclock’s stability is due to the stabilizing effect of the PPC. V. THE FULL MODEL WITH VOLUME GROWTH AND BINOMIAL PARTITIONING
Living cells constantly grow and divide, and proteins thus have to be synthesized to balance dilution. In themain text, we argued that the principal effect of dilution is to introduce an effective degradation rate set by thecell doubling time. Here, we show that this is indeed the case: we study a model in which growth, cell divisionand binomial partitioning of the proteins upon cell division are modeled explicitly [13], and show that its qualitativebehavior is similar to the model of the main text, in which the volume is held constant and protein degradation occursat a rate that is constant in time.The model we consider here is the combined PPC-TTC model presented in the main text, but with the degradationreaction, Eq. 9, replaced by a scheme in which the bacterial volume V grows exponentially as V ( t ) = V e t ln 2 T d , (S23)where T d denotes the doubling time after which the volume reaches twice its minimum V and cell division is triggered.Division includes dividing the volume by two, partitioning the proteins binomially [13], and deleting events on thequeue of the delay associated with the KaiC production reaction with a probability of 0.5 for each daughter cell.To compare the results of this model with those from the main text, we take T d = ln 2 /µ , where µ is the proteindegradation rate of the model in the main text; if proteins were to decay only by dilution in a cell with a doubling time T d , then µ would be the effective protein degradation rate; if proteins are also degraded actively, then µ = ln 2 /T d isa lower bound on the actual degradation rate.Fig. S5 shows time traces of the total KaiC concentration and the KaiC phosphorylation level for this refined model.It is seen that the oscillations of the total KaiC concentration are more noisy than those in the model in which the Kaiproteins are produced and degraded with rates that are constant in time (Fig. 4A of the main text). Clearly, binomialpartitioning is a major source of noise, with the random removal of items from the queue associated with the KaiCsynthesis reaction being the largest source of noise. Nonetheless, the oscillations of the KaiC phosphorylation level (A) C on c en t r a t i on [ µ M ] P ho s pho r y l a t i on R a t i o Time [h]Total KaiCPhosphorylation Ratio (B) C on c en t r a t i on [ µ M ] C op y nu m be r Time [h]KaiC concentrationVolume [a.u.]KaiC copy number
FIG. S5: Time traces of the full model, given by Eqs. 1-9 of the main text, modified to take into account cell division andbinomial partitioning. Here, the cell divides when the volume reaches V m = 1 . µ m , which occurs every 20 h, as indicatedby the arrows above the graph; a cell doubling time of 20 h corresponds to an effective degradation rate of µ = 0 . / h . Theaverage volume is V = 1 µ m . (A) Time traces of the phosphorylation level and the total KaiC concentration. (B) Time tracesof the total KaiC concentration, the KaiC copy number and the volume. Please note that the time trace of the total KaiCconcentration is hardly affected when cell division happens to occur during the degradation phase, while it has a relativelylarge effect when cell division happens to occur during the KaiC production phase; this is because the removal of items fromthe queue associated with the KaiC synthesis reaction effectively lowers the synthesis rate; indeed this explains the change inslope in the KaiC concentration during the production phase. are much less affected, with the correlation number of cycles being n / = 88. Indeed, while this model combininga TTC with a PPC is fairly robust, an oscillator with exponential volume growth and binomial partitioning builtupon a TTC alone, is not stable. This supports our statement in the main text that a PPC can strongly enhance therobustness of a TTC. In a future publication, we will systematically study the effect of bursts in gene expression andbinomial partitioning. VI. ALTERNATIVE TTC MODELS
We not only studied the TTC model discussed in the main text, but also the simplest possible TTC model, namelyone in which KaiC directly represses its own synthesis (with a delay): ∅ βK / ( K +[C] ) = ⇒ τ ; σ τ C , C µ → ∅ . (S24)The result of this model is shown in Fig. S6. It is seen that this TTC model behaves very similarly to that of the maintext (compare Fig. 5B): it shows robust oscillations only for high degradation rates. Also including enzyme-mediated,zero-order protein degradation [14] does not qualitatively change this behavior. This supports our argument that theTTC has an intrinsic difficulty in generating oscillations when the doubling time is long and there is no strong activedegradation.In the simulations of Fig. 5B of the main text and Fig. S6 of the SI we vary the protein synthesis rate β c andthe protein-production delay τ such that the average KaiC concentration remains constant at the in vitro level [5, 6]and the period remains constant at 24 hours. This procedure means that as the degradation rate is increased, notonly the synthesis rate has to be increased, but also the delay. One may wonder what would happen if the delaywere kept constant and only the synthesis rate was adjusted, such that only the period remains constant at 24 hours;the average KaiC concentration is thus allowed to change. To address this, we show in Fig. S7 the results of adeterministic calculation based on the mean-field, chemical rate equations. Panel A shows as a function of the proteindegradation rate µ , the KaiC synthesis rate β that keeps the oscillation period at 24 hours, for different values of thedelay. It is seen that if the delay would be kept constant at 5-6 hours, as in the full model, the synthesis rate wouldhave to be increased dramatically when the degradation rate is increased, to keep the period at 24 hours. This wouldraise not only the average KaiC concentration, but, more importantly, also the protein synthesis rate ; this wouldincrease the energetic cost of sustaining a TTC dramatically, as shown in Fig. S7B. In the main text, we keep notonly the clock period constant, but also the average KaiC concentration. This means that when the degradation rateis increased, not only the synthesis rate, but also the delay is increased; however, even in this scenario, the cost ofsynthesizing proteins still increases markedly (Fig. S7B). -3 -2 -1 C o rr e l a t i on N u m be r o f C yc l e s Degradation Rate µ [1/h]PPC TTC
Full
FIG. S6: A comparison of the performance of the PPC-only model (Eqs. S1 - S7), the full model including both PPC andTTC (Eqs. S1 - S5 and Eqs. S13 - S15), and the simplest possible TTC model, namely one in which KaiC directly repressesits own synthesis (see Eq. S24). Compare to Fig. 5B of the main text.
VII. PERIOD AS A FUNCTION OF CELL VOLUME AND PROTEIN DEGRADATION RATE
Fig. S8 shows the period of the oscillation of the KaiC phosphorylation level in the full model of the main text(Eqs. 1–9) as a function of the cell volume and the protein degradation rate. As before, the protein synthesis rates areadjusted such that the average protein concentrations are constant and similar to those used in the in vitro experiments[5, 6]. It is seen that the dependence of the oscillation period on the cell volume and protein degradation rate is ratherweak. This is because the rhythm of the clock is dictated by the PPC, which is insensitive to the absolute rates of (A) P r odu c t i on R a t e [ µ M / h ] Degradation Rate µ [1/h] τ = 6h τ = 7h τ = 8h τ = 9h τ = 10h (B) P r odu c t i on C o s t [ µ M / h ] Degradation Rate µ [1/h] τ = 6h τ = 7h τ = 8h τ = 9h τ = 10h FIG. S7: The effect of changing the delay and the protein synthesis rate in a simple TTC based on Eq. S24. (A) The KaiCproduction rate β that keeps the oscillation period at 24 hours, as a function of the KaiC degradation rate, for different valuesof the delay τ . (B) The average synthesis rate ( i.e. the average number of KaiC molecules produced per clock period), which isa measure for the cost associated with protein production during one period, as a function of the degradation rate. The modelis a simple TTC in which KaiC directly represses its own synthesis, given by Eq. S24. The calculations were performed usingthe deterministic, mean-field chemical rate equations, because performing stochastic simulations with high production rates isprohibitively expensive. The crosses denote the points where not only the oscillation period is 24 hours, but also the averageKaiC concentration equals that used in the main text, which is the KaiC concentration used in the in vitro experiments [5, 6].The lines end for low degradation rates because the oscillations cease to exist. Panel A shows that to keep the clock period at24 hours, the synthesis rate has to increase with the degradation rate, if the delay is kept constant. This not only raises theaverage KaiC concentration, but also increases the energetic cost of sustaining the TTC, as shown in panel B. -1 Volume V [ µ m ]10 -3 -2 -1 D eg r ada t i on R a t e µ [ / h ]
20 22 24 26 28 30 P e r i od [ h ] FIG. S8: Period of the oscillation in the KaiC phosphorylation level in the full model of the main text as a function of cellvolume and protein degradation rate. When the degradation rate is varied, the protein synthesis rates are adjusted such thatthe average protein concentrations are constant and similar to those used in the in vitro experiments [5, 6]. The period isessentially independent of the volume, and exhibits only a weak dependence on the degradation rate µ . protein synthesis and decay. We note here that the period of the oscillation, as well as its amplitude, would changeif the ratio of the concentrations of the Kai proteins were changed. While the dependence of both the amplitude andthe period of the in vitro PPC on the ratio of the concentrations of the Kai proteins has been characterized in detail[5, 6], the dependence of the in vivo oscillator on their ratio has not been studied experimentally.
VIII. RHYTHMS OF
KAIBC
EXPRESSION WHEN
KAIA
IS OVEREXPRESSED
Kitayama et al. have shown that kaiBC expression oscillates with a circadian period in the presence of an excess ofKaiA, although it is not clear whether these oscillations are sustained or damped [15]. Our PPC-TTC model of themain text, which is described by Eqs. S1—S5 and Eqs. S13—S15 above, generates oscillations in kaiBC expressionwith a period of 24 hours when kaiA is overexpressed threefold, as shown in Fig. S9A. This figure shows that thephosphorylation level also exhibits weak oscillations, which are not seen experimentally; this is due to the fact thatour PPC model neglects phosphorylation of inactive KaiC, which is known to occur at high KaiA concentrations [7].To rectify this, we have extended our model to include KaiA-stimulated phosphorylation of inactive KaiC, using thesame reactions as those used for active KaiC; specifically, we add to the reactions of Eqs. S1—S5 and Eqs. S13—S15the following reactions: A y B x e C i + A ˜ k Af i ⇄ ˜ k Ab i A y B x e C i A ˜ k pf → A y B x e C i +1 + A , (S25)for each phosphorylation level i ; the rate constants equal those of the corresponding reactions of active KaiC, exceptthat the KaiA-KaiC association rate is reduced by a factor 100. All other rate constants are as in Table S1. We alsoinclude auto-activation of RpaA via e R k ma → R , (S26)with k ma = 25hr − . Auto-activation of RpaA becomes necessary because the concentrations of the KaiC phosphoformsthat activate RpaA become very low when KaiA is in excess. (The freshly injected KaiC hexamers don’t make itto the bottom of the phosphorylation cycle because of the excess KaiA.) This model not only matches the in vitro observation that, when an excess of KaiA is added during the dephosphorylation phase, the phosphorylation level ofKaiC rises immediately [7], but also reproduces the in vivo oscillations of the total amount of KaiC when KaiA isoverexpressed [15], as shown in Fig. S9B.0 (A) K a i C c on c en t r a t i on [ µ M ] Total KaiC0 24 48 72 96 120 144 168 192 216 t P ho s pho r y l a t i on r a t i o Phosphorylation ratio (B) K a i C c on c en t r a t i on [ µ M ] Total KaiC0 24 48 72 96 120 144 168 192 216 t P ho s pho r y l a t i on r a t i o Phosphorylation ratio
FIG. S9: The PPC-TTC model predicts oscillations in kaiBC expression when kaiA is overexpressed. (A) The original PPC-TTC model of Eqs. S1—S5 and Eqs. S13—S15 (B) The modified PPC-TTC model, given by the original model of Eqs. S1—S5and Eqs. S13—S15 but extended to include KaiA-stimulated phosphorylation of inactive KaiC (Eq. S25) and auto-activationof RpaA (Eq. S26). Note that the total amount of KaiC exhibits pronounced, albeit somewhat arhythmic, oscillations, whilethe phosphorylation level is high and fairly constant. For both panels, the protein degradation rate is µ = 0 . − and thevolume is V = 2 . µ m ; kaiA is overexpressed by a factor of 3. IX. MEASURING THE ROBUSTNESS
In this section, we discuss how we calculate the correlation number of cycles n / for our various models. We beginwith some theoretical background: Consider a phase variable ϕ ( t ) that increases with an average frequency ω and isalso subject to noise. Its time evolution can be written as dϕ ( t ) dt = ω + ξ ( t ) with h ξ ( t ) ξ ( t ′ ) i = σ δ ( t − t ′ ) . (S27)Here ξ ( t ) is Gaussian white noise of strength σ , and h◦i indicates averaging over different realizations of the noise.Integrating the equation with the initial condition ϕ (0) = 0 yields ϕ ( t ) = ωt + W ( t ) , (S28)where W ( t ) is a Gaussian random variable with mean zero and variance (cid:10) ( ϕ ( t ) − ωt ) (cid:11) = σ t. (S29)From this, we can construct an oscillating signal x ( t ) = x + a · sin ϕ ( t ) (S30)with mean x and amplitude a . The autocorrelation function of (S30) is then C ( t ′ ) = h δx ( t ) δx ( t + t ′ ) ih δx ( t ) i = e − σ | t ′ | · cos( ωt ′ ) , (S31)where δx ( t ) ≡ x ( t ) − x is the deviation of the signal from its average value x . We thus expect that the correlationfunction is a sinusoid modulated by a single exponential decay. A. Incorporating Amplitude Noise
Because the phase is the only “soft” direction of the dynamical system, one generically expects any noisy oscillatorat long times to act similarly to the simple model of a phase oscillator just described. Nonetheless, a more realisticdescription would also include a fluctuating amplitude. We incorporate this behavior phenomenologically by includinga time dependence in the amplitude: a → a ( t ) = a + ξ ( a ) ( t ) ⇒ x ( t ) = x + (cid:0) a + ξ ( a ) ( t ) (cid:1) · sin (cid:0) ωt + W ( t ) (cid:1) (S32)1where ξ ( a )( t ) denotes a Gaussian white noise process of strength σ a and therefore neglects correlations in the amplitudefluctuations. We can again calculate the correlation function analytically. By definition, we have C (0) = 1 for t = 0,but because a ( t ) now contains a white noise term, C ( t ) jumps discontinuously to a smaller value for any t >
0. Onefinds C ( t ′ ) = h δx ( t ) δx ( t + t ′ ) ih δx ( t ) i = a a + σ a | {z } = ν e − σ | t ′ | · cos( ωt ′ ) , t ′ > . (S33)Considering the more natural case of a finite correlation time in the amplitude noise, the picture will only changeslightly: instead of jumping from 1 to ν at τ = 0, the envelope will undergo a smooth transition for this envelopeinvolving two time scales: a short time scale of order the correlation time of the amplitude fluctuations and a muchlonger time scale associated with the phase diffusion. In practice, we found that including amplitude fluctuations didnot significantly change our estimate for n / . B. Computing the Correlation Number of Cycles
To calculate the correlation number of cycles n / , we begin by using our simulation results to estimate the corre-lation function C ( t ). After an initial equilibration phase of 500 hr, we do simulations for 50,000 hr. From these, weextract N = 500 ,
000 values x i for the time trace x ( t ) at equidistant points in time, which we can then use to calculatethe correlation function at times separated by an interval ∆ t = 0 . x i = x ( t + i · ∆ t ) , C ( i · ∆ t ) = 1( N − i ) · h δx i N − i X j =0 δx j δx j + i . (S34)Here x ( t ) is either the phosphorylation level p ( t ) or the total concentration of KaiC hexamers. p ( t ) is used for allcases except for models where there is only a TTC and a phosphorylation level cannot be defined.Once C ( t ) has been computed, it can be fitted to the form (S33) to determine the free parameters ν , σ , and ω . Inpractice we perform the fit using gnuplot , which in turn implements a Marquardt-Levenberg algorithm.Finally, it remains to translate the fitted parameter values into an estimate of n / . Eguchi et al. define n / asthe number of cycles after which the standard deviation of the phase is π [16]. Thus, we have q(cid:10) [ ϕ ( L · n / ) − ωt ] (cid:11) = π ⇒ n / = π Lσ , (S35)where L = 2 π/ω denotes the period and we have used Eq. (S29) to solve the left equation for n / . The value of n / is then obtained by substituting the fitted values for σ and ω . Using this method we can reliably measure n / from1 to 1,000. The upper bound is given by fact that we only compute the correlation function C ( t ) up to t = 3 ,
000 hand therefore correlation functions which decay on much longer timescales are difficult to detect.Concerning the error bar on the computed n / , it should be noted that there are two distinct sources of error: onedue to statistical fluctuations, and one due to systematic errors, e.g. that the correlation function cannot be fitted toEq. S33. To estimate the former, for certain parameter values we repeated the procedure described above 32 times, i.e. we performed 32 independent simulations and computed the mean and the standard deviation of the set of 32independently estimated values for n / . For V = 1 µ m and µ = 0 . − , this gave n / = 168 ±
31 (mean ± std.deviation) for the combined PPC-TTC model. To address the second type of error, we also computed n / using twodifferent methods. One is the method of Eguchi et al. , which is based on examining the times at which the oscillationreaches its maximum in each cycle, and thus does not assume any particular form like Eq. (S30) for the oscillatingsignal [16]. The other method involves computing the width of the dominant peak in the power spectrum of thetime traces. The three methods gave similar error bars and values for n / that agree within the error bar. Whileall methods gave the same result, we found the method based on the correlation function (Eq. S33) more robust fornoisy oscillations at low cell volumes. [1] van Zon JS, Lubensky DK, Altena PRH, ten Wolde PR (2007) An allosteric model of circadian KaiC phosphorylation. Proc Natl Acad Sci USA [2] Imai K, Nishiwaki T, Kondo T, Iwasaki H (2004) Circadian rhythms in the synthesis and degradation of a master clockprotein KaiC in cyanobacteria. J Biol Chem
Simulating, Analyzing, and Animating Dynamical Systems: A Guide to Xppaut for Researchers andStudents (SIAM, Philadelphia).[4] Gillespie DT (1977) Exact stochastic simulation of coupled chemical reactions.
J. Phys. Chem.
Mol Cell
FEBS Lett .[7] Rust MJ, Markson JS, Lane WS, Fisher DS, O’Shea EK (2007) Ordered phosphorylation governs oscillation of a three-protein circadian clock.
Science
EMBO J
Nature
Nat Gen
Cell
Biophys J
Proc.Natl. Acad. Sci. USA
Phys Rev Lett
Genes Dev