Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks

Cell-autonomous and self-sustained molecular oscillators drive circadian behavior and physiology in mammals. From rhythms recorded in cultured fibroblasts we identified the dominant cause for amplitude reduction as desynchronization of self-sustained oscillators. Here, we propose a general framework for quantifying luminescence signals from biochemical oscillators, both in populations and individual cells. Our model combines three essential aspects of circadian clocks: the stability of the limit cycle, fluctuations, and intercellular coupling. From population recordings we can simultaneously estimate the stiffness of individual frequencies, the period dispersion, and the interaction strength. Consistent with previous work, coupling is found to be weak and insufficient to synchronize cells. Moreover, we find that frequency fluctuations remain correlated for longer than one clock cycle, which is confirmed from individual cell recordings. Using genetic models for circadian clocks, we show that this reflects the stability properties of the underlying circadian limit-cycle oscillators, and we identify biochemical parameters that influence oscillator stability in mammals. Our study thus points to stabilizing mechanisms that dampen fluctuations to maintain accurate timing in peripheral circadian oscillators.


Introduction
Cell-autonomous molecular pacemakers coordinated by pacemaker neurons in the suprachiasmatic nucleus (SCN) drive circadian rhythms in mammals (Reppert and Weaver, 2002;Schibler and Sassone-Corsi, 2002). Whereas central pacemakers were shown to elicit robust circadian firing rhythms with negligible damping (Welsh et al, 1995;Liu et al, 1997), oscillations of mRNA levels in peripheral organs and cell cultures showed marked damping (Yamazaki et al, 2000;Yoo et al, 2004). This raised the question whether SCN oscillations were qualitatively different from those in peripheral tissues, even though the underlying molecular clocks seem to utilize the same genetic components (Yagita et al, 2001;Schibler and Sassone-Corsi, 2002).
Several studies using luciferase and GFP reporters have now clarified the issue. For instance, tissue explants generate robust oscillatory bioluminescence signals that damp out after B7 days, whereas rhythms in cultured cells derived from the same tissues persist for B20 days . Repeating the experiment in SCN lesioned mice does not abolish rhythms but induces phase shifts, indicating the existence of organ-specific synchronization cues probably masked by the SCN under normal conditions. Although this suggested that peripheral clocks can tick independent of a functional SCN, the issue was addressed directly in a series of studies that combined single cells and population assays. Recordings in immortalized (Nagoshi et al, 2004) and primary (Welsh et al, 2004) mouse fibroblasts showed that single cells generated cell-autonomous rhythms. Furthermore, the effect of a serum shock was to resynchronize randomly phased oscillators rather than jumpstarting arrested oscillators (Nagoshi et al, 2004). Mathematical analysis of longer bioluminescence recordings reinforced these observations and confirmed that dephasing was the dominant cause for amplitude reduction. In embryonic cell lines from the zebrafish that were transfected with a zfperiod4luciferase reporter, Carr and Whitmore (2005) observed the resynchronization of single cellular oscillators by light. Importantly, by monitoring individual cells during 6 days, the authors observed that successive periods within a single cell drift in time (Supplementary Figure S2). Consensus from several models thus strongly supports that damping seen in cultures at the population level reflects desynchronization rather than damping of individual oscillators. Moreover, the default unsynchronized state of a cell population consists of individual oscillating cells with random phase distributions so that populations appear globally non-oscillating. Serum pulses merely resynchronize the phases without starting new oscillations. Finally, all referred work hypothesized that interoscillator coupling was absent in cell cultures, a result that was consistent with relatively short coculture experiments (Nagoshi et al, 2004) and also found in cyanobacterial colonies (Mihalcescu et al, 2004). Nevertheless, the relevance of intercell interaction in peripheral clocks remains open.
Here, we present a general methodology to study biomolecular oscillators and deduce informative parameters from population recordings or individual cells. We investigate the combined effects of limit-cycle stability, intrinsic cellular fluctuations, and oscillator coupling using a compact stochastic mathematical model. Specifically we study the consequences of noisy frequencies and phase coupling on the collective phase dynamics in populations of peripheral circadian oscillators. Using two independent bioluminescence data sets from Nagoshi et al (2004) and Welsh et al (2004), we show that our low-dimensional model captures the data nicely. Our formulation also allows estimating the intercellular coupling strength; we find that whereas the coupling strength is insufficient for synchronization, phase crosstalk between cells can occur at a low rate. Furthermore, we predict a new time scale of about 1 day describing the stiffness of individual circadian frequencies, a quantity that also directly probes the stability of the autonomous oscillator. Finally, we identify biochemical parameters that influence oscillator stability in two models of mammalian circadian clocks.

Results and discussion
Model for interacting noisy-phase oscillators Because they contain a low number of relevant parameters, phase oscillators have been useful to study collective synchronization, phase shifting, and entrainment properties of circadian oscillators (Winfree, 1967;Garcia-Ojalvo et al, 2004;Mihalcescu et al, 2004;Roenneberg et al, 2005). Our previous model for populations of phase oscillators assumed that each cell has a randomly chosen static frequency. In addition, each oscillator could damp out so that its peakto-trough amplitude ratio would decrease exponentially with a time scale T. Using bioluminescence recordings from whole cell cultures, this study showed that the primary cause for amplitude loss was detuning owing to the frequency dispersion s f and not the decay of individual oscillators. This was indicated by the long decay time T (18.8 days) found to be comparable to the experiment duration (Nagoshi et al, 2004).
The frequency drifts reported by Carr and Whitmore (2005), together with our observation that the static frequency model underestimated the frequency dispersion measured in individual cells with a YFP reporter (Nagoshi et al, 2004), prompted us to extend our model to individual oscillator frequencies that drift in time. Additionally, we explicitly model intercellular coupling instead of assuming that it can be neglected. Briefly, the phase derivative for each cell is taken as a stochastic time-dependent frequency plus coupling term ( Figure 1A). Specifically, the frequencies are modeled as an Ornstein-Uhlenbeck (OU) process. The latter is commonly used in the cellular context (Garcia-Ojalvo et al, 2004;Suel et al, 2006) where fluctuations are expected to be correlated in time owing to finite half-lives of other cell components. OU processes are the simplest generalization of Gaussian white noise that introduce exponentially decaying time correlations with a  and frequencies f i (t) describes coupled circadian phase oscillators. The total luminescence signal s(t) is the sum of a population of initially N 0 oscillators each contributing a cosine signal centered around a constant A with relative amplitude B. Cell death follows a Poisson process with time constant t reflected by the indicator variable y i t (t) taking value 1 before (and 0 after) cell i has died. The timedependent frequencies and phases of the individual oscillators are subject to a stochastic differential equation (cf. Materials and methods and Supplementary information). (B) Sample frequency trajectory; g and s f 2 are free constants representing the inverse memory of the frequency trajectories and the frequency dispersion, respectively. (C) Parameter listing. K describes the intercellular coupling between the phases and is taken as all-to-all. More realistic coupling geometries are considered in Figure 3. characteristic time scale g ( Figure 1B; for introduction see Lemons, 2002). Only few parameters are introduced. First s f , the frequency dispersion, is a measure of noise. More precisely, it describes how the limit cycle is susceptible to noise sources, for example, intrinsic noise. On the contrary, g is a property of the deterministic system (when Z¼0 in Figure 1A), which reflects the stiffness of the frequencies or, more generally, the stability of the oscillator. By varying g the model interpolates smoothly between static frequencies (small g), whereas for larger g the frequencies change rapidly and the phase dynamics resembles a diffusion process (Supplementary Figure S2), as in Mihalcescu et al (2004). The coupling among phases is described by the parameter K. With coupling (K40), the model becomes far more complicated, but assuming all-to-all coupling, we derived an expression for the critical coupling value K c above which the population synchronizes (Rougemont and Naef, 2006). The behavior of the synchronization threshold is recapitulated in Supplementary Figure S3, and reflects that rapidly drifting frequencies are easier to synchronize than stiff frequencies.
Other fluctuations could influence the biolumiscence signals. For example, amplitude fluctuations have been considered by Mihalcescu et al (2004). However, these do not affect the estimation of the dephasing parameters g and s f from population-averaged signals. Namely, if A i (t) describes the time-dependent amplitude of cell i, we only need to assume that its fluctuations are independent of the time of cell death (described by the random variable y i t (t) in equation (2); Figure 1) and the phases j i (t). Moreover, a sufficiently large number of cells is required (cf. Supplementary information), which is verified empirically by the good fit in Figure 2A (inset). Moreover, we estimated that 5000 cells contribute to the signal at the end of the recording time (Nagoshi et al, 2004). The same arguments hold for fluctuations in the relative amplitude B. However, relative amplitude fluctuations play a role in the analysis of autocorrelation from single-cell recordings ( Figure 2D).

Frequency dynamics in cell-autonomous oscillators
To establish whether the new model accurately describes bioluminescence signals, we analyzed two independent data sets (from Figure 3B in Nagoshi et al (2004) and Figure 3C, luminometer track in Welsh et al, 2004). Hereafter, we refer to theses data as D1 and D2, respectively. Both used cultured fibroblast cells; however, the first were from the immortalized NIH3T3 line, whereas the second were dissociated from tails of knock-in mPer2 LuciferaseÀSV40 mice. The 19-day bioluminescence recording D1 ( Figure 2A) uses a Bmal1 luceriferase reporter. To connect model and theory, we fit the detrended signal Z(t) to the predicted population-averaged signal ( Figure 2B). Fitting the two data sets leads to nicely compatible values for the parameters s f and g describing the individual oscillators ( Figure 2B, Supplementary Figure S5 and Table I). Note that all parameters can be estimated reliably, and the error bars indicate that the model does not overfit the data. We discuss the results for D1 in some detail. The new model estimates a frequency dispersion of 0.1 per day, which converted to hours leads to a standard deviation in the periods of 2.4 hrs, which is very close to the 2.9 h measured in single cells (Nagoshi et al, 2004) and more accurate than the previous estimate (0.93 h) based on static frequencies. Furthermore, the estimated g¼0.6470.17, reflecting a frequency damping time of 1.56 days, is consistent across data sets. This implies that frequency disturbances take longer than a period length to decay; in other words the initiation of a new cycle does not fully erase the previous cycle.
Whereas s f and g are comparable in both data sets, the cell half-lives and oscillatory amplitude B differ significantly. For example, it appears that the primary cell cultures from D2 are longer lived, with an estimated decay significantly longer than the 3.2 days from D1. A probable explanation is that the primary culture was grown under much richer serum (10%) condition than the immortalized fibroblasts (0.5%), and this drastically affected cell lifespan (U Schibler, personal communication). On the other hand, the oscillatory amplitude is larger in D1 (B¼0.9 versus 0.26). Here, it is likely that synthesis and degradation kinetics of the reporter transcript and protein play a role, as these clearly determine circadian amplitude. For example, the protein half-life of a rhythmically transcribed gene must be short enough for rhythmic protein accumulation to be detected. It is possible that the mRNA half-life of the fusion protein in the mPer2 LuciferaseÀSV40 mouse is longer than in the Bmal-luciferase reporter ( Figure 2B). For instance, the mRNA amplitudes reported by Gachon et al (2004) (Figure 1D and F) in the liver are approximately 5-fold for per2 and 10fold for Bmal1, which is similar to our estimates, but the stability of the luciferase fusion gene should be predominantly determined by the luciferase 3 0 UTR rather than the per2 3 0 UTR.

Oscillators exchange subthreshold phase signals
Our approach also allows to estimate intercell coupling. Previous work addressed the phase coupling among bacterial colonies using a model for two coupled phase oscillators (Mihalcescu et al, 2004); here we use our recent result for coupling in oscillator populations (Rougemont and Naef, 2006). From D1 and D2, we find that the coupling strength is subthreshold and therefore cannot synchronize the cells, consistent with previous coculture (Nagoshi et al, 2004) and transplantation experiments (Guo et al, 2006). How could this small phase crosstalk be mediated in cell cultures? Whereas coupling in SCN neurons depends on synaptic transmission (Liu and Reppert, 2000;Yamaguchi et al, 2003;Ohta et al, 2005;Maywood et al, 2006), coupling in peripheral circadian clocks could occur via paracrine signaling, for example, a broad class of signaling cues were shown to elicit circadian rhythms in fibroblasts (Balsalobre et al, 2000a,b). It is possible that some residual form of such cues could be active in culture. The measured coupling (K¼0.05 day À1 ) signifies that neighboring cells are able to shift the phase of an oscillator by half a cycle in 10 days. In comparison, synchrony would develop if the same delay could be mediated in less than about 5 days, assuming fixed s f and g ( Figure 2C, left). Alternatively, synchrony would occur if the stiffness was reduced (g increased) or the frequency dispersion was reduced by about twofold (Figure 2C, right).  Table I (g¼0.9 day À1 and s f ¼0.1 day À1 ). Amplitude fluctuations were modeled as a correlated process with mean B, g B ¼5g, and s B /B¼0.4, leading to a rapidly decreasing initial transient in the envelope (exact prediction in blue; cf. Supplementary information). The approximation for large g B used to fit panel A is shown is cyan. The short (dephasing) and long-time (phase diffusion) regimes are indicated in red and green, respectively. In the calculations presented so far, we have assumed that each cell was evenly coupled to every other cell (all-to-all coupling). This choice reflected mathematical convenience, as we can then derive expressions for the critical coupling strength. To investigate the validity of this approximation, we assume that synchrony is mediated through diffusing peptides. With sizes of about 1 nm, we can estimate using the Einstein-Stokes formula that 10 cell shells can be reached within 12 min, which is small compared with the oscillatory period. Therefore, if signaling peptides are secreted, each cell effectively receives temporally coherent synchronization signals from a large number of neighboring cells. To assess whether this is sufficient for synchronization, we have assumed realistic cell geometries generated from random Voronoi tilings ( Figure 3A). In these, each cell has an average of six neighbors, but the environment for each cell is slightly different. Simulations show that adding cell shells leads to a synchronization behavior that converges to the all-to-all model ( Figure 3B), indicating that our estimates of coupling are reliable. In fact the coupling estimated from the all-to-all model is a lower bound for the real coupling strength.

Frequency dynamics from individual cell recordings
Next we analyze the frequency dynamics in single-cell luciferase recordings (cf. Figures 2B and 3C in Welsh et al, 2004) using autocorrelation of the signals. A property of autocorrelations is that these are unaffected by coupling as long as K is below the synchronization threshold (section 1.5, Supplementary information). This highlights an interesting difference between single cell and population recordings: subthreshold coupling cannot be estimated from single cells, whereas it can be estimated from the dephasing dynamics of initially synchronized populations, as in Figure 2C. One salient feature of individual cell bioluminescence signals is that their amplitude fluctuates significantly. The short time behavior of autocorrelations shows abrupt drop in the autocorrelation within the first period (arrow in Figure 2D), which is captured by a simple model, assuming independent amplitude and phase fluctuations ( Figure 2E and Supplementary information, section 1.3). Even though we had few cells at our disposal, fitting the mode to the data leads to an estimated frequency dispersion s f within a factor two of the population estimate, whereas the drift parameter g is within the error bars of the population estimate. This good agreement (cf. overlay in Supplementary Figure S6) from different approaches thus supports that correlated frequency fluctuations are an essential signature of peripheral circadian oscillators.

Origins of frequency fluctuations
Phase models, which have been popular in circadian biology since the work of Winfree (1967), postulate that molecular mechanisms generate sustained oscillations without describing the molecular interactions among clock components. However, the wealth of biochemical data about circadian pathways has allowed the development of rate equation models that show how limit-cycle oscillation can arise (Forger and Peskin, 2003;Leloup and Goldbeter, 2003;Locke et al, 2005) and resist noise in clock circuits (Barkai and Leibler, 2000;Gonze et al, 2002;Vilar et al, 2002;Forger and Peskin, 2005;Gonze and Goldbeter, 2006). Several of these models rely on the mutual feedback of an activator and repressor pair that triggers relaxation oscillations (Barkai and Leibler, 2000;Vilar et al, 2002), whereas others use delayed feedback (Gonze et al, 2002). Under physiological conditions, such chemical reaction networks face two types of fluctuations: (i) those that follow from the finite numbers of molecules, DNA, mRNA, or proteins (often termed intrinsic) and (ii) fluctuations that act globally on all parts of the network like thermal fluctuations or changes in cell volume in dividing cells (termed extrinsic). Recent stochastic simulations probing the role of intrinsic noise in two mammalian clock models showed that the period variability is directly affected by the number of clock proteins. Specifically, the variance in circadian period length is inversely proportional to the number of molecules (Leloup and Goldbeter, 2003;Forger and Peskin, 2005). Whereas intrinsic fluctuations thus lead to period variability, the magnitude of these fluctuations depends on the stability of the deterministic limit cycle. This is reflected in our phase model by the relation s f 2 ¼s Z 2 /(2g), that is, for fixed noise s Z 2 , more stable limit cycles (g large) are less susceptible to noise and thus have smaller frequency dispersion. As g measures the stability of the phase oscillator, it is natural to ask how this inverse timescale relates to a canonical measure of limit-cycle stability, namely the negative logarithm of the leading Floquet multiplier m F (Eckmann and Ruelle, 1985;Strogatz, 2000). In Supplementary Figures S7-S9, we use g F ¼Àlog(m F ) to measure limit-cycle stability. To address the above question we apply Floquet analysis to a generic circadian model describing a delayed negative feedback loop (Gonze and Goldbeter, 2006). This model includes three dynamical variables and describes mRNA transcription, translation, and nuclear translocation of an autoregulatory clock protein. In parallel we simulate intrinsic noise through a master equation and use the Gillespie (1976) algorithm to simulate trajectories (as in Gonze and Goldbeter, 2006) from which we estimate the variability in frequency. By varying the number of molecules (O) and transcription rate, we found that population averages of stochastic trajectories could be well approximated by the form predicted for the phase model (Supplementary Figure S9), with parameter values g¼g F and s f very close to the estimate from the trajectories. This correspondence shows that the Floquet stability can be probed experimentally using our method, and that we can use powerful analytical tools to study how limit-cycle stability depends on model parameters. For illustration, we vary each of the nine parameters independently while monitoring stability and period of the limit cycle (Supplementary Figure S7). We find that stability can be increased most efficiently by raising the transcription rate (v s ) of the mRNA or by reducing the half-max parameter (K m ) for mRNA degradation. Meanwhile, increasing the translation rate of the protein (k s ) or the nuclear translocation rate (k 1 ) reduces period length as expected, whereas period lengthens with increased v s .
Finally, we apply stability analysis to a detailed sixteendimensional model for the mammalian circadian clock (Leloup and Goldbeter, 2003). In the mammalian clock, the principal activators are the Bmal1, Clock, and Npas2 transcription factors, whereas repression is mediated by Per1, Per2, Cry1, Cry2, and RevErba (Schibler and Naef, 2005). The model by Leloup and Goldbeter (2003) (LG) is based on a merged Per gene, a merged Cry, and Bmal1. It also describes protein phosphorylation, for example, Per is phosphorylated by casein kinase 1 (CK1), complex formation, and transport from the cytoplasm to nucleus where mRNAs are transcribed. The model has 53 parameters, but we restricted ourselves to four groups of parameters: transcription rates, phosphorylation rates of the Per and Cry proteins, complex formation between Per and Cry, and nuclear translocation rates (Supplementary Figure S8). Sensitivity analysis for the period length in this model was detailed by Leloup and Goldbeter (2004). Among the changes that most affect period length, increased complex formation rate and nuclear entry rate shorten the period as expected, whereas increased Bmal1 transcription or Period phosphorylation lengthens the period (Supplementary Figure  S8). Floquet analysis shows that only few of the parameters tested strongly influence limit-cycle stability, and most of the changes lead to less stable oscillators. Note that the value of g F (0.87 inverse periods) for the nominal parameters are in the range of our measured g ( Table I), suggesting that the LG model has realistic limit-cycle stability properties. In contrast, the above three-variable model has a more stable limit cycle with g F ¼2.65 inverse periods. Interestingly, the parameter that most increases the stability is the phosphorylation rate of Per (Supplementary Figure S8, star). Thus, this predicts that within the LG model, the observed frequency dispersion could be reduced maximally by overexpressing CK1 such that the phosphorylation rate of Per would be increased by a factor of 3.

Conclusion and outlook
We showed how circadian bioluminescence signals recorded in peripheral clock cells can be analyzed to provide insight into three essential aspects of circadian clocks: limit-cycle stability, their susceptibility to fluctuations, and intercellular coupling. Analyzing independent bioluminescence recordings leads to consistent parameter for frequency dispersion, frequency stiffness, and coupling between cells. Furthermore, estimates from populations and single cells were in good agreement. Interestingly, our study predicted that oscillator stability is such that frequencies in individual cells remain correlated beyond one circadian cycle length. Additionally, we estimated phase crosstalk in cell cultures, which indicated that the coupling strength was nonzero but only about half of that needed to synchronize the cells. We showed that collective synchronization would occur if the frequency dispersion would be tighter, for example, if it could be reduced by a factor of about 2 ( Figure 2C). One possible role of residual but subthreshold coupling in peripheral circadian oscillators is that it facilitates entrainment by systemic cues, which are restricted to time windows shorter than the period (supported by simulations, data not shown). Thus, our phase model shows how dynamical stability, cellular noise level, and intercellular coupling shape collective and individual bioluminescence rhythms in peripheral circadian oscillators. We ended by linking the stability of limit cycles to the properties of biochemical oscillator models and pointed toward molecular mechanisms that are predicted to increase the stability of the circadian clock in mammals.
Biochemical oscillators commonly serve to coordinate cellular processes occurring over a wide range of timescales. We showed that peripheral circadian oscillators are weakly coupled, but others interact strongly to elicit collective synchronization. The latter include the respiratory cycle in yeast grown under constant condition (Henson, 2004;Klevecz et al, 2004;Locke et al, 2005) or the somite clock, as measured in tissue explants (Horikawa et al, 2006;Masamizu et al, 2006). The approach presented here is naturally suited to study fluctuations in these collectively synchronized biological timekeepers. As luminescence reporters are becoming widely used in chronobiology, we expect that it will provide a compact framework to compare the dephasing dynamics in a broad class of molecularly distinct oscillators operating in fluctuating cellular environments.

Data sets
We analyze two independent circadian bioluminescence data sets (from Figure 3B in Nagoshi et al (2004) and Figure 3C, luminometer track in Welsh et al (2004)). We refer to these data as D1 and D2, respectively.

Phase model
Drifting frequencies with the properties of a fixed mean E[f]¼m and variance var[f]¼s f 2 are generated through an Ornstein-Uhlenbeck process described by the stochastic differential equation in which Z(t) is a Gaussian white noise source with variance parameterized as s Z 2 ¼2gs f 2 . Here, g and s f 2 are free constants representing the decay rate and amplitude of frequency fluctuations, respectively. The mathematical aspects of synchronization in this model are found in Rougemont and Naef (2006). Additional results used in this paper, for example, the form for the autocorrelations, are derived from Supplementary information.
Simulations were performed with the discrete dynamical updates (Lemons, 2002): f i ðt þ dtÞ ¼f i ðtÞe Àgdt þ m f ð1 À e Àgdt Þ þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À e À2gdt p z; f i ðt þ dtÞ ¼f i ðtÞ þ dtðf i ðtÞ þ K X j2Ni sinðf j À f i ÞÞ where N i denotes the neighbors of cell i and z is a random number drawn from a Gaussian distribution with mean 0 and variance s f 2 . Initial conditions f i (t¼0)¼0 reflect synchronization by the serum pulse and f i (t¼0) was drawn for a Gaussian with mean m f and variance s f 2 . For Figure 2B and C, the number of oscillators (N¼5000) was sufficient to prevent finite size effects from influencing parameter estimation in the time span from 0 to 10 days.

Regression and detrending
To study the effects of phase dynamics, we work with the detrended variable ZðtÞ ¼ sðtÞ À ANðtÞ ANðtÞ which renormalizes the signal for cell death and amplitude A (Nagoshi et al (2004); Supplementary information). N(t) denotes the population size at time t and Z(t) lies between À1 and þ 1. All regression and associated statistical tests were performed with the R language for statistical computing and graphics (http://cran. r-project.org) using the lm and nlm routines. Simulations were implemented in R and C. For the parameter estimations in Table I, a routine, that simulated the averaged signal from N¼5000 oscillators was passed to the nlm routine to fit the detrended signal Z(t). Posterior likelihoods of the parameters in Figure 2 were generated using a grid of parameters in the (K, g), (K, s f ), (g, s f ) planes that were passed to the same simulation routine.

Voronoi tilings
Voronoi tilings for randomly seeded cell nuclei were used to define neighbors in 2D cell arrangements. First and second neighbors were considered for the coupling term in Figure 3. Details about the geometrical algorithm are given in Section 4 of Supplementary information.

Floquet analysis
We compute the Floquet multipliers for the Gonze and Leloup models using the continuation software AUTO (Doedel et al, 2001). To find the limit cycles, initial transcription rates were set to zero, for which the trivial fixed point with all concentrations equal to 0 is stable. Parameters were then increased one by one to their nominal values (those used in the original articles) whereas limit cycles solutions were tracked along with their periods and Floquet multipliers.

Supplementary information
Supplementary information is available at the Molecular Systems Biology website (www.nature.com/msb).