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 ; Liu ), oscillations of mRNA levels in peripheral organs and cell cultures showed marked damping (Yamazaki ; Yoo ). 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 ; 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 ∼7 days, whereas rhythms in cultured cells derived from the same tissues persist for ∼20 days (Yoo ). 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 ) and primary (Welsh ) 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 ). 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 zfperiod4-luciferase 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 ) and also found in cyanobacterial colonies (Mihalcescu ). 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 and Welsh , 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 ; Mihalcescu ; Roenneberg ). 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 peak-to-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 σ 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 ).

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 ), 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 ; Suel ) 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 characteristic time scale γ (Figure 1B; for introduction see Lemons, 2002). Only few parameters are introduced. First σ, 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, γ is a property of the deterministic system (when η=0 in Figure 1A), which reflects the stiffness of the frequencies or, more generally, the stability of the oscillator. By varying γ the model interpolates smoothly between static frequencies (small γ), whereas for larger γ the frequencies change rapidly and the phase dynamics resembles a diffusion process (Supplementary Figure S2), as in Mihalcescu . The coupling among phases is described by the parameter K. With coupling (K>0), the model becomes far more complicated, but assuming all-to-all coupling, we derived an expression for the critical coupling value Kc 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.

Figure 1

Stochastic phase model with drifting frequencies and intercellular phase coupling. (A) An extended Kuramoto model for the oscillator phases φ(t) and frequencies f(t) describes coupled circadian phase oscillators. The total luminescence signal s(t) is the sum of a population of initially N0 oscillators each contributing a cosine signal centered around a constant A with relative amplitude B. Cell death follows a Poisson process with time constant τ reflected by the indicator variable θτ(t) taking value 1 before (and 0 after) cell i has died. The time-dependent 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; γ and σ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.

Other fluctuations could influence the biolumiscence signals. For example, amplitude fluctuations have been considered by Mihalcescu . However, these do not affect the estimation of the dephasing parameters γ and σ from population-averaged signals. Namely, if Ai(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 θτ(t) in equation (2); Figure 1) and the phases ϕ(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 ). 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).

Figure 2

Analysis of bioluminescence data by Nagoshi et al (referred to as D1). (A) Raw data reproduced from Nagoshi . Inset: logarithmic scale emphasizes the exponential signal decrease reflecting cell death with half-life 3.22 days. (B) Maximum likelihood fit of the detrended signal Z to our model. The data were detrended using band-pass filtering as detailed by Nagoshi . (C) Posterior likelihoods of the parameters. Projections for each pair of model parameters γ, σ, and K are shown: red indicates high probability; standard errors around the maximum likelihood parameters are indicated (cf. Table I). The critical coupling lines (black) with fixed third parameter indicate that the coupling should be increased for synchrony (first two panels), or alternatively the frequency dispersion should be reduced (third panel). (D, E) Frequency drifts from bioluminescence signal in individual cells from the autocorrelation analysis of 10 individual cells (from Welsh , Figures 2B and 3C). (D) Each color is the autocorrelation for one cell. For convenience the absolute value is plotted and the maxima are fit to the envelope

predicted from the model (cf. panel E and Supplementary information) which leads to a best fit (black line) with C=0.72, σ=0.059±0.003 day−1, and γ=0.39±0.11 day−1. (E) A single cell was simulated for 30 days with parameters similar to those for the data D2 in Table I (γ=0.9 day−1 and σ=0.1 day−1). Amplitude fluctuations were modeled as a correlated process with mean B, γ=5γ, and σ/B=0.4, leading to a rapidly decreasing initial transient in the envelope (exact prediction in blue; cf. Supplementary information). The approximation for large γ 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.

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 and Figure 3C, luminometer track in Welsh ). 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 mPer2Luciferase−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 σ and γ 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 ) and more accurate than the previous estimate (0.93 h) based on static frequencies. Furthermore, the estimated γ=0.64±0.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.

Figure 3

Comparison of synchronization behavior between the all-to-all model and random 2D cell arrangements. (A) 2D cell culture geometries are modeled as random Voronoi tessellations of the square. First and second neighbors of the central cell (in yellow) are indicated in green and red, respectively. (B) Synchronization transition for cell cultures with local coupling. The synchronization parameter R∞ is computed as a function of the coupling strength K for three different geometries: all-to-all coupling (red), nearest-neighbor coupling (blue), and coupling extending to second nearest neighbor (purple). For each coupling strength and geometry, 10 000 cells were simulated. Mean and standard error are represented for five independent cell arrangements (details in Supplementary information, section 4).

Table 1

Parameter estimates in two independent data sets

	Half-life (days)	K (day−1)	B	Period 1/μf (h)	γ (day−1)	σf (day−1)
Nagoshi et al (D1)	3.22±0.01**	0.05±0.02*	0.90±0.01**	25.75±0.02**	0.64±0.17**	0.1±0.007**
Welsh et al (D2)	55.2±1.3**	0.02±0.08	0.26±0.003**	25.48±0.03**	0.89±0.64	0.1±0.01**

Nonlinear expression for the predicted Z(t) (cf. Supplementary information) was fit to the luminescence signals (Figure 2C and Supplementary Figure S5) to estimate K, B, γ, μ, and σ. Exponential trend (cell half-life) was estimated from the linear regression in Figure 2A (inset). * and ** refer to fit parameters with P<0.01 and P<0.001, respectively. Standard errors in the parameter estimates are also indicated.

Whereas σ and γ 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 mPer2Luciferase−SV40 mouse is longer than in the Bmal-luciferase reporter (Figure 2B). For instance, the mRNA amplitudes reported by Gachon (Figure 1D and F) in the liver are approximately 5–fold for per2 and 10-fold for Bmal1, which is similar to our estimates, but the stability of the luciferase fusion gene should be predominantly determined by the luciferase 3′UTR rather than the per2 3′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 ); 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 ) and transplantation experiments (Guo ). 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 ; Ohta ; Maywood ), 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 , 2000b). 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 σ and γ (Figure 2C, left). Alternatively, synchrony would occur if the stiffness was reduced (γ increased) or the frequency dispersion was reduced by about twofold (Figure 2C, right).

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 ) 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 σ within a factor two of the population estimate, whereas the drift parameter γ 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 ) and resist noise in clock circuits (Barkai and Leibler, 2000; Gonze ; Vilar ; 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 ), whereas others use delayed feedback (Gonze ). 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 σ2=ση2/(2γ), that is, for fixed noise ση2, more stable limit cycles (γ large) are less susceptible to noise and thus have smaller frequency dispersion. As γ 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 μF (Eckmann and Ruelle, 1985; Strogatz, 2000). In Supplementary Figures S7–S9, we use γF=−log(μ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 (Ω) 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 γ=γF and σ 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 (vs) of the mRNA or by reducing the half-max parameter (Km) for mRNA degradation. Meanwhile, increasing the translation rate of the protein (ks) or the nuclear translocation rate (k1) reduces period length as expected, whereas period lengthens with increased vs.

Finally, we apply stability analysis to a detailed sixteen-dimensional 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 RevErbα (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 γF (0.87 inverse periods) for the nominal parameters are in the range of our measured γ (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 γ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 ; Locke ) or the somite clock, as measured in tissue explants (Horikawa ; Masamizu ). 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.

Materials and methods

Data sets

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

Phase model

Drifting frequencies with the properties of a fixed mean E[f]=μ and variance var[f]=σ2 are generated through an Ornstein–Uhlenbeck process described by the stochastic differential equation

in which η(t) is a Gaussian white noise source with variance parameterized as ση2=2γσ2. Here, γ and σ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):

where N denotes the neighbors of cell i and ζ is a random number drawn from a Gaussian distribution with mean 0 and variance σ2. Initial conditions φ(t=0)=0 reflect synchronization by the serum pulse and f(t=0) was drawn for a Gaussian with mean μ and variance σ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

which renormalizes the signal for cell death and amplitude A (Nagoshi ; 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, γ), (K, σ), (γ, σ) 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 ). 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.