Feedback Loops in the Mammal Body Clock Form a Repressilator Circuit

To study the complex gene regulatory network of the mammalian circadian oscillator, we constructed a mathematical model with only the key components as explicit variables. For example, transcriptional profiles reveal clear redundancies in the network of core clock genes [1, 4] with RORE-binding activators (Rorα,-β,-γ) exhibiting opposite phases as the RORE-binding inhibitors (Rev-erb-α,-β). This allows to describe the regulatory actions by a single term controlled by Rev-erb-α levels, while the systems behaviour remains the same. The additional effects by Ror-genes and Rev-erb-β can be taken into account by changes of parameters describing the strength of Rev-erb-α regulation. Analogously, we combine the regulations via D-boxes into one term. The Dbp-gene represents the combined effects of the activators Dbp, Hlf and Tef and the inhibitor E4bp4. Transcriptional regulation via E-boxes is particularly complex [13]. In our model, Bmal1 quantifies the positive regulation after dimerization with Clock or Npas2, while the genes Per2 and Cry1 represent early and late E-box driven genes, respectively. The essential role of a rather late Cry1 phase has been demonstrated in detail elsewhere [14, 15].

Overall, we designed a regulatory network consisting of five variables only [8]. Fig 1A shows that even this core clock network exhibits multiple negative and positive feedback loops. Importantly, our model successfully describes published phase relations, amplitudes and waveforms of clock gene expression profiles (Fig 1B). A detailed comparison with experimentally measured profiles is described in [8].

Our gene regulatory network model contains 34 kinetic parameters. Since quantitative details of post-transcriptional processes including phosphorylations, complex formations and nuclear translocation are not known, we represent these processes by explicit delays taken from experimentally determined phase-differences between peak expression of mRNA and protein [4]. Degradation rates were taken from large scale studies of mRNA decay [16, 17] and protein measurements [18–20]. Exponents in transcriptional regulation terms are derived from the number of experimentally verified cis-regulatory elements [4, 21].

The remaining parameters describe the kinetics of transcriptional regulation, which is not known in quantitative detail. Thus, we applied global optimization techniques to fit parameters to carefully normalized qPCR data sets from mouse liver and adrenal gland [8]. For both tissues data from light-dark cycles (LD) and constant darkness (DD) were available. Interestingly, all four expression profiles show clear similarities and thus we fitted a consensus model to these four data sets. The complete set of equations and parameters is provided in (S1 Appendix).

The resulting data-driven gene regulatory network model allows to address the following questions: (i) Which kinetic parameters are most relevant for 24h rhythm generation? (ii) Are oscillations of all five genes necessary for self-sustained rhythms? (iii) What are the most essential regulatory interactions required for rhythm generation? We will answer these questions in the next sections by systematically varying parameters and clamping gene expression levels to their mean values. Thereby, we identify design principles in the network necessary and sufficient for generating circadian oscillations.

Our set of default parameters has been fitted to mRNA expression profiles of circadian clock genes from mouse liver and adrenal gland. It is conceivable that the chosen kinetic parameters are different among tissues and also depend on the specific physiological conditions. In order to test which parameters are most relevant for 24h oscillations, we varied all parameters by two orders of magnitude around the default values. Fig 2 represents the results for four particularly interesting parameters. The periods are plotted for parameter values where self-sustained oscillations occur.

It turns out that oscillations persist for wide ranges of kinetic parameters supporting the robustness of the model. The period increases with the delay between Per2 transcription and its function as an inhibitor (Fig 2A). Indeed, the FASPS mutation of PER2 reduces protein life-time, leading to a faster turn-over and hence to shorter delays [22], thereby implying a shorter period and much earlier phases [23]. Increasing the degradation rate of Cry1 mRNA leads to period shortening as expected (Fig 2B) and consistent with the shorter period of the Cry1^-/- knockout mice [24].

In addition, appropriate interactions of Bmal1, Rev-erb-α, Per2 and Cry1 are required to generate self-sustained rhythms (Fig 2C and 2D). Variations of kinetic parameters associated with transcriptional regulations have minor effects on the period near their default values, consistent with the observation that the clock is resilient to changing transcription rates [25].

The most surprising observation, however, are period jumps for somewhat larger parameter changes (Fig 2). In particular, the detection of long and short periods within a very narrow parameter range suggests that multiple mechanisms might co-exist which can generate self-sustained rhythms. Indeed, the systematic analysis described below allows us to attribute oscillations with different period to specific loops in the model.

For example upon increasing the Per2 delay, the period falls down to 15h after rising up to more than 30h (Fig 2A). Here, the period jump occurs, when the explicit delay is very large (above 10.5h) and influences a subsequent cycle of the oscillations rather than the current one. The rhythm-generating loop, however, remains the same. Further increase of the delay then leads again to an increase of the period until the phase pattern and period before the jump emerges again. An animation of this gradual parameter variation is provided as S1 Video. Interstingly, there is hysteresis near the period jumps, indicating coexisting limit cycles (also termed “birhythmicity” [26]).

A period jump also occurs upon variation of the Cry1 degradation rate (Fig 2B). Here, the jump to a short period of < 10h is associated with a change of the rhythm-generating loop: Cry1 self-inhibition generates these oscillations. Since the self-inhibitory loop exhibits a rather short delay of τ_Cry1 = 3.13h, the resulting period is consequently quite small (comprehensive list of feedback loops and delays in the S4 Appendix). We show in S2 Appendix that in the transition region two rhythms persist (termed “torus”).

If we increase Per2 inhibition by Cry1, oscillations vanish via a supercritical Hopf bifurcation. At much larger parameter values another Hopf bifurcation leads to a limit cycle governed by a double-negative feedback loop involving Per2, Dbp and Cry1.

Taken together, even a relatively small network of just five genes can establish multiple mechanisms generating oscillations, some with periods in the circadian range. While it is generally believed [2] that the negative feedback loop via Per/Cry is the primary driver of circadian oscillations, these multiple regulatory mechanisms even within a relatively small network raise the possibility that the underlying key mechanism for circadian rhythm generation is more complex.

To investigate, which network nodes (genes) are essential for circadian rhythm generation we systematically studied all possible sub-networks under default parametrization.

Our gene network with 7 positive and 10 negative regulations exhibits multiple feedback loops (Fig 1A). Delayed negative feedback loops constitute the basic elements of self-sustained oscillations [27, 28]. Often these negative feedback loops are interlinked with positive feedbacks ensuring robust and tunable rhythms [29–33]. Thus, we focused on which sub-networks forming feedback loops are able to generate sustained rhythms for physiologically plausible parameters.

To this end, we systematically clamped all possible combinations of gene-subsets to their respective oscillation mean values. The mean values are obtained from simulations of the complete network. Clamping the level of Dbp or Bmal1, for example, retains the corresponding positive regulations, but excludes Dbp and Bmal1 as drivers or transmitters of oscillations, thereby focusing on the remaining genes. The clamping to mean values ensures that the system remains near the carefully tuned and physiologically reasonable default configuration. Clamping of genes corresponds to constitutively expressed genes using non-rhythmic promoter constructs [15, 34–36]. Compared to knockout studies, our clamping protocol is less invasive and keeps the system close to its physiological ranges.

There are 5 combinations of 4 genes resulting from clamping only one single gene. For all of the resulting networks there exist certain parameter configurations with oscillatory solutions (blue bars in Fig 3). Clamping Rev-erb-α, Per2 or Cry1 has strong effects: Using the default parameters of the complete network, oscillations vanish. In order to explore the rhythm-generating capabilities of the sub-networks more extensively around the default parameter set, we varied each parameter of the system in a range from 5-fold reduction to 5-fold increase in repeated simulations with 200 points on a log scale. For every simulation we tested, whether or not the sub-network oscillates (Fig 3A). It turns out that in principle all sub-systems of 4 genes are capable of generating oscillations with reasonable periods.

Interestingly, clamping Dbp (first blue bar, Fig 3A) or Bmal1 (last blue bar, Fig 3A) sustains oscillations in about 90% of parameter combinations with a median period close to 24h. This is in line with experimental findings showing that Bmal1 cycling is not necessary for circadian rhythms [36, 37], Dbp^-/- knock-out mice are still rhythmic [38] and triple-knockouts of D-box regulators have only minor effects [39]. Thus, both experimental evidence and our modelling results underline that Bmal1 and Dbp cycling is not essential for sustaining oscillations.

Simultaneously clamping two genes to their mean values results in (53)=10 sub-networks of 3 genes. We find that 5 of these networks are capable of generating self-sustained oscillations, when allowing up to 5-fold adjustments of single parameters. Interestingly, Rev-erb-α is present in most of these oscillatory sub-systems (as an example, see Fig 4B).

Simultaneously clamping 3 genes leads to (52)=10 sets of only 2 remaining genes. Surprisingly, 3 of these pairs are still able to oscillate for appropriate parameter adjustments (Fig 3B). Notably, the negative feedback loop involving Bmal1 and Rev-erb-α oscillates with a period of about 24 hours after only a minor parameter change (compare Fig 4A). It turns out that some of the previously identified oscillations in larger sets of 3 and 4 genes can be traced back to this simple loop. This finding confirms earlier observations that the feedback loop via nuclear receptors can serve as a possible mechanism for rhythm generation [3, 4, 7].

In the previous section gene levels (nodes) were clamped to their mean values, allowing sub-networks to be identified as possible rhythm generators. Now we expand our approach to combinatorial clamping of regulatory interactions (edges in the network graph in Fig 1A) allowing the identification of sub-networks on a process-level. Thereby, network motifs most essential for the generation of 24h rhythms can be identified.

In our model, transcriptional regulations are described by products of activating and inhibiting terms corresponding to the influence of regulating genes [4]. If the expression value of a regulating gene is set constant to its mean value in the term of only one specific target-gene, we call the corresponding interaction “clamped”. For more details on the method, see (S3 Appendix).

Since the gene network contains altogether 17 regulatory interactions, there are 2¹⁷ = 131,072 combinations, or ON/OFF configurations, if OFF means clamping. For all these combinations we tested in detail, whether or not oscillations persist, but did not consider additional variation of kinetic parameters. We found that 14,125 (about 10%) of all network configurations oscillate.

In order to evaluate the importance of specific regulatory interactions we calculated for each interaction the relative frequency of inclusion in the set of oscillatory network configurations. Among all possible configurations any given process is ON or OFF in 50% of the cases. Thus, considering the set of oscillatory ON/OFF configurations, an edge that is not part of the essential loop would still occur in one-half of the cases.

Indeed, the analysis of all oscillatory ON/OFF configurations reveals that most of the processes occur in 50% of the oscillating configurations as expected for a non-essential process. However, a distinct set of regulatory interactions turned out to be present in almost 100% of the oscillating network configurations.

To our surprise, only 3 of the 17 regulatory interactions are exceptionally important to keep the network rhythmic, occurring in almost every oscillating configuration (marked in Fig 5 by thick red lines). All other regulations can be clamped to prevent them from transmitting rhythms: The remaining 3 regulatory interactions still retain oscillations. While the period generated by this 3-gene sub-network in isolation is somewhat longer upon default parametrization, the full network compensates this by fine-tuning through other regulations, including a feedforward loop [40]. Interestingly, the three regulations are all inhibitory: Per2 inhibits Rev-erb-α, Rev-erb-α inhibits Cry1 and Cry1 inhibits Per2.

Such a symmetric triangular motif of inhibitory interactions has been introduced as a paradigm of synthetic oscillators termed “repressilator” [9].

In contrast to most models of the circadian clock, which are essentially based on variations of the Goodwin model [22, 41], the repressilator comprises three subsequent inhibitions rather than a single negative feedback. It is known that classical Goodwin-based models need strong negative cooperativity (minimal Hill coefficient of 8—probably unrealistic biochemically) and long balanced degradation times to obtain self-sustained oscillations [28, 33, 42]. Within the repressilator, however, the delay and the required non-linearities can be distributed over the three inhibitions.

To test the compatibility of the repressilator, we performed a robustness analysis of two simple prototypic models with a single feedback loop, one with only one inhibition and one with the repressilator motif (for details see S5 Appendix). In particular, we generated random parameter sets for both models and compared the frequencies of self-sustained oscillations and the minimal Hill coefficients necessary to generate oscillations (see Fig. 3 in S5 Appendix). We found that the repressilator model has a higher fraction of oscillations and can oscillate with fairly low Hill coefficients of about 2. Note, that modified Goodwin oscillators with additional nonlinearities allow reductions of the Hill coefficient as well [43, 44]. Generally, systems with multiple nonlinearities and delayed feedbacks allow robust oscillations with reasonable Hill coefficients [32, 45, 46]. The repressilator motif allows to distribute nonlinearities and delays.

The repressilator motif is represented as a serial inhibition of Cry1 via Rev-erb-α, of Rev-erb-α via Per2 and of Per2 via Cry1. The two activators Bmal1 and Dbp can be clamped to their mean values without loosing oscillations.

It has indeed been reported that constant Bmal1 levels can sustain rhythms [36, 37] and that triple-knockouts of D-box regulators have only minor effects on circadian rhythmicity [39]. In contrast, double-knockouts of Cry, Per and Rev-erb genes lead to arrhythmicity [5, 24, 47] supporting our finding that circadian rhythms are not generated by a single negative feedback loop, but by a gene regulatory network with a repressilator as a core motif.

Double-knockouts induce behavioral arrhythmicity. Since core clock genes oscillate in surprisingly similar phase relationships in almost all tissues [8, 59, 61], it is very likely that the KO experiments imply also tissue arrythmicity. Indeed, studies of double-knockouts include data on arrhythmic tissues and cells [5, 47, 62].

In previous studies, models have been adapted to available mutant phenotypes [3, 32, 63]. Since our variables group together genes with similar regulatory effects, a comparison with knockout data is not easy. Our clamping protocols resemble constitutive expression or overexpression, and thus we discuss related experiments. It has been shown that constitutive expression or overexpression of Per genes impairs rhythms [34, 35, 64, 65]. Similarly, constitutive or out-of-phase expression of Cry1 impairs rhythmicity [15] and overexpression of Cry1 leads to arrhythmicity [58]. Knockouts and knockdowns of Cry1 lead to arrhythmicity in tissues and cells [62, 66], even though the coupling within the SCN can rescue rhythmicity [62] corresponding to a short-period phenotype of Cry1 knockouts [24]. Interestingly, knockouts and knockdowns of Cry2, an early E-box target not regulated by Rev-erb-α, stay rhythmic with large amplitudes [62, 66, 67]. The essential role of Rev-erb-α inhibition of Cry1 is demonstrated by the removal of the intronic ROR-elements leading to early phases of Cry1 and vanishing amplitudes in single cells [14]. In summary, there is strong experimental evidence that the cycling of the 3 repressilator genes is of central importance for a cellular clock.

Our 5-gene model is based on carefully normalized qPCR data of liver and adrenal gland [8]. More recently, expression profiles of 14 different tissues have been published [59]. In all of these tissues the repressilator genes are oscillating with significant amplitudes and with serially ordered phases consistent with the repressilator mechanism (see S6 Appendix). Similar observations were reported by Yamamoto et al. [61].

In addition to mRNA rhythms protein oscillations are relevant to understand regulatory processes. Unfortunately, liver proteome studies could not quantify core clock protein rhythms due to limited resolution [68, 69]. A recent quantification of clock proteins confirms early protein expression of REV-ERBα, followed by peaks of PER2 and CRY1 [70]. Recent ChIP-Seq experiments allow the estimation of binding phases of regulatory clock proteins [5, 6, 13, 71]. It turns out that REV-ERBα binds early (Circadian Time CT = 6–10), followed by PER2 binding around CT16 and CRY1 binding at around CT24. These subsequent binding peaks are fully consistent with the proposed repressilator mechanism.

Our starting point was a gene-regulatory model based on expression profiles of core clock genes in mouse liver and adrenal gland. As shown in Fig 5, the repressilator is the dominant motif of this gene-regulatory network.

However, Figs 3 and 4 illustrate that also other negative feedback loops are capable of generating oscillations. Furthermore, positive feedback loops are known to support rhythm generation [33]. A comprehensive list of loops within our gene regulatory network is given in (S4 Appendix), showing the interrelations and coherence of loops. Our results suggest, that multiple loops support the generation of circadian oscillations, while the repressilator constitutes an essential core mechanism: While the pure repressilator generates oscillations with increased periods, the addition of other regulations including a feedforward loop [40] from Cry1 to Per2 via Dbp tune the period to values of about 24h.

In summary, our comprehensive analysis of a data-driven core-clock model suggests that the synergy of multiple regulatory motifs allows robust and tunable self-sustained oscillations. We further propose, that a series of subsequent inhibitions known as repressilator constitutes a core motif of the mammalian circadian clock gene-regulatory network.