A stochastic model for circadian rhythms from coupled ultradian oscillators

Abstract

Background

Circadian rhythms with varying components exist in organisms ranging from humans to cyanobacteria. A simple evolutionarily plausible mechanism for the origin of such a variety of circadian oscillators, proposed in earlier work, involves the non-disruptive coupling of pre-existing ultradian transcriptional-translational oscillators (TTOs), producing "beats," in individual cells. However, like other TTO models of circadian rhythms, it is important to establish that the inherent stochasticity of the protein binding and unbinding does not invalidate the finding of clear oscillations with circadian period.

Results

The TTOs of our model are described in two versions: 1) a version in which the activation or inhibition of genes is regulated stochastically, where the 'unoccupied" (or "free") time of the site under consideration depends on the concentration of a protein complex produced by another site, and 2) a deterministic, "time-averaged" version in which the switching between the "free" and "occupied" states of the sites occurs so rapidly that the stochastic effects average out. The second case is proved to emerge from the first in a mathematically rigorous way. Numerical results for both scenarios are presented and compared.

Conclusion

Our model proves to be robust to the stochasticity of protein binding/unbinding at experimentally determined rates and even at rates several orders of magnitude slower. We have not only confirmed this by numerical simulation, but have shown in a mathematically rigorous way that the time-averaged deterministic system is indeed the fast-binding-rate limit of the full stochastic model.

Background

We are concerned with mechanisms that can account for circadian rhythms at the cellular level. Although circadian oscillators exist in complex multicellular organisms as well as in single-cell organisms, it is thought that most occur in single cells [1–3]. We have previously [4] described a model for circadian oscillations in which ultradian oscillators, which have been widely observed to occur in living systems, are coupled to produce circadian periods. The model was based, as is much of the related literature, on so-called transcriptional-translational oscillators (TTOs), in which genes are activated or inhibited for transcription by protein products of the oscillating system itself (transcriptional activators or suppressors, respectively). Several models for interactions between more than one oscillator to generate a circadian one have been described [5–7], but ours differs in positing coupling between the protein products of independent ultradian oscillators. We argued that our model provides a plausible evolutionary origin for circadian oscillators across a range of organisms, since it allows existing ultradian oscillations to be co-opted as components of circadian oscillators without disturbing their primary functions.

A challenging feature of TTOs is the fact that in cells, a given transcribed gene is present in one, or at most a small number, of copies, and its interaction with a transcriptional regulator is not correctly modeled by deterministic differential equations as used in [4]. Rather, because the number of copies of an expressed gene, and at some times the numbers of transcription factor molecules in a cell is small, such interactions are more accurately described by stochastic equations, and this has been done for a number of existing models [6, 8–10], using a classical algorithm due to Gillespie [11]. In some cases this results in shorter autocorrelation times [8] or random fluctuations [6]. Typically, as for example in reference [9], the effect of the stochasticity is to degrate the circadian oscillations, but for fast enough binding rates, the circadian oscillations are maintained. Our objective here is to apply the stochastic approach to a model similar to the one described in [4]. To do this, we have estimated the rates of association and dissociation of transcription factors from their DNA binding sites. We have then incorporated these rates, together with parameters previously used [4], into the new version of the model, in which the DNA binding steps have been treated as stochastic processes. The subsequent steps of translation and turnover of protein and mRNA have been left as deterministic ones, since the numbers of molecules in these processes are large.

We suggest that if the model is well-behaved with the critical DNA-binding step as a stochastic process, then the remaining steps can be left as deterministic without compromising the reliability of the model. Three quite different time scales arise in the model. The binding and dissociation of the transcription factors to DNA sites occur on a fast time scale, as discussed below. We introduce an (artificial) parameter ε with dimensions of time to adjust the time scales for these events and to explore the limit ε → 0. In our Numerical Tests section we vary ε; for several numerical simulations we use a value that corresponds to a relatively high rate of binding and dissociation, as explained in the Model section below. Under these conditions the results are essentially indistinguishable from the simulations for a time-averaged deterministic model which is obtained in the limit ε → 0. We subsequently show that the model is well-behaved even for binding rates that are at least 1000-fold slower.

The second significant time scale is given by the periods of the individual ultradian oscillators, which are of the order of a few hours. The critical parameters for these oscillations are those describing the half-lives of mRNA, proteins, and protein complexes. Following our numerical tests, we conduct a brief exploratory analysis of the range of periods of our "primary" oscillators.

The third time scale is, of course, the circadian rhythm time scale, which in our model arises from an interaction of two of the simpler ultradian oscillators of slightly different frequencies. Natural selection could explain why pairs of frequencies leading to the right "beats" have emerged in the course of evolution. In fact, the common occurrence of ultradian oscillators would make it easy for evolution to produce circadian rhythms out of different components in different organisms, as is actually observed [4]. This mechanism has the added advantages of robustness and easy adaptability (the period of the beat will change with minor adjustments of the frequency ratio between the two primary oscillators, but this ratio could stay quite stable even if the parameters involved varied with external conditions such as temperature). A power spectrum analysis presented below demonstrates the robustness of the model with respect to the parameter ε. We mention that power spectra could be used to analyze observational data for a potential validation of the model. First steps in this direction were taken in [12].

The model

Our model involves TTOs contained in a single cell. As described in [4], the model comprises two ultradian "primary" oscillators whose protein products are coupled to drive a circadian rhythm. For simplicity, the two coupled primary oscillators are essentially identical, with only their frequencies different, since the critical feature is the ability to couple TTOs through known molecular processes (formation of transcriptional-regulatory protein heterodimers). Therefore, the key question regarding the ability of a stochastic process to describe stable circadian oscillators can be addressed in terms of one primary oscillator. In this system, two genes (DNA sites) are transcribed into mRNA, and this process is the origin of the following chemical dynamics.

Transcription by gene 1 occurs when site 1 (its regulatory region) is unoccupied. Its state is given by a random variable X₁, so that

X₁ = 0 if site 1 is empty; X₁ = 1 if site 1 is occupied by D₂ (see below)

When gene 1 is active it produces mRNA (measured in molecules per cell, R₁) at a constant rate k₁₃. These molecules undergo first-order decay with a rate constant k₁₄.

The mRNA molecules are translated into protein P₁, which: (a) decays at rate constant k₁₆, (b) forms homodimers D₁ at rate k₁₇, and (c) forms heterodimers D₁₃ with proteins P₃ from a third gene (see below) with a rate constant k₆₁.

The homodimer D₁ binds to site 2, and thereby activates the transciption of gene 2. The state of gene 2 is given by the value of a random variable Y₁ so that

Y₁ = 0 if site 2 is empty, and Y₁ = 1 if site 2 is occupied by D₁.

Transcription of gene 2 and translation of its mRNA into protein P₂, which forms homodimer D₂, which in turn feeds back to inhibit gene 1 (above). In addition, the P₂ molecules decay with a certain (biological) half-life.

These linked reactions generate a TTO for an appropriate choice of parameters. The parameters used in our subsequent calculations are listed in Table 1. Our model entails gene 1 being inhibited by homodimer D₂ and gene 2 being activated by homodimer D₁. This is the mechanism leading to primary oscillations.

Table 1 Parameters. The dimensions are [k₁₃] = hr^-1, [k₁₄] = (nr. × hr)^-1, [k₁₇] = (nr.² × hr)^-1, [r] = 1 etc. We assign [ε] = hr, so that r, s,... become dimensionless

We denote by R _i , P _i , D _i , i = 1, 2 the concentrations of the mRNA, the translated protein and the homodimer produced by site i. The above scenario is then summarized in the following system of stochastic differential equations (only two of the equations contain the random variables X₁ and Y₁ explicitly, but all dependent variable are then random variables of necessity). The parameters k₁₃ etc. have the same meaning as in Ref. [4], and we have kept the notation used there; this explains the unconventional numbering (some of the equations from the reference, and hence some of the parameters, are no longer needed).

${R^{\prime }}_1=k_{13}(1-X_1)-k_{14}{R}_1\left(1\right)$

${P^{\prime }}_1=k_{15}{R}_1-k_{16}{P}_1-2k_{17}{P}_1^2+2k_{18}{D}_1-k_{61}{P}_1{P}_3+k_{62}{D}_{13}\left(2\right)$

${D^{\prime }}_1=k_{17}{P}_1^2-k_{18}{D}_1\left(3\right)$

${R^{\prime }}_2=k_{13}{Y}_1-k_{14}{R}_2\left(4\right)$

${P^{\prime }}_2=k_{25}{R}_2-k_{16}{P}_2-2k_{27}{P}_2^2+2k_{28}{D}_2\left(5\right)$

${D^{\prime }}_2=k_{27}{P}_2^2-k_{28}{D}_2\left(6\right)$

The last two terms in the second equation reflect the combination of proteins P₁ and P₃ (which is produced by the second primary oscillator) to form the heterodimer D₁₃. This heterodimer in turn breaks down into pairs P₁ and P₃ at rate constant k₆₂.

The second primary oscillator is given by a nearly identical set of equations, except that the periods of the oscillations are slightly different. This can, of course, be achieved by changing the parameters in many ways, but the simplest method is to have the two TTOs identical in nature but with different time scales. To do this we simply multiply each right hand side by a fixed constant δ > 0, where δ is close (but not identical) to one. For example, the first equation of the second oscillator will read

$${R^{\prime }}_3$$

(1)

= δ(k₁₃(1 - X₂) - k₁₄R₃).

The parameters chosen reflect, where available, reasonable choices of known molecular processes. The critical ones for establishing the periods of the primary oscillators are the decay times of the mRNAs and proteins. For the former, a half-life of 13–17 minutes and for the latter, 4–17 minutes generate ultradian oscillations in the model. The values used in the simulation are given in Table 1.

The coupling between the two sites communicating in each oscillator is, of course, provided by the random variables X _i , Y _i . The times for which these random variables stay constant are assumed to be exponentially distributed. For example,

Prob{X₁ = 0 in (t, t + h)|X₁(t) = 0} = exp (-D₂(t)h/ε) + o(h),

Prob{X₁ = 1 in (t, t + h)|X₁(t) = 1} = exp (-rh/ε) + o(h)

while

Prob{Y₁ = 0 in (t, t + h)|Y₁(t) = 0} = exp (-D₁(t)h/ε) + o(h),

Prob{Y₁ = 1 in (t, t + h)|Y₁(t) = 1} = exp (-st/ε) + o(h).

Here ε is a time scaling parameter, introduced for convenience to exploit the fact that the binding and unbinding of the homodimers occurs on a faster time scale than the remaining processes. The constants r and s measure, relative to the scale ε, the average times for which the sites will remain occupied. As this is an internal parameter of the site it should not depend on the states of the rest of the system (like, for example, the dimer concentrations).

We use ε to gauge the rate constant for binding of the transcriptional-regulatory proteins (D1, D2) to the binding sites on the relevant genes. Experimental work has shown that the second-order rate constant for the binding of transcription-regulating proteins to DNA can be 100 to 1000 times greater than the maximum rate predicted for three-dimensional diffusion [13, 14]. With transcription-regulating protein concentrations measured in molecules/nucleus, using the experimental rate constant for binding of the lac repressor to its cognate DNA [10^-10(Msec)^-1], and assuming that a small eukaryotic nucleus has an effective volume of 40% of its total volume, this suggests a value for ε of 0.10 seconds (2.8 × 10^-5 hours). This can be interpreted as the time required for a binding event when Dl or D2 is present at 1 molecule/nucleus. At higher concentrations (of D1 or D2), this time will shorten proportionately. The average "free" time of the binding site for D₂ is thus ε/D₂, and the average "occupied" time is ε/r. Their quotient is independent of ε, but will change with the homodimer concentration D₂. Similar interpretations apply for X₂ and Y₂ and the random variables associated with the second primary oscillator. We have used the value ε = 0.1 sec for producing most of the numerical simulations in our Numerical Tests Section below (Figures 1, 2, 3, 4). However, as shown in Figures 5 and 6, an ε of 1000 times greater value (corresponding to a 1000-fold slower rate of binding) yields effectively the same power spectrum for the circadian model. This is comparable to the observation by Forger and Peskin that in their model for mammalian circadian rhythms the on/off times need to be in the order of seconds.

Figure 1

The time evolution of the proteins P₁ and P₃ according to the time-averaged model.

Figure 2

The time evolution of the heterodimer D₁₃ and the homodimer D₅ according to the time-averaged model.

Figure 3

The time evolution of the proteins P₁ and P₃ according to the stochastic model.

Figure 4

The time evolution of the heterodimer D₁₃ and the homodimer D₅ according to the stochastic model.

Figure 5

Power spectra for D₅ when ε = 2.8 × 10^-5 (no smoothing).

Figure 6

Power spectra for D₅ when ε = 2.8 × 10^-2 (smoothed with a Daniell filter of length 11).

The average times for which a dimer stays bound (ε/r, ε/s, etc.) are independent of the state of the system. In contrast, the "free" times are inversely proportional to the concentration of the attaching homodimer. In one of our simulations we use r = 25 and ε = 10^-1sec (which corresponds to $\frac{\epsilon }{r}=\frac{1}{250}$ sec, or an average of 900,000 binding events per hour). We shall see that the corresponding stochastic simulation compares well with a limiting scenario for which ε = 0. Before we describe this limiting scenario in detail we present the remaining equations making up the complete oscillatory system.

As stated earlier, the protein products P₁ and P₃ of the first and second primary oscillators combine to produce the heterodimer D₁₃. As formulated in the model, this heterodimer binds to the regulatory site of a fifth gene and activates it for transcription (other constructs, involving other heterodimeric products of the two primary oscillators, and either stimulation or inhibition of transcription of the fifth gene, could also be used). Transcription, translation, and dimerization of the protein product of gene 5 yields the product D₅, which is the primary circadian output of the model (although all variables show circadian behaviour to a greater or less extent, as seen in the graphical results).

The corresponding system is

${D^{\prime }}_{13}=k_{61}{P}_1{P}_3-k_{62}{D}_{13}\left(7\right)$

${R^{\prime }}_5=k_{53}{X}_3-k_{54}{R}_5\left(8\right)$

${P^{\prime }}_5=k_{15}{R}_5-k_{16}{P}_5-2k_{57}{P}_5^2+2k_{58}{D}_5\left(9\right)$

${D^{\prime }}_5=k_{57}{P}_5^2-k_{58}{D}_5,\left(10\right)$

and

$\begin{array}{c}Prob\{X_3=0\text{ in }(t,t+h)|X_3(t)=0\}=\exp \left(\frac{-D_{13}(t)}{\epsilon }h\right)+o(h),\\ Prob\{X_3=1\text{ in }(t,t+h)|X_3(t)=1\}=\exp \left(\frac{-q}{\epsilon }h\right)+o(h).\end{array}$

The time-averaged deterministic model

We employ renewal reward theory (see [15]) to derive a system of ordinary differential equations which replaces (1–6) by a "time-averaged" system in the limit ε → 0. To this end, note first that if D₂ were independent of time, the time average of X₁(t) over "macroscopic" time intervals (i.e., intervals of scale much larger than ε) is $\frac{D_2}{r+D_2}$. The corresponding average of 1 - X₁(t) is then $\frac{r}{r+D_2}$. Renewal reward theory implies that this intuition is mathematically accurate.

Specifically, define a cycle to consist of a period of unoccupied time followed by a period of occupied time. The cycle ends with detachment. The period of unoccupied time is exponentially distributed with mean ε/D₂. Suppose, in the language of renewal reward theory, that no reward is received during this time. The following occupied part of the cycle is exponentially distributed with mean ε/r, and we assume that the reward associated with this period is exactly equal to the amount of occupied time. Then, by renewal reward theory, the long-term average reward (i.e., the proportion of occupied time) is with probability 1 equal to E(R)/E(L) where E(R) is the expected reward during a cycle and E(L) is the expected length of a cycle. In the case under consideration

E(R) = ε/r, E(L) = ε/r + ε/D₂,

so the long-term time average of X₁(t) is D₂/(r + D₂), i.e., lim_ε→0X_1ε(t) = $\frac{D_2}{r+D_2}$ (here, we denote the random variables X _i as X _iε to emphasize the dependence on ε). This time average will hold over any time interval over which D₂ is constant or changes sufficiently slowly. In this time-averaged system Eqns. (1,4) then become

${R^{\prime }}_1=k_{13}\frac{r}{r+D_2}-k_{14}{R}_1\left(11\right)$

${R^{\prime }}_2=k_{13}\frac{D_1}{s+D_1}-k_{14}{R}_2\left(12\right)$

and the remaining equations stay the same. Similarly, Equation (8) becomes

${R^{\prime }}_5=k_{53}\frac{D_{13}}{(q+D_{13})}-k_{54}{R}_5.$

This intuitive argument is not rigorous. As is transparent from the equations for the primary oscillators, all the dependent variables are random variables with time fluctuations at time scale ε. In particular, D₁ and D₂ (and likewise D₃ and D₄) experience stochastic fluctuations in their third derivatives (${R^{\prime }}_1$ experiences random jumps, as does ${P^{″}}_1$, and as does ${D^{‴}}_1$). The integration process involved in the computation of D _i , (i = 1, 2) will average out these fluctuations, so that D _i will indeed vary more slowly than, say, R _i . An argument based on the Arzelà-Ascoli Theorem can be used to translate these observations into a mathematical proof.

To this end we denote by R_1ε, P_1ε, D_1εetc. the solution of (1–6) for some ε > 0 and given initial values R₁(0), P₁(0),..., and denote by R₁, P₁, D₁ etc. the solution of Eqns. (11, 12) ff. for the same initial values.

We prove

Proposition 1 Almost surely for all t > 0,

$\begin{array}{l}\underset{\epsilon \to 0}{\lim }{R}_{1\epsilon }(t)=R_1(t)\\ \underset{\epsilon \to 0}{\lim }{P}_{1\epsilon }(t)=P_1(t)\\ \underset{\epsilon \to 0}{\lim }{D}_{1\epsilon }(t)=D_1(t)\\ \underset{\epsilon \to 0}{\lim }{R}_{2\epsilon }(t)=R_2(t)\end{array}$

etc.

Proof.

Step 1. Consider an arbitrary but fixed time interval [0, T] and let (ε _n ) be a sequence such that ε _n → 0 as n → ∞. For each n we consider a realization, again denoted by R_1εetc., of the initial value problem (1–6) ff. with the given fixed initial data.

The resulting functions $R_{1{\epsilon }_n}$, $P_{1{\epsilon }_n}$, $D_{1{\epsilon }_n}$,... all remain bounded and have (uniformly in ε) bounded first derivatives on [0, T]. By the Arzelà-Ascoli Theorem, there is a convergent subsequence of ε _n , denoted again by ε _n . We denote the limits by ${\tilde{R}}_1$, ${\tilde{P}}_1$,.... What we show next is that these limits are solutions of the deterministic limit equations (11,12) ff.

Step 2. We write ε rather than ε _n to simplify the notation. Observe that

$R_{1\epsilon }(t)=R_1(0)e^{-k_{14}t}+k_{13}{\displaystyle {\int }_0^t(1-X_{1\epsilon })(\tau )e^{k_{14}(\tau -t)}d\tau }$

and

$R_1(t)=R_1(0)e^{-k_{14}t}+k_{13}{\displaystyle {\int }_0^t\frac{r}{r+D_2(\tau )}}{e}^{k_{14}(\tau -t)}d\tau .$

The central step of our proof is showing that ${\tilde{R}}_1$ and ${\tilde{D}}_2$ are also related by (11). This will follow if we can show that for any differentiable function f = f(τ) and any fixed time interval [s, t]

$\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_s^t(1-X_{1\epsilon })(\tau )f(\tau )}d\tau ={\displaystyle {\int }_s^t\frac{r}{r+{\tilde{D}}_2(\tau )}}f(\tau )d\tau .$

To this end consider a partition {s, s +Δ,s + 2Δ,..., s + n Δ = t} of [s, t], where Δ = $\frac{t-s}{n}$. Then

${\displaystyle {\int }_s^t(1-X_{1\epsilon })(\tau )f(\tau )}d\tau ={\displaystyle \sum _{k=0}^{n-1}{\displaystyle {\int }_{s+k\Delta }^{s+(k+1)\Delta }(1-X_{1\epsilon })(\tau )f(\tau )}}d\tau .$

On [s + k Δ, s + (k + 1)Δ] we have

f(τ) = f(s + k Δ) + O(Δ)

so

${\displaystyle {\int }_{s+k\Delta }^{s+(k+1)\Delta }(1-X_{1\epsilon })(\tau )f(\tau )}d\tau =[f(s+k\Delta )+O(\Delta )]{\displaystyle {\int }_{s+k\Delta }^{s+(k+1)\Delta }(1-X_{1\epsilon })(\tau )}d\tau$

Because of the equicontinuity we have uniformly in ε

D_2,ε(τ) = ${\tilde{D}}_2$(s + k Δ) + O(Δ) + μ(ε),

Where μ(ε) → 0 as ε → 0 Hence, by the renewal reward result quoted earlier [15] we have almost surely

$\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{s+k\Delta }^{s+(k+1)\Delta }(1-X_{1\epsilon })(\tau )f(\tau )}d\tau =\frac{r\Delta }{r+{\tilde{D}}_2(s+k\Delta )+O(\Delta )}f(s+k\Delta )+O({\Delta }^2)$

so

$\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_s^t(1-X_{1\epsilon })(\tau )f(\tau )}d\tau ={\displaystyle \sum _{k=0}^{n-1}\frac{r}{r+{\tilde{D}}_2(s+k\Delta )}}\Delta f(s+k\Delta )+O(\Delta ),$

and in the limit Δ → 0 the right hand side converges to

${\displaystyle {\int }_s^t\frac{r}{r+{\tilde{D}}_2(\tau )}f(\tau )}d\tau .$

Step 3. The argument in step 2 and similar (but simpler) reasonings for the other dependent variables show that the R _i , P _i and D _i , i = 1, 2 and the ${\tilde{R}}_i$, ${\tilde{P}}_i$ and ${\tilde{D}}_i$ are both solutions of the same initial value problem. By unique solvability it follows that these solutions are identical, so for example ${\tilde{R}}_1$(t) = R₁(t) for all t. This uniqueness also implies (by a standard argument) that the passage to a subsequence of the ε _n made earlier is not necessary, but that in fact lim_ε→0R _iε (t) = R _i (t) and likewise for all other dependent variables.

This completes the proof.

Remark. This result is only a first step in a possible more complete analysis of the whole process. Specifically, we intend to study the partial differential equations governing the probabilities that the stochastic variables R _i , P _i , D _i assume values in certain ranges, derive the deterministic model given earlier as a set of equations for the first moments of these variables, and proceed to study fluctuations. The nonlinear coupling in our equations makes this a challenging program.

Numerical tests

Here we present some results of simulations performed with the XPPAUT package (see [16, 17]). The chosen parameters are those from Table 1. Figure 1 shows the time course of the proteins P₁ and P₃ for the deterministic model, which oscillate with a period of about 3 hours but differ slightly in their periods. A slight circadian variation is seen; it is much more promiment in Figure 2, where the responses of the protein products of the fifth DNA site are shown; note the time lag of D₅ with respect to D₁₃.

In Figures 3 and 4 the same calculation was done for the stochastic model. This calculation used Gillespie's method [11], where the ε was chosen as 2.8 × 10^-5hrs. The results are essentially identical to the ones for the time-averaged model.

As a control measure we performed some calculations with larger ε, for example ε = 2.8 × 10^-3 hrs and ε = 0.028 hrs. For the former case, especially, the results were close to the time-averaged simulations. For the latter case, deviations from the time-averaged simulations became noticable: the amplitude of the circadian oscillations in D₅ fluctuated stochastically and their period decreased slightly.

Despite these more significant stochastic effects with larger ε, the integrity of the circadian period is remarkably robust in our model with respect to the choice of ε. We demonstrate this by computing Fourier power spectra of D₅ time series generated by simulations with ε = 2.8 × 10^-5 and ε = 2.8 × 10^-2 (see Figures 5 and 6). The former was calculated from a time series of 7447 data points at intervals of 1 minute, representing 124.1 hours of real time. The latter was calculated from a time series of 9920 data points at intervals of 10 minutes, representing 1653.2 hours of real time. We chose to integrate for a longer time in the latter case because the circadian oscillations were less regular. The power spectrum is shown in decibels (decibels = 10 log₁₀(power), where power = |X _i |² for X _i , the i^thfrequency component of the Fourier transform of the time series {x _k }). The frequencies of the primary oscillators show up clearly in the power spectra at close to 8 and 9 cycles per day respectively, and the circadian oscillations are clearly overwhelmingly dominant at close to (but not exactly) 1 cycle per day in both cases. Even after 65 "days" with ε = 0.028, the stochastic oscillator remained in phase with the circadian period; the wave form appeared to persist indefinitely.

Remarks on the frequencies of the primary oscillators

The fundamental idea of our model is that circadian oscillations can easily be achieved via coupling of faster oscillators. We now address the question of whether the primary oscillators could attain circadian periods without need for coupling within reasonable ranges of parameter values based on known biochemistry. To this end we investigated which (if any) intrinsic limitations there are on the periods of the primary oscillators introduced earlier. We first explored (randomly) variations of the growth parameters k₁₃, k₁₅, k₁₇, etc., and the unbinding rates r and s to see how they would affect the periods of the time-averaged single primary oscillator

$\begin{array}{c}\frac{d{R}_1}{dt}=k_{13}\frac{r}{r+D_2}-k_{14}{R}_1\\ \frac{d{P}_1}{dt}=k_{15}{R}_1-k_{16}{P}_1-2k_{17}{P}_1^2+2k_{18}{D}_1\\ \frac{d{D}_1}{dt}=k_{17}{P}_1^2-k_{18}{D}_1\\ \frac{d{R}_2}{dt}=k_{13}\frac{D_1}{s+D_1}-k_{14}{R}_2\\ \frac{d{P}_2}{dt}=k_{25}{R}_2-k_{16}{P}_2-2k_{27}{P}_2^2+2k_{28}{D}_2\\ \frac{d{D}_2}{dt}=k_{27}{P}_2^2-k_{28}{D}_2\end{array}$

Initially we kept the decay parameters k₁₄, k₁₆, k₁₈ fixed and just varied k₁₃. This had a modest effect on the period; the longest which was observed was 3.5 hrs. Random experiments of this nature did not produce periods of circadian length.

For a systematic investigation of the dependence of the periods on the parameters, we then set k₂₅ = k₁₅, k₂₇ = k₁₇, k₂₈ = k₁₈ and linearized the system about its unique positive equilibrium (R_1E, P_1E, D_1E, R_2E, P_2E, D_2E). The linearization yields the 6-by-6 matrix

$A=\left(\begin{array}{cccccc}-k_{14}& 0& 0& 0& 0& \frac{-k_{13}r}{{(r+D_{2E})}^2}\\ k_{15}& -k_{16}-2k_{17}{P}_{1E}& 2k_{18}& 0& 0& 0\\ 0& 2k_{17}{P}_{1E}& -k_{18}& 0& 0& 0\\ 0& 0& \frac{k_{13}s}{{(s+D_{1E})}^2}& -k_{14}& 0& 0\\ 0& 0& 0& k_{15}& -k_{16}-2k_{17}{P}_{2E}& 2k_{18}\\ 0& 0& 0& 0& 2k_{17}{P}_{2E}& -k_{18}\end{array}\right)$

Its eigenvalues satisfy det (A - λI) = 0. This yields the characteristic equation

(k₁₄ + λ)² (k₁₈ + λ)² (k₁₆ + 2k₁₇P_1E+λ) (k₁₆ + 2k₁₇P_2E+λ) + $c_0^2$ = 0.     (14)

Here,

$c_0^2=\frac{4k_{15}^2{k}_{17}^2{k}_{13}^{2}rs}{{(s+D_{1E})}^2{(r+D_{2E})}^2}.$

To identify solutions with longer periods we look for a pair of eigenvalues with positive real part and small imaginary parts. Observe that $c_0^2$ = 0 in (14) produces 6 real and negative eigenvalues (eigenvalues are counted with their multiplicity). If we now increase $c_0^2$, one pair of eigenvalues approaches and eventually crosses the imaginary axis (Hopf bifurcation), producing the oscillations. However, the only way to force the crossing of the imaginary axis at small imaginary value is to move a pair of eigenvalues closer to the imaginary axis to begin with (i.e., when $c_0^2$ = 0).

To achieve this, we first modified the parameter k₁₇ governing the rate of homodimer formation. However, decreasing k₁₇ turns out to increase P_1E, counter-acting attempts to move the crossing pair closer to the real axis.

Finally, the actual rate constant of homodimer decay, k₁₈, is not known, although it is unlikely to be smaller than 1 per hour. Choosing it to be exactly 1 per hour (earlier it was set to 15 per hour) we increased the periods up to 9 hours. Setting k₁₈ this low is probably not reasonable, but given no a priori firm bounds as to how small k₁₈ can actually be (a comment that applies to k₁₄ and k₁₆ as well), no simple predictions on the size of the periods of the primary oscillators can be made.

The following set of parameters produces a wavelength of about 22 hours:

k₁₃ = 1000, k₁₄ = k₁₆ = 1, k₁₅ = 400, k₁₇ = 10^-5, k₁₈ = 0.25, r = 1, s = 9000. Thus almost circadian periods can be obtained, but only by stretching parameters beyond biochemically reasonable values.

Conclusion

We have shown that TTOs in both their stochastic and time-averaged versions produce stable ultradian oscillations for reasonable parameter choices. Although the effect of the stochasticity is to degrade the circadian rhythms as in other models like that of Forger and Peskin [9], these oscillations are nevertheless robust in our model with respect to the scaling parameter governing the dimer-driven stochastic activation or inhibition of the relevant gene sites. Couplings of such TTOs with slight variations in their periods offer a simple mechanism to explain the emergence of circadian rhythms as "beats". This explanation has the added desirable feature of making circadian rhythms readily adaptable to evolutionary pressures.