Circadian (∼24 h) timekeeping is essential for the lives of many organisms. To understand the biochemical mechanisms of this timekeeping, we have developed a detailed mathematical model of the mammalian circadian clock. Our model can accurately predict diverse experimental data including the phenotypes of mutations or knockdown of clock genes as well as the time courses and relative expression of clock transcripts and proteins. Using this model, we show how a universal motif of circadian timekeeping, where repressors tightly bind activators rather than directly binding to DNA, can generate oscillations when activators and repressors are in stoichiometric balance. Furthermore, we find that an additional slow negative feedback loop preserves this stoichiometric balance and maintains timekeeping with a fixed period. The role of this mechanism in generating robust rhythms is validated by analysis of a simple and general model and a previous model of the Drosophila circadian clock. We propose a double‐negative feedback loop design for biological clocks whose period needs to be tightly regulated even with large changes in gene dosage.
An accurate mathematical model of the mammalian circadian clock provides novel insights into the mechanisms that generate 24‐h rhythms. A double‐negative feedback loop design is proposed for biological clocks whose period needs to be tightly regulated.
A 1–1 stoichiometric balance and tight binding between activators (PER–CRY) and repressors (BMAL1–CLOCK/NPAS2) is required for sustained rhythmicity.
Stoichiometry is balanced by an additional negative feedback loop consisting of a stable activator.
Our detailed model can explain more experimental data than previous models.
Mathematical analysis of a simple model supports our claims.
Circadian (∼24 h) clocks time many physiological and metabolic processes. When these clocks were first discovered, three basic properties were identified (Dunlap et al, 2004). (1) Rhythms need to be autonomous. (2) Rhythms need to be capable of adjusting in response to external signals. (3) Rhythms need to persist over a wide range of temperatures. More recently, the biochemical mechanisms of circadian timekeeping have been identified (Ko and Takahashi, 2006). In particular, interlocked transcription–translation feedback loops (TTFLs) have been discovered as the basic mechanism of rhythm generation in many organisms (Novak and Tyson, 2008). With this discovery, recent experimentation has identified another property of circadian rhythms in higher organisms. Circadian rhythms persist with a 24‐h period even in the presence of large changes in the expression of the components of these TTFLs (Ko and Takahashi, 2006; Dibner et al, 2009). While mechanisms for rhythm generation with a flexible period have been identified (Stricker et al, 2008; Tsai et al, 2008; Tigges et al, 2009), mechanisms for this robustness of period to gene dosage remain unexplained, even by mathematical models (Dibner et al, 2009).
Two interlocked negative feedback loops have been identified in the TTFL networks generating circadian rhythms in higher organisms (Figure 1) (Blau and Young, 1999; Glossop et al, 1999; Benito et al, 2007; Liu et al, 2008). A ‘core’ negative feedback loop consists of repressors (PERIOD and TIMELESS in Drosophila or PERIOD1–3 and CRYPTOCHROME1–2 in mammals), which inactivate activators (CYCLE and CLOCK in Drosophila and BMAL1–2 and CLOCK in mammals) of their own transcription. An additional negative feedback loop controls the expression of the activators, which inactivate their own transcription through Vrille (Drosophila) or the Rev‐erbs genes (Mammals) (Blau and Young, 1999; Preitner et al, 2002). While other feedback loops have also been identified, these two negative feedback loops seem to predominate (Blau and Young, 1999; Glossop et al, 1999; Benito et al, 2007; Liu et al, 2008; Bugge et al, 2012; Cho et al, 2012).
Near 24‐h oscillations persist even when the components of the TTFLs of the circadian clock are over or under expressed. Heterozygous mutations of clock genes never abolish rhythmicity, and their period phenotypes are either indistinguishable from the wild‐type (WT) phenotypes or much smaller than mutations that affect post‐translational modifications (Baggs et al, 2009; Etchegaray et al, 2009; Lee et al, 2009). Abolishing rhythmicity through single gene knockout is surprisingly difficult (Baggs et al, 2009; Ko et al, 2010). Moreover, the mammalian circadian clock is also resistant to global changes in transcription rates (Dibner et al, 2009). These results all suggest that gene dosage may not be important for circadian timekeeping in higher organisms.
Gene dosage, however, is not completely unimportant for timekeeping. Knockdown of clock genes causes increased expression in similar components through paralog compensation, which may help restore gene dosage and indicates that gene dosage needs to be tightly regulated (Baggs et al, 2009). Population rhythmicity in mouse embryonic fibroblasts shows much lower amplitude than in liver, which might be due to the fact that the ratio of repressors to activators is significantly lower in fibroblasts than that found in liver (Lee et al, 2001, 2011). A 1–1 stoichiometric binding occurs between the activators and repressors driving rhythms in Drosophila (Menet et al, 2010), although not in Neurospora (He et al, 2005; Huang et al, 2007).
Here, we propose a mechanistic explanation for the robustness to gene dosage in the circadian clock of higher organisms through mathematical modeling. We develop the most detailed mathematical model of the mammalian circadian clock available, which should be useful in many future studies. Our model reproduces a surprising amount of experimental data on the mammalian circadian clock including the time courses and relative concentrations of key transcripts and proteins, the effects of mutations of key clock genes, and the effects of changes in gene dosage. With this model, we show that proper stoichiometric balance between activators (BMAL–CLOCK/NPAS2) and repressors (PER1–2/CRY1–2) is key to sustained oscillations. Furthermore, we find that an additional slow negative feedback loop, in which activators indirectly inactivate themselves, improves the regulation of the stoichiometric balance and sustains oscillations with a nearly constant period over a large change in gene expression level. Tight binding between activators and repressors is also predicted to be crucial for rhythm generation. These mechanisms are also validated by mathematical analysis of a simplified mathematical model of the mammalian circadian clock, and simulations of a previously published Drosophila model. We here propose a novel design for biological oscillators where maintaining period is crucial: a core negative feedback loop with repression by protein sequestration, with an additional negative feedback loop, which controls a relatively stable activator.
Mathematical modeling of the mammalian circadian clock
We develop a new mathematical model of the intracellular mammalian circadian clock. This model contains key genes, mRNAs and proteins (PER1, PER2, CRY1, CRY2, BMAL1/2, NPAS2, CLOCK, CKIε/δ, GSK3β, Rev‐erbα/β) that have been found to be central to mammalian circadian timekeeping (Figure 1A). While greatly expanded, the model is largely based on our previous model, which has made surprising predictions about mammalian timekeeping that have been subsequently verified experimentally (Forger and Peskin, 2003; Gallego et al, 2006; Ko et al, 2010; Yamada and Forger, 2010). Modifications and extensions of the model are described in the Materials and methods, Supplementary information and Supplementary Tables 1 and 2. The parameters of the model are estimated using experimental data and a simulated annealing method (a global stochastic parameter searcher) (Gonzalez et al, 2007) (see Materials and methods, Supplementary information and Supplementary Table 3 for details). In particular, we incorporated experimentally determined rate constants (Supplementary Table 3) (Kwon et al, 2006; Siepka et al, 2007; Chen et al, 2009; Suter et al, 2011), fit the time courses of both mRNA and proteins (Figure 2A and B) (Lee et al, 2001; Reppert and Weaver, 2001; Ueda et al, 2005) and fit the relative abundance of proteins (Figure 2C) (Lee et al, 2001).
Our model accurately predicts the phenotype of known mutations of genes in the central circadian clock (suprachiasmatic nuclei, SCN) (Yoo et al, 2005; Baggs et al, 2009; Ko et al, 2010), which other models do not predict (Table I) (Forger and Peskin, 2003; Leloup and Goldbeter, 2003; Mirsky et al, 2009; Relógio et al, 2011). Interestingly, our model shows opposite phenotypes for Cry1−/− and Cry2−/− matching experimental data (Liu et al, 2007). There are two differences between CRY1 and CRY2 in our model. First, Cry1 transcription is delayed through repression by Rev‐erbα and Rev‐erbβ (Preitner et al, 2002; Liu et al, 2008; Ukai‐Tadenuma et al, 2011). Additionally, Cry1 mRNA is more stable than Cry2 mRNA and CRY1 protein is more stable than CRY2 protein (Busino et al, 2007; Siepka et al, 2007; Chen et al, 2009). Since a longer half‐life causes rhythms to be delayed, and delayed rhythms cause a longer period (Forger, 2011; Ukai‐Tadenuma et al, 2011), removing CRY1 shortens the period and removing CRY2 lengthens the period. The opposite phenotypes of Clock−/− (null mutation) and ClockΔ19/+ (dominant‐negative mutation) are also correctly simulated in the model for the first time (Vitaterna et al, 1994; Herzog et al, 1998; Debruyne et al, 2006). Moreover, our model also predicts the mutant phenotypes of isolated SCN neurons, which are different from the SCN slices (Liu et al, 2007). We note that SCN slices have significantly higher gene expression of per1 and per2 through CREB/CRE pathway than isolated SCN neurons (Yamaguchi et al, 2003). Interestingly, when we reduced per1 and per2 expression about 60% in our model, our model was able to accurately reproduce the phenotypes of isolated SCN neurons (Table II).
We also conducted a sensitivity analysis to look at what parameters determine the period of our model. Four of the top five high parameters, in our sensitivity analysis, were also in the top five found in a previous sensitivity analysis with the original Forger and Peskin model and which was used to conclude that PER2 plays a dominant role in period determination (Wilkins et al, 2007) (see Supplementary Figure 1).
Proper stoichiometric balance between activators and repressors is crucial to sustained rhythms
Since our mathematical model can accurately predict the phenotype of known mutations of the mammalian circadian clock, we next looked for a mechanism that could explain why some phenotypes were rhythmic, while others were not. We found that stoichiometry plays a key role in determining which mutations showed rhythmic phenotypes. Here, we define stoichiometry as the average ratio between the concentration of repressors (all forms of PER and CRY in the nucleus) to that of activators (all forms of BMAL–CLOCK/NPAS2 in the nucleus) over a period. Moreover, we specifically refer to repressors and activators of E/E’‐boxes when discussing stoichiometry. We found that mutations that caused the stoichiometry to be too high or too low, yielded arrhythmic phenotypes (Figure 3A). So long as the mutations allowed the stoichiometry to be around a 1–1 ratio, relatively high amplitude oscillations were seen. Thus, we predict that stoichiometry provides a unifying principle to determine the rhythmicity of mutations of the mammalian circadian clock. To further test this principle, we constitutively expressed either the Per2 gene (the dominant repressor gene) or the Bmal and Clock genes (the dominant activator genes) at different levels. Interestingly, within a range centered near a 1–1 stoichiometry, the model shows sustained oscillations with high amplitude (Figure 3B). However, if the stoichiometry was too high or too low, rhythms are dampened or completely absent (Figure 3B). This matches a recent experimental study showing that the amplitude and sustainability of population rhythms increase when the level of PER–CRY is increased closer to that of BMAL1–CLOCK in mouse fibroblasts (Lee et al, 2011).
We defined the stoichiometry as the average ratio between the total concentrations of repressors to that of activators over a period. However, recent work has shown that CRY1 has stronger repressor activity than CRY2. The underlying biochemical mechanisms for this result have not been fully identified (Khan et al, 2012). If the difference is due to a different post‐translational mechanism (e.g., binding between PER and CRY, which could affect the repressor concentration in the nucleus), the current definition of stoichiometry can be kept. Otherwise, a more sophisticated definition of stoichiometry may be needed (e.g., one that gives more weight to concentration of CRY1 than that of CRY2).
How stoichiometry generates rhythms
To test the role of stoichiometry in sustaining oscillations, we developed a simple model by modifying the well‐studied Goodwin model (Goodwin, 1965) to include an activator (A), which becomes inactive when bound by a repressor (P) (Figure 3C). Transcription is proportional to the fraction of free activator that is not bound by the repressor, f(P, A, Kd) (Buchler and Cross, 2009), matching experimental data from the mammalian circadian clock (Supplementary Figure 2) (Froy et al, 2002). mRNA (M) is translated to a repressor protein (Pc). The protein enters the nucleus (P) and binds and inhibits the activator (A). This generates a single‐negative feedback loop (SNF) since the activator is constitutively expressed. The model is similar to a previously published mathematical model (Francois and Hakim, 2005); however, we allow for both association and dissociation of the activator and repressor (through a defined Kd), which turns out to be crucial for understanding the effects of stoichiometry. By nondimensionalization and setting the clearance rates of all species to be equal (to increase the chance of oscillations, see Forger, 2011), only two parameters remain: the activator concentration (A) and the dissociation constant (Kd) (see Supplementary information).
When we changed the activator concentration, which changed the stoichiometry (average ratio between the level of repressor (P) to the level of activator (A)), sustained oscillations were only seen at around a 1–1 stoichiometry similar to our detailed model (Figure 3D). As the other parameter (Kd) decreased (indicating tight binding), the range of stoichiometry that permitted oscillation increased (Figure 3D). Interestingly, if the binding was too weak, the rhythms did not occur. The tight binding between activators and repressors is also found in the detailed model, and in the mammalian circadian clock (Lee et al, 2001; Froy et al, 2002; Sato et al, 2006). This indicates that the sustained rhythms require tight binding as well as balanced stoichiometry in the circadian clock.
Many previous studies have argued that ultrasensitive responses (e.g., a large change in transcription rate for a small change in repressor or activator concentration) can cause oscillations in feedback loops (Kim and Ferrell, 2007; Buchler and Louis, 2008; Novak and Tyson, 2008; Forger, 2011). A previous study showed that an ultrasensitive response can be generated by tight binding of activators and repressors in a synthetic system (Buchler and Cross, 2009). Taken together, this provides a potential mechanism of rhythm generation. That is, when the total concentration of repressor is higher than that of activators, the repressor sequesters and buffers activator and inhibits transcription completely (Buchler and Louis, 2008). As the repressor is depleted, the excess free activators are no longer sequestered by repressors and are free to turn on the transcription. At this threshold, transcription of repressor shows an ultrasensitive response to the concentration of repressor or activator. Ultrasensitive responses amplify rhythms and prevent rhythms from dampening (Forger, 2011). In both our simple and our detailed model, we found ultrasensitive responses around a 1–1 stoichiometry (Supplementary Figure 3A). When the stoichiometry was not around 1–1, an ultrasensitive response was not seen, and both models did not show sustained rhythms.
Over the course of a day, as levels of repressor and activator change, the stoichiometry and also sensitivity change as well. We found that the 1–1 average stoichiometry is required to generate the ultrasensitive response, which causes rhythms through mathematical analysis, confirming our simulation results (Figure 3D). That is, via both local and global stability analysis, we derived an approximate range of the stoichiometries (〈S〉) that permit oscillations
(see Supplementary information). In agreement with our simulations shown in Figure 3D, this mathematical analysis also suggests that (1) oscillations are seen around a 1–1 stoichiometry; (2) the stoichiometry needs to be >8/9 for sustained rhythmicity; (3) as the binding between activators and repressors becomes tighter, the upper bound on stoichiometry increases; (4) if the binding is too weak (e.g., Kd=10−3), sustained oscillations do not occur.
An additional negative feedback loop improves the regulation of stoichiometric balance
If stoichiometry is key to sustained oscillation, are there mechanisms within circadian clocks that keep the stoichiometry of components balanced? Does the additional negative feedback loop of the negative–negative feedback loop (NNF) structure, found in circadian clocks, help balance stoichiometry? To test this structure, we added an additional negative feedback loop into our simple model (Figure 4A). Previously, other studies suggested that an additional positive, rather than negative, feedback loop could sustain intracellular clocks (Barkai and Leibler, 2000; Stricker et al, 2008; Tsai et al, 2008; Tigges et al, 2009). We tested these structures by including an additional protein R (Rev‐ERBs or RORs in the mammalian circadian clocks) that is transcribed in a similar way to P. R then represses (as in the Rev‐erbs) or promotes (as in the Rors) the production of A in the negative–negative feedback loop (NNF) or the positive–negative feedback loop (PNF) structure, respectively (Figure 4A).
We studied how the SNF, NNF and PNF structures effectively maintain the stoichiometric balance when model parameters (e.g., transcription rate) are changed. With both simulation and steady‐state analysis, we found that the NNF structure is best at keeping stoichiometry balanced while the PNF structure is worst at keeping stoichiometry balanced, regardless which parameters are perturbed (see Supplementary information, Figure 4B and Supplementary Figure 4A–C). Moreover, our detailed model, which also follows the NNF structure, also carefully balanced the stoichiometry by controlling the expression of repressors and activators. Knockdown of the repressor Cry1 leads to higher expression of the repressors, which are controlled by E‐boxes, and lower expression of the activators, which are controlled by a ROREs (Figure 4C). Opposite effects are seen when the activator CLOCK is removed (Figure 4C). This active control of repressors and/or activators via the NNF structure regulates the stoichiometric balance tightly (Supplementary Figure 4D) and matches experimental data on gene dosage (Baggs et al, 2009). Moreover, the detailed model (with the NNF structure) also correctly predicts the change of clock gene expression after the removal of the additional negative feedback loop (Rev‐erbα,β−/−) (Figure 4D) (Liu et al, 2008; Bugge et al, 2012; Cho et al, 2012). In particular, knockout of the Rev‐erbα,β decreases PER expression, but increase CRY1 expression. For our nominal set of parameters, oscillations are still possible when this additional negative feedback is removed. However, for other sets of parameters, where stoichiometry is not as well balanced, removal of this additional negative feedback stops rhythmicity (see below). This could explain the phenotype of the Rev‐erbα,β−/−, which show some indications of rhythmicity (Bugge et al, 2012; Cho et al, 2012). Our model predicts that rhythm generation remains in cell types that have a near balanced stoichiometry, and a lack of rhythms in cell types without a balanced stoichiometry.
A slow additional negative feedback loop improves the robustness of rhythms
Our central hypothesis is that, as stoichiometry is more tightly regulated, oscillations will occur over a wider range of parameters. To confirm this, we varied the transcription rate of the activator (or activator concentration in the SNF) and the transcription rate of the repressor to determine which sets of parameters yielded oscillations. While the SNF, NNF and PNF structures have almost the same behavior with their nominal parameters (mean stoichiometry, amplitude and period, see Supplementary Figure 5A), the NNF structure oscillated over the widest range of parameters and the PNF oscillated over the narrowest range of parameters in the simple model (Figure 5A; Supplementary Figure 5C). Interestingly, as the activator becomes more stable (i.e., the additional negative feedback becomes slower), the NNF structure allows sustained oscillations over a wider range of parameters (Supplementary Figure 5D). Indeed, the clearance rate of the activators is significantly slower than other circadian clock components (Supplementary Table 4) (Kwon et al, 2006).
We also checked the role of the NNF structure in our detailed mammalian clock model. We modified the NNF structure of the detailed model to that of an SNF by fixing the activator (BMAL, CLOCK and NPAS2) concentration to the average value found in their WT simulations. We also constructed the PNF structure by converting the repressor (REV‐ERBs) to an activator (e.g., the RORs) in the NNF structure. This did not significantly change the rhythms in the core feedback loop (Supplementary Figure 5B), matching previous studies that showed that the loss or change in rhythms in the activators had little effect on the circadian rhythms (Liu et al, 2008). It is tempting to conclude that the additional feedback loops controlling activators are not important in the circadian clocks. However, when we changed the transcription rate of the repressor (Per) and activator (Bmal, Clock and Npas2), the original model (with an NNF structure) had the widest range of parameters where oscillations occur while the PNF structure had the narrowest range of parameters (Figure 5B). Interestingly, experiments have shown that REV‐ERBs play a more dominant role than the RORs, indicating that our proposed mechanism may play an important role in in vivo timekeeping (Liu et al, 2008). Thus, the choice of the additional feedback greatly affected the range of parameters where oscillations are seen.
We also examined the role of the additional negative feedback loop in a mathematical model of the Drosophila circadian clock (Smolen et al, 2002). The original study that developed the model concluded that the NNF and SNF structures were equally likely to show oscillations. However, their study only changed transcription rates by 20%. With a larger perturbation of parameters, we found that the additional negative feedback loop significantly extends the range of parameters that yield oscillations (Figure 5C).
A network design for cellular timekeeping where maintaining a fixed period is crucial
The PNF structure can create a robust biological oscillator that has a tunable period when the additional positive feedback loop is fast (i.e., the activator degrades quickly) (Stricker et al, 2008; Tsai et al, 2008; Tigges et al, 2009) (Figure 6A). Consistent with these findings, our simple model with the PNF structure has a tunable period for changes in gene expression levels (Figure 6C). However, the simple model with the NNF structure has a nearly constant period in the presence of large changes in gene expression levels (Figure 6B and D). Furthermore, this NNF structure becomes more robust as the additional negative feedback loop slows (i.e., the activator degrades more slowly) (Supplementary Figure 5D) in contrast to the fast positive feedback of the tunable clocks (Stricker et al, 2008; Tsai et al, 2008; Tigges et al, 2009). Consequently, our results propose two different designs for robust biological oscillators. The NNF structure (Figure 6B) is suitable for biological clocks in which the maintenance of a fixed period is crucial (e.g., circadian clocks). The PNF structure (Figure 6A) is suitable for the biological oscillators that need to tune their period (e.g., cell cycle or pacemaker in the sino‐atrial node) (Tsai et al, 2008). This is also supported by mathematical analysis of the simple model (for more details, see Supplementary information).
Our work identifies several key mechanisms that allow 24‐h rhythms in the circadian clocks of higher organisms: (1) Proper stoichiometric balance between the activators and the repressors, (2) tight binding between activators and repressors, (3) the NNF structure and (4) longer half‐life of activators than repressors. These mechanisms synergistically generate rhythms with periods robust to gene dosages (Figure 6D). The range of the stoichiometry where the rhythms occur widens as binding between activators and repressors tightens (Figure 3D). Moreover, the NNF structure regulates the expression of activators as well as repressors to balance stoichiometry (Figure 4B and C). For instance, increased stoichiometry (elevated repressor concentrations) strengthens the repression in the core negative feedback loop and reduces the expression of the repressors (e.g., Pers and Crys) and Rev‐erbs. The decreased expression of Rev‐erbs weakens the additional negative feedback and increases the expression of activators (Bmal1 and Npas2), which lowers the stoichiometry (Supplementary Figure 4D). When this is done on a slower timescale, so that the basics of the 24‐h timekeeping are unaffected, the robustness of the rhythms is enhanced (Supplementary Figure 5D).
Relation to previous experimental data
Many experimental observations could be interpreted as mechanisms by which the mammalian circadian clock balances stoichiometry. When the repressor (CRY) is overexpressed or the repressor (PER) is removed, the activator (BMAL1) concentration is found to increase or decrease, respectively (Shearman et al, 2000; Fan et al, 2007). When a repressor's expression is reduced, the expression of other repressors is increased and the expression of activators is decreased (Baggs et al, 2009). Knockdown of activators yields opposite effects (Baggs et al, 2009). Both our detailed and simplified NNF models confirm these results (Figure 4B and C; Supplementary Figure 4D). Additionally, the rhythms of the mammalian circadian clock persist even after the transcription of all clock genes are reduced significantly (Dibner et al, 2009). In agreement with these data, both the detailed and the simple model oscillate after significant reduction of the transcription rates of both activators and repressors because their stoichiometry is maintained (Supplementary Figure 5E). Our study also suggests an underlying mechanism (ultrasensitive response) for a previous experimental observation showing that the robustness of circadian rhythms is enhanced by making the level of PER–CRY closer to that of CLOCK–BMAL1 in mouse fibroblasts (Supplementary information; Supplementary Figure 3A) (Lee et al, 2011).
Experimental data also support the role of the slow additional negative feedback loop in regulating circadian timekeeping in higher organisms. The time course of the activator (BMAL1 in mammal or CLK in Drosophila) seems to be controlled mainly by the additional negative feedback loop (Rev‐Erbs or Vrille) (Blau and Young, 1999; Benito et al, 2007; Liu et al, 2008). The elimination of additional positive feedback loop has little effect on circadian clocks in contrast to other cellular clocks based on the PNF structure (Benito et al, 2007; Kim and Ferrell, 2007; Liu et al, 2008; Tsai et al, 2008). Furthermore, a key step, the clearance rate of the activators, which governs the timescale of the additional feedback loop, is significantly slower than other circadian clock components (Supplementary Table 4) (Kwon et al, 2006). Removing the slow additional negative feedback loops in the mammalian clock (Rev‐erbα−/−) yields timekeeping where the period is not as well‐maintained (Preitner et al, 2002). Moreover, recent studies have confirmed a pivotal role for the additional negative feedback loop for regulating the circadian rhythms via double knockout of Rev‐erbα and Rev‐erbβ (Bugge et al, 2012; Cho et al, 2012). Thus, our proposed mechanism of robust circadian timekeeping matches known data on the mammalian circadian clock. Further comparison with known experimental data is shown in Supplementary Figure 7.
Relation to previous modeling work
Our study is the first circadian modeling study that shows the importance of a balanced stoichiometry in rhythm generation. Our results for the SNF structure match a previous model based on the protein sequestration (Francois and Hakim, 2005), which focuses on other mechanisms, for example, slow RNA dynamics, that do not play a role in circadian clocks. We have identified a basic mechanism of tight binding and protein sequestration for generating high sensitivity, similar to what has been proposed in the cell cycle and synthetic studies (Buchler and Cross, 2009), as the key rhythm generating mechanism in our model. Previous circadian clock models do not use this mechanism, and a careful justification, based on experimental data from higher organisms, of the mechanisms for generating high sensitivity and, consequently, oscillations, in these models has yet to be performed (Yoo et al, 2005). In fact, several of these mechanisms have been called into question (Forger and Peskin, 2003).
Previous models have used different mechanisms for rhythm generation (e.g., high‐Hill coefficients) and have proposed different roles for the additional negative feedback loop. They have proposed that the additional negative feedback loop is capable of independent oscillations, even when the core negative feedback loop was removed (Leloup and Goldbeter, 2004; Relógio et al, 2011). However, despite much experimental study, no oscillations have yet been found from this additional feedback loop in isolation (Sato et al, 2006) and the known phenotypes of knockout of genes in this additional feedback loop had not been correctly predicted (Preitner et al, 2002; Relógio et al, 2011). Moreover, other previous studies argued that the additional negative feedback loop is not important (Becker‐Weimann et al, 2004), which does not match with recent experimental data on the mammalian circadian clock (Bugge et al, 2012; Cho et al, 2012). We claim that the additional negative feedback loop is not an independent oscillator, nor ancillary, but acts to regulate stoichiometry.
Interestingly, the predictions of previous modeling studies (Griffith, 1968; Becker‐Weimann et al, 2004) match experimental data from the Neurospora circadian clock, in which a 1–1 stoichiometry is not important and the additional negative feedback loop seems to not play an important role (Baker et al, 2012). Our predictions match experimental data from circadian clocks in higher organisms (Supplementary Figures 7 and 8).
Proposed experiments based on model predictions
Our most important prediction may be the following: when the stoichiometry between activators and repressors is within a fixed range, oscillations are sustained, and outside this range oscillations are damped (Figure 3). This can be tested by measuring the relative concentration of activators and repressors in many tissues and in the presence of several possible mutations that lead to damped or sustained rhythms. This has been done in WT fibroblasts and liver (Lee et al, 2001, 2011), but has not been done in other tissues or mutants. Moreover, we note that these previous experiments were done in population cell assays, whereas single‐cell measurements may be needed to determine whether damped oscillations are the result of damped rhythms in single cell, or greater population desynchrony (Welsh et al, 2004; Leise et al, 2012).
The behavior of isolated SCN neurons is similar to fibroblasts in that mutations of circadian genes can easily lead to arrhythmicity (Liu et al, 2007). We note that intercellular coupling in the SCN not only synchronizes SCN neurons, but also increases transcription of per1 and per2 (Yamaguchi et al, 2003), which may balance stoichiometry and help sustain rhythms when repressors are effectively removed (Tables I and II). Thus, we predict that increasing transcription of per1 and/or per2 could enhance rhythmicity in isolated SCN neurons similar to what is seen in fibroblasts (Lee et al, 2001). Moreover, our model predicts that cells with low stoichiometry (e.g., isolated SCN neurons) shows larger phase shifts in response to light than cells with 1–1 stoichiometry (e.g., SCN slices) (data not shown). It would be interesting future work to see whether different cell types have different PRCs depending on their stoichiometry.
We also predict that tight binding between activators and repressors is required for rhythmicity (Figure 3D). Several studies have identified binding sites for PER and CRY on BMAL1 and CLOCK (Sato et al, 2006; Langmesser et al, 2008; Ye et al, 2011). Point mutations in binding sites can generate different binding affinities between PER–CRY and BMAL1–CLOCK. Comparing the experimentally measured binding affinities of these mutants, with the resultant rhythms, or lack thereof, would directly test this prediction.
Loss of the additional negative feedback loop (e.g., in the Rev‐erbs−/−, constitutive expression of Rev‐erbs or constitutive expression of BMAL) is predicted to cause the intracellular circadian clock to oscillate over a much narrower range of conditions (Figure 5). It would be interesting to test whether these cells would have less temperature compensation or would lose rhythms more easily when other genes are knocked out (e.g., Cry2−/−, Per1−/−). Moreover, we predict that in the Rev‐erbs−/−, rhythms persist in cell types with a balanced stoichiometry, but not in poorly balanced cells (Figure 5). It would be interesting future work to investigate whether SCN and peripheral clocks have different phenotypes of Rev‐erbs−/− depending on their stoichiometry. We also predict that Rev‐erbs−/− cells show a wider period distribution than WT (Figure 6).
Our modeling and analysis also predict that relatively stable activators (e.g., BMAL1 and CLOCK) in the additional negative feedback loop allow rhythmicity over a wide range of conditions (Supplementary Figure 5D). These activators can be destabilized with point mutations (Sahar et al, 2010). Simply destabilizing the activators might lead to lower activator concentrations and unbalance stoichiometry, which is also predicted to reduce rhythmicity. However, we predict a loss of rhythmicity when these activators are destabilized, even when the overall activator concentrations are controlled for.
Perhaps the most direct way to test our model is to build the clock described in our simple NNF model using the tools of synthetic biology. Other synthetic clocks have been built, and the design we propose is not more complex than what has been previously built (Stricker et al, 2008; Tsai et al, 2008; Tigges et al, 2009). Validation could first be done in an analog electric circuit, even though this might be much less convincing. Building a synthetic clock would be of particular importance since it would be the first synthetic clock predicted to have a tightly regulated period.
Further work should explore the role of the NNF structure in the presence of molecular noise (Forger and Peskin, 2005; Li and Li, 2008). Here, we studied the role of an additional negative feedback loop controlling the activators of the circadian clocks of higher organisms. Future work could consider the functions of the additional negative feedback loops in other organisms. In particular, the plant circadian clock has a different feedback loop structure than the mammalian or Drosophila circadian clocks (Pokhilko et al, 2012). It would be interesting to see if our ideas carry over to other organisms and other cellular clocks. Furthermore, other types of feedback loops in the circadian clocks of higher organisms could be explored. Here, we found that balancing stoichiometry properly might be a universal principle of biological timekeeping. This finding not only is in agreement with experiment data from the circadian clocks in higher organisms, but even in agreement with the circadian clock in cyanobacteria as well (Rust et al, 2007). It would be interesting to test the role of stoichiometry in other cellular clocks, such as developmental clocks.
Materials and methods
Modifications and extensions of the detailed model
Detailed modeling of additional feedback loops: The new model includes the secondary feedback loops, which regulate transcription of genes with a RORE in their promoters, including Bmals, Npas2 and Cry1 (Preitner et al, 2002; Debruyne et al, 2006; Liu et al, 2008).
Updated mechanisms of BMAL–CLOCK/NPAS2 repression: Matching recent findings, we updated the mechanisms by which the repressor (PER/CRY) inhibits the activator (BMAL–CLOCK/NPAS2) (Kondratov et al, 2006; Dardente et al, 2007; Chen et al, 2009; Ye et al, 2011).
Accounting for the heterogeneity of different genes with E‐boxes: We introduced three different types of E‐boxes for Per1/Per2/Cry1, Cry2 and Rev‐erbs, matching experimental data (Ueda et al, 2005; Lee et al, 2011).
Improved description of the effect of light on the circadian clock: We included a previous model of the effect of light on the circadian clock (Kronauer et al, 1999).
Variables and equations of the detailed model
The monomer proteins considered in our model are PER1, PER2, CRY1, CRY2, BMALs, CLOCK/NPAS2, REV‐ERBs, CKI and GSK3β. Although only 10 monomers are considered in the model, they can produce many complexes depending on the state of binding, phosphorylation and subcellular locations. To describe these all complexes, 181 variables are needed (Supplementary Tables 1 and 2): 159 variables are for protein complexes, 12 variables are for mRNAs, 8 variables are indicator of the promoter activity and 2 variables are for light effect and GSK3β activity. The reactions between these variables are described by ODE systems using explicit mass kinetics as in the original model (Forger and Peskin, 2003). More details of the detailed equations are provided in the Supplementary information.
Parameter estimation of the detailed model
While the original model used 36 parameters, the new model has the 75 parameters due to the extensions and modifications of the model. Despite the increased number of parameters, we could get tighter restriction on the range of parameters with newly published data (listed below). Over these ranges, parameters are estimated by fitting to more various types of data: time courses of gene expressions and proteins, abundance of proteins and mutation phenotypes.
We choose 14 parameters (degradation rate of mRNAs and proteins) matching published experimental data. These parameter values were allowed to vary up to 50% from the experimentally determined values to account for experimental error and cellular heterogeneity.
The dissociation constant between BMALs–CLOCK and CRY is set greater than that between BMALs–CLOCK and PER (Chen et al, 2009).
The ratio between cytoplasm and nucleus volume are limited to between 1 and 3.5 (Miller et al, 1989).
The other parameters are also restricted into a biologically reasonable range (see Supplementary Table S3).
Within these restrictions, a simulated annealing method (SA, a global stochastic parameter searcher) (Gonzalez et al, 2007) was used to estimate the parameters in two steps. First, we found parameters that provides a good fit with mRNA and protein time profiles measured in mouse SCN (Reppert and Weaver, 2001; Ueda et al, 2005) and relative abundance of clock proteins measured in mouse liver (Lee et al, 2001) and fibroblast (Lee et al, 2009, 2011) (Figure 2A–C). In this fitting, we used a similar cost function to that used in estimating the parameters of the original model (Forger and Peskin, 2003)
Here, j runs through 6 mRNAs and 4 proteins. nj is the number of data points (12 for mRNA and 13 for protein). sij and eij are simulated time courses and experimentally measured time courses, respectively. sij are normalized, in the same way as was done in the experimental data (see Figure 2 for details). wij=5 when eij=1 and wij=1 otherwise, so that the cost function has more weight at the peak time than other times. pmk and pk are maximum value of protein abundance, respectively. pmk and pk are normalized, so that the maximum abundance of the CRY1 protein is 1.
After the first round of SA, we found several parameter sets qualitatively matching with experimental data on phenotypes of mutations of mice (WT, short, long and AR) (Table I). Then we used these parameter sets as initial parameter sets for another round of SA to get the final parameter set, which shows a quantitatively good fit with knockout mutation phenotype as well as time profiles (Supplementary Table 3). The cost function used for the second round is followed.
mpl and ml are simulated period and experimentally measured period of rhythmic phenotypes of mutations, respectively. man are simulated relative amplitude of arrhythmic phenotypes of mutation (e.g., Per2−/− or Bmal1−/−).’
Validation of the detailed model with experimental data
Time profile of mRNA and proteins.
As previously mentioned, we fit the model simulations with time profiles of clock mRNA and proteins in SCN to estimate the unknown parameters. We followed the same experimental procedures used to measure time profiles. The model was entrained under 12 h–12 h light/dark cycle with 500 lux light strength for 20 days. Then, the concentrations of mRNAs were measured during the following 48 h in darkness and measured time courses were compared with experimental data (Ueda et al, 2005) (Figure 2A). In the same way, the simulated protein time profiles are also fit with the data (Reppert and Weaver, 2001) (Figure 2B).
Relative abundance of proteins.
Relative abundance among core clock proteins were compared with liver (Lee et al, 2001) and fibroblast data (Lee et al, 2011) because SCN data has not yet been reported (Figure 2C).
Knockout mutation phenotype.
We also tested whether our model could predict the phenotype of mutations of clock genes (Table I). Overall, the model simulation well matches with SCN or behavioral phenotype. Homozygous and heterozygous knockouts were simulated by reducing transcription rates by 100% and 50%, respectively. To simulate the Rev‐erbα−/−, we also reduced the transcription rate of the Rev‐Erbs by 50%, which represented both Rev‐erbα and Rev‐erbβ in our model. To model the Bmal1−/−, we reduced transcriptional rate of Bmals by 95%, which accounts for the low levels of Bmal2 when compared with Bmal1 (Ko et al, 2010). For the ClockΔ19/+, the half of WT CLOCK proteins were mutated to be transcriptionally inactive, yet still competed with the remaining WT CLOCK proteins. For the CK1εtau/tau, we increase the CK1 phosphorylation rate for PER1 and PER2 by four times.
Simulation of the detailed model
All the simulations and parameter search were done with 150 × 8 Ghz CPU using MATHEMATICA 8.0 (Wolfram Research).
Model description and analysis of the simple model
The simple model is available in SBML, Mathematica, Matlab, C++ and XPPAUT format from the ModelDB (Access code: 145800) (Hines et al, 2004; Mendes et al, 2009). The detailed model is available in Mathematica, Matlab and XPPAUT format from the ModelDB (Access code: 145801). (Schmidt and Jirstrand, 2006). See Supplementary information.
We thank Choogon Lee, John Hogenesch, Saad Kahn, Jihwan Myung and Justin Blau for comments on this manuscript; and Justin Dunmyre for help on XPP‐AUTO program. This work was supported by NSF DMS‐1026317 for computing resources and by AFOSR Grant FA 9550‐11‐1‐0165.
Author contributions: JK and DF designed the research; JK performed the experiments; JK and DF analyzed the data; JK and DF did mathematical analysis and JK and DF wrote the paper.
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary Figures S1–8, Supplementary Tables S1–4 [msb201262-sup-0001.pdf]
Simple model [msb201262-sup-0002.zip]
Detailed model [msb201262-sup-0003.zip]
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2012 EMBO and Macmillan Publishers Limited