1 Introduction
Most living organisms, from cyanobacteria to insects, plants and mammals, have developed the capability of generating autonomously sustained oscillations with a period close to 24 h. These oscillations, known as circadian rhythms, are endogenous because they can occur in constant environmental conditions, e.g., constant darkness [1,2]. Experimental studies during the last decade have shed much light on the molecular mechanism of circadian rhythms [3]. Initial studies pertained to the fly Drosophila [4,5] and the fungus Neurospora [3]. Molecular studies of circadian rhythms have since been extended to cyanobacteria, plants and mammals [6,7]. In all cases investigated so far, the molecular mechanism of circadian oscillations relies on the negative autoregulation exerted by a protein on the expression of its gene [3–8]. Thus, in Drosophila, the proteins PER and TIM form a complex that indirectly represses the activation of the per and tim genes, while in Neurospora it is the FRQ protein that represses the expression of its gene frq [3,6]. The situation in mammals resembles that observed in Drosophila, but instead of TIM it is the CRY protein that forms a regulatory complex with a PER protein to inhibit the expression of the per genes [7]. Light can entrain circadian rhythms by inducing degradation of the TIM protein in Drosophila, and expression of the frq and per genes in Neurospora and mammals, respectively [3–7].
A number of mathematical models for circadian rhythms have been proposed [9–15] on the basis of these experimental observations. These models are of a deterministic nature and take the form of a system of coupled ordinary differential equations. The models predict that in a certain range of parameter values the genetic control network undergoes sustained oscillations of the limit cycle type corresponding to the circadian rhythm, whereas outside this range the gene network operates in a stable steady state. The question arises as to whether deterministic models are always appropriate for the description of circadian clocks [16]. Indeed, the number of molecules involved in the regulatory mechanism producing circadian rhythms at the cellular level may well be reduced. This number could vary from a few thousands down to hundreds and even a few tens of protein or messenger RNA molecules in each rhythm-producing cell. At such low concentrations it is more appropriate to resort to a stochastic approach to study the molecular bases of the oscillatory phenomenon.
In a previous publication [17], we compared the stochastic and deterministic versions of a core molecular model for circadian oscillations based on the negative regulation exerted by a protein on the expression of its gene. Stochastic simulations were performed by means of the Gillespie algorithm [18,19] after decomposing the deterministic model into elementary steps. We studied the effect of molecular noise by assessing the robustness of circadian oscillations as a function of the number of interacting molecules. We showed that robust circadian rhythmicity could already occur when the maximum numbers of mRNA and clock protein molecules are in the tens and hundreds, respectively. Cooperativity of repression and periodic forcing by light-dark cycles enhance the robustness of circadian oscillations. In subsequent work [20], we compared two stochastic versions of this core model, one fully developed into elementary steps, and the other non-developed. We showed that stochastic treatment of these two versions of the model for circadian rhythms yields similar results.
The purpose of the present paper is to extend our comparison of deterministic and stochastic models for circadian oscillations, by considering a more detailed model for the Drosophila circadian clock incorporating the formation of a complex between the PER and TIM proteins. Although this model does not take into account other proteins such as CLOCK and CYC involved in the circadian oscillatory mechanism, it is nevertheless more realistic than the core model and accounts for a larger number of experimental observations. Moreover, this extended model predicts the possibility of autonomous chaotic behaviour [21]. Such a property will allow us to assess the effect of molecular noise not only on periodic but also on chaotic oscillations.
We first present in Section 2 the deterministic and stochastic versions of the molecular model for circadian rhythms. In the stochastic version we introduce molecular noise without decomposing the deterministic mechanism into detailed reaction steps. The results of stochastic simulations performed by means of the Gillespie method [18,19] are presented in Section 3 for the case of periodic behaviour. We assess the role of fluctuations by determining the effect of the number of mRNA and protein molecules on circadian rhythmicity. In Section 4 we examine how the proximity from a bifurcation point influences the robustness of the oscillations with respect to molecular noise. How such a noise affects chaotic behaviour is examined in Section 5. Section 6 is devoted to a comparative discussion of deterministic versus stochastic simulations.
2 Deterministic and stochastic versions of the molecular model for circadian oscillations
2.1 Deterministic model for circadian oscillations
We consider a ten-variable model previously proposed for circadian oscillations of the PER and TIM proteins and of per and tim mRNAs in Drosophila [9,10]. The model, schematised in Fig. 1, is based on the negative feedback exerted by the complex between the nuclear PER and TIM proteins on the expression of their genes. For each of these proteins the gene is first transcribed in the nucleus into messenger RNA (mRNA). The latter is transported into the cytosol where it is degraded and translated into the protein P0 (T0). The protein PER (TIM) undergoes multiple phosphorylation, from P0 into P1 (T0 into T1) and from P1 into P2 (T1 into T2). These modifications, catalysed by a protein kinase, are reverted by a phosphatase. The fully phosphorylated form of the proteins is marked up for degradation and forms a complex (C), which is transported into the nucleus in a reversible manner. The nuclear form of the PER–TIM complex (CN) represses the transcription of the per and tim genes. Recent experiments indicate that repression is in fact of indirect nature: the CLOCK and CYC proteins promote the expression of the per and tim genes and are prevented from exerting this activation when forming a complex with PER and TIM [3–7]. In the model, the variables are the concentrations of the mRNAs (MP and MT), of the various forms of the PER and TIM proteins (P0, P1, P2, T0, T1, T2), and of the cytosolic (C) and nuclear (CN) forms of the PER–TIM complex. The temporal evolution of these concentration variables is governed by a system of 10 kinetic equations that are listed in Appendix A (see [9,10] for further details).
The deterministic PER–TIM model accounts for the occurrence of sustained oscillations in continuous darkness. When taking into account the control exerted by light on the maximum protein degradation rate vdT, the model also accounts for entrainment of circadian oscillations by light-dark cycles and for their phase shifting by pulses of light. Circadian oscillations have also been obtained in more detailed models based on indirect repression and involving additional clock gene products such as CLOCK and CYC ([14,15], J.-C. Leloup and A. Goldbeter, submitted for publication).
As indicated above, besides periodic oscillations, the model of Fig. 1 is also capable of generating autonomous chaos in conditions corresponding to continuous darkness. Although this behaviour is probably not of physiological significance [21], it provides us with the rare opportunity of testing the effect of molecular noise on chaotic behaviour in a realistic model based on genetic regulation.
2.2 Stochastic version of the model for circadian oscillations
To assess the effect of molecular noise, we describe the reaction steps as stochastic birth and death processes [22]. Numerical simulations of the temporal evolution of the genetic control system are performed by means of the Gillespie method [17,18]. Besides other approaches [23–25], this method has been used to determine the dynamics of chemical [24,25] biochemical [26] or genetic systems [27] in the presence of molecular noise. The Gillespie method associates a probability with each reaction step; at each time step the algorithm randomly determines the reaction that takes place according to its probability, as well as the time interval to the next reaction step. The numbers of molecules of the different reacting species as well as the probabilities are updated at each time step. In this approach [18,19], a parameter denoted controls the number of molecules present in the system. Using the Gillespie method we performed stochastic simulations of the model described in section 2.1.
Our previous analysis of a core molecular model for circadian rhythms showed [17] that similar results are obtained when decomposing the deterministic model into elementary steps (developed model) or when the deterministic model is not decomposed into such steps and non-linear kinetic functions are simply included in the probabilities associated with the global reaction steps (non-developed model). Because it is much easier to handle, we shall therefore restrict our stochastic analysis to such a non-developed version of the model schematised in Fig. 1.
The non-linear terms appearing in the kinetic equations (A.1) listed in Appendix A represent compact kinetic expressions obtained after application of quasi-steady-state hypotheses on enzyme–substrate or gene–repressor complexes. The resulting expressions are of the Michaelis–Menten type for enzyme reaction rates, or of the Hill type for cooperative binding of the repressor (here, the nuclear PER–TIM complex) to the gene promoter. We attribute to each linear or non-linear term of the kinetic equations a probability of occurrence of the corresponding reaction. These reactions and their associated probability are listed in Table 1 in Appendix A. Thus reaction (1) corresponds to the transcription of the per gene into per mRNA, MP; the occurrence of this reaction with a probability w1 results in increasing by 1 the number of molecules of MP. Reaction (4) results in increasing by one the number of P1 molecules and decreasing by one the number of molecules of P0.
Stochastic model for circadian rhythms, corresponding to the mechanism schematised in Fig. 1
Reaction number | Reaction step | Probability of reaction |
1 | ||
2 | ||
3 | w3=ksP×MP | |
4 | ||
5 | ||
6 | ||
7 | ||
8 | ||
9 | w9=k4×C | |
10 | ||
11 | ||
12 | ||
13 | w13=ksT×MT | |
14 | ||
15 | ||
16 | ||
17 | ||
18 | ||
19 | w19=k1×C | |
20 | w20=k2×CN | |
21 | w21=kd×MP | |
22 | w22=kd×P0 | |
23 | w23=kd×P1 | |
24 | w24=kd×P2 | |
25 | w25=kd×MT | |
26 | w26=kd×T0 | |
27 | w27=kd×T1 | |
28 | w28=kd×T2 | |
29 | w29=kdC×C | |
30 | w30=kdN×CN |
In contrast to the treatment presented in our previous work [17], here we do not decompose the binding of the repressor to the gene promoter into successive elementary steps, and rather retain the Hill function description for cooperative repression by the nuclear PER–TIM complex CN (steps (1) and (11)). A similar global approach is taken for describing degradation of per and tim mRNAs (reactions (2) and (12)), translation of mRNAs into proteins P0 and T0 (reactions (3) and (13)), phosphorylation of P0 and T0 into P1 and T1 (reactions (4) and (14)) and of P1 and T1 into P2 and T2 (reactions (6) and (16)), as well as dephosphorylation of P1 and T1 into P0 and T0 (reactions 5 and 15) and of P2 and T2 into P1 and T1 (reactions (7) and (17)), enzymatic degradation of P2 and T2 (reactions (10) and (18)), reversible formation of the complex C (reactions 8 and 9), and reversible transport of complex C into and out of the nucleus (reactions (19) and (20)). Steps (21)–(30) relate to non-specific degradation of the various mRNA or protein species. Reactions (2), (4)–(7), (10), (12), and (14)–(18) are of the Michaelian type; reaction (8) is of bimolecular nature, while reactions (3), (9), (13) and (19)–(30) correspond to linear kinetics.
3 Effect of molecular noise on periodic oscillations
Before dealing with the effect of molecular noise, let us recall the predictions of the deterministic model governed by equations (A.1). Typical circadian oscillations predicted by this model are shown in Fig. 2A (middle panel). These oscillations correspond to the evolution toward a limit cycle (Fig. 2A, left panel) shown here as a projection onto the plane formed by the concentrations of per mRNA (MP) and nuclear PER–TIM complex (CN). Because the behaviour is periodic and deterministic (i.e. oscillations occur in the absence of molecular noise), the histogram of periods yields a single line corresponding to a circadian period close to 24 h (Fig. 2A, right panel).
Turning to the effect of molecular noise, we now consider the dynamic behaviour of the stochastic version of the PER–TIM model for circadian rhythms. Shown in panels B–D in Fig. 2 are the limit cycles (left panels), sustained oscillations of per mRNA (middle panels), and histograms of periods (right panels) obtained by stochastic simulations for (B), 100 (C) and 10 (D). For these values of , the numbers of molecules of nuclear PER–TIM complex and per mRNA vary in the range 500–2500 and 0–3000, 50–300 and 0–300, and 0–60 and 0–50, respectively.
The data presented in Fig. 2 show that the circadian oscillatory behaviour predicted by the deterministic model is recovered when using the stochastic model for circadian rhythms. The mere effect of molecular noise is to increase the effective thickness of the limit cycle. The results further indicate that robust circadian oscillations are still produced by the stochastic model when the maximum numbers of mRNA and protein molecules are in the order of hundreds. It is only when these numbers decrease down to a few tens that noise begins to overcome rhythmicity, even though oscillations still subsist (Fig. 2D); the period histogram is then much wider but still presents a maximum close to a circadian value.
Another way to illustrate the effect of molecular noise is to draw a histogram of frequencies of passage through the different points in the phase plane. In the deterministic case, such a plot would yield the same frequency of passage through all points of the limit cycle, because these are all visited the same number of times. In the presence of noise, some regions of the phase space are visited more often than others where dispersion is strong. This results in the occurrence of peaks in the frequency of passage, as shown in Fig. 3A. The bottom part of Fig. 3A shows a contour plot in the plane (per mRNA molecules versus nuclear PER–TIM complexes). The contour plot is obtained by projecting onto this plane the intersections of the histogram with 50 parallel planes, corresponding to equidistant frequencies. For the sake of clarity the contour plot is rotated by 90° and shown enlarged in Fig. 3B. The most frequently visited regions in the phase plane correspond to the decrease in nuclear PER–TIM levels and rise in per mRNA. Largest dispersion occurs when mRNA decreases and PER–TIM levels rise.
4 Influence of the vicinity of a bifurcation point
The question arises as to whether the proximity from a bifurcation point may influence the robustness of circadian oscillations with respect to molecular noise. To address this issue, we first construct in Fig. 4 a bifurcation diagram showing the onset of limit cycle oscillations as a function of the maximum rate of TIM degradation, vdT, for the deterministic model governed by equations (A.1). The state of the system is represented in Fig. 4 by a single variable, the concentration of per mRNA. Below , the system evolves toward a stable steady state. Above this critical bifurcation value, limit cycle oscillations occur. The diagram shows the envelope of the oscillations: the upper and lower branches yield the maximum and minimum values of per mRNA concentration as a function of vdT in the course of oscillations.
The results of stochastic simulations performed with increasing values of parameter vdT are shown in Fig. 5 in the form of dynamics in the phase plane (right panels) or time series (left panels). The steady state or limit cycle predicted by the deterministic model for the corresponding parameter values is shown in the phase plane by a white dot or a white curve, respectively. The four panels in Fig. 5 correspond to the four vdT values indicated by dashed vertical lines in Fig. 4. For , the deterministic system evolves to a stable steady state far from the bifurcation point; stochastic simulations show low-amplitude fluctuations around the deterministic steady state (Fig. 5A). For , the deterministic system evolves to a stable steady state close to the bifurcation point; stochastic simulations show fluctuations of larger amplitude around the deterministic steady state (Fig. 5B). For , the deterministic system undergoes limit cycle oscillations of reduced amplitude just beyond the bifurcation point; stochastic simulations show small-amplitude, noisy oscillations around the deterministic limit cycle (Fig. 5C), which resemble the fluctuations shown in Fig. 5B. Finally, for , the deterministic system undergoes limit cycle oscillations of large amplitude far from the bifurcation point; stochastic simulations show large-amplitude, relatively less noisy oscillations around the deterministic limit cycle (Fig. 5D).
The data in Fig. 5 indicate that stochastic simulations allow us to recover the dynamics predicted by the deterministic model. Below a critical parameter value, the system displays low-amplitude fluctuations around a stable steady state, while above this value sustained oscillations occur. Moreover, the effect of fluctuations is enhanced by the proximity from the bifurcation point.
5 Effect of molecular noise on chaotic behaviour
We have previously reported the occurrence of autonomous chaos in the deterministic model of Fig. 1 [21]. Thus, for some parameter values, in conditions corresponding to continuous darkness, sustained aperiodic oscillations occur in this model, which correspond to the evolution toward a strange attractor in the phase space (see Fig. 6A). This phenomenon provides us with the rare opportunity to assess the effect of molecular noise on chaotic behaviour in a realistic biochemical model. Shown in panels B and C of Fig. 6 for and 100, respectively, are the results of stochastic simulations performed for parameter values corresponding to those producing chaos in the deterministic model in Fig. 6A.
The results indicate that chaos persists in the presence of noise, but the structure of the strange attractor begins to be blurred when the number of molecules decreases and the amplitude of fluctuations rises. Nevertheless the small curler-like substructure that characterises the strange attractor in the PER–TIM model (Fig. 6A, left panel) is still visible in the attractors obtained by stochastic simulations (left panels in Fig. 6B and C).
Because the PER–TIM model can produce periodic as well as chaotic oscillations, we can use this model to compare the effect of molecular noise on the two types of dynamic behaviour. A convenient tool for such a comparison is provided by Poincaré sections. For the case of periodic oscillations, the trajectory in the (MP, MT, CN) plane takes the form of a limit cycle. Intersection of this trajectory with a plane corresponding to a given value of MT generally yields two points, one for which MT is on the rise and the other for which MT is on the decline. Shown in Fig. 7A is the point intersection obtained for the deterministic model when MT passes the value 1.5 nM upward. Panels B and C in Fig. 7 show the Poincaré sections obtained by stochastic simulations with and 100, respectively, for the corresponding value of the number of MT molecules (i.e., 1.5 ). Instead of a single point, we obtain a cloud of points surrounding the deterministic Poincaré section; the smaller the number of molecules, the more scattered the cloud.
In the case of chaos, instead of a single point the intersection in the deterministic model takes the form of an open continuous curve (Fig. 7D). Stochastic simulations produce a cloudy version of this curve; the scattering of the points in this cloud increases as the number of molecules decreases (compare panels E and F in Fig. 7). Interestingly, in Fig. 7F, we can observe a denser part of the cloud in the left upper part of the figure; this dense region corresponds to noise-induced visits to the curler-like sub-region of the strange attractor (see Fig. 6, left panels).
6 Discussion
Circadian rhythms originate from the negative autoregulation of gene expression. Models based on such regulatory processes have been proposed for circadian rhythms in Drosophila and Neurospora. These models are of a deterministic nature. When molecular noise becomes significant, however, the question arises as to whether genetic control mechanisms can give rise to coherent circadian oscillations at the cellular level. In the presence of reduced numbers of molecules of the mRNA and protein species involved in the circadian clock mechanism, stochastic simulations are needed to address this issue. In a previous publication [17], we studied the stochastic version of a core model for circadian oscillations based on negative autoregulation of a clock gene by its protein product. For stochastic simulations, this deterministic five-variable model was first decomposed into a series of 30 elementary reaction steps. Stochastic simulations performed by means of the Gillespie method [18,19] indicated that robust circadian oscillations, comparable to those predicted by the deterministic approach, can occur in this developed stochastic model even when the maximum numbers of mRNA and clock protein molecules are in the tens and hundreds, respectively. The robustness of circadian rhythms, quantified by the period distribution and the half-time of the autocorrelation function [17], is enhanced when the cooperativity of the repression process increases and when the numbers of mRNA and protein molecules involved in the oscillatory mechanism become larger. Formation of complexes between regulatory proteins is another stabilising factor (D. Forger and C. Peskin, personal communication). Stochastic simulations indicate, moreover, that forcing by a light-dark cycle stabilises the phase of the oscillations [17].
In a follow up study [20] we compared two stochastic versions of the core model for circadian oscillations. In one version (the developed stochastic model), described above, all reactions were decomposed into elementary steps, while in the second version (non-developed model), the reactions were not decomposed into elementary steps and non-linear rate expressions were simply included in the probability of occurrence of each reaction. A comparison of the two stochastic versions of the core model indicates that they both yield similar results; the two approaches are equivalent and corroborate the predictions of the deterministic model. It is only when the maximum numbers of mRNA and protein molecules are reduced to very low values, of the order of a few tens, that molecular noise begins to prevent the emergence of coherent circadian oscillations. These numerical results are corroborated by a recent analytical study [28].
Here we have extended the comparison of deterministic and stochastic models for circadian rhythms by considering a more extended molecular model proposed for circadian oscillations of the PER and TIM proteins and their mRNAs in Drosophila. The ten-variable model is based on the negative autoregulation of the per and tim genes by a complex between the PER and TIM proteins. Although additional proteins are known to be at work in the molecular mechanism of the circadian clock in Drosophila and mammals, the PER–TIM model appears as more realistic than the core mechanism that we previously considered for deterministic [9] and stochastic simulations [17]. Because the non-developed and developed stochastic versions of the core model for circadian rhythms gave similar results, we focused here on the comparison of the deterministic PER–TIM model with its non-developed stochastic version, which is easier to handle than the developed one. An advantage of the PER–TIM model is that, in addition to periodic oscillations, it also admits, for appropriate parameter values, chaotic behaviour. This allowed us to assess the effect of molecular noise both on periodic and chaotic oscillations.
Stochastic simulations of the PER–TIM model corroborate the conclusions reached for the simpler core mechanism for circadian rhythms. Sustained circadian oscillations predicted by the deterministic model are recovered by stochastic simulations. Robust circadian oscillations closely related to those obtained with the deterministic model can already occur when the maximum numbers of protein and mRNA molecules are in the hundreds or tens, respectively (the mean numbers of these molecules in the course of oscillations are smaller, because the levels of both protein and mRNA can go down to very low values at the trough of the oscillations). Robustness increases with the numbers of molecules involved, as reflected by the period distribution (see the period histograms in Fig. 2). It is only at very low numbers of molecules of protein and mRNA (reached, for example, for ) that noise begins to obliterate circadian rhythmicity [17]. In both the deterministic and stochastic versions of the PER–TIM model, critical parameter values separate the evolution toward a stable steady state from the evolution toward sustained oscillations. The effect of noise is merely to induce fluctuations around the stable steady state or around the stable limit cycle predicted by the deterministic model. As shown here (see Fig. 5), and elsewhere on the core model [20], the effect of molecular noise is enhanced in the vicinity of a bifurcation point.
A previous study of models based on negative autoregulation of gene expression reported a lack of robustness of circadian oscillations with respect to molecular noise [16]. The difference with respect to the results presented here and elsewhere [17,20] is likely due to the lower values considered for bimolecular rate constants characterising the association of the repressor to the gene promoter [17,29]. When these parameter values remain too small, the analysis of the deterministic model indicates that the steady state becomes stable but displays the property of excitability: in the stochastic model, instead of sustained oscillations, fluctuations give then rise to irregular, large-amplitude excursions away from the stable steady state.
Autonomous chaos has previously been described in the deterministic PER–TIM model for circadian rhythms in Drosophila [21]. As for periodic circadian oscillations, chaotic behaviour is recovered when stochastic simulations are performed for parameter values corresponding to chaos in the deterministic model. Beyond a noisy appearance due to fluctuations in the presence of reduced numbers of mRNA and protein molecules, the structure of the strange attractor remains discernible (Fig. 6), in agreement with results obtained in models for chemical oscillations [30,31]. The use of Poincaré sections (Fig. 7) allowed us to place on firmer ground the comparison of periodic and chaotic behaviour in the presence of molecular noise. The results lead to a clear distinction between noisy limit cycle oscillations and chaotic oscillations subjected to noise. The agreement between the stochastic and deterministic approaches therefore holds for periodic oscillations as well as for chaotic behaviour.
Acknowledgements
This work was supported by Grant No. 3.4607.99 from the Fonds de la recherche scientifique médicale (FRSM, Belgium). J.-C.L. is Chargé de recherches du FNRS (Belgium). Support from the Fondation David & Alice Van Buuren to D.G. is also acknowledged. We wish to thank F. Baras, G. Dupont, D. Forger, and P. Gaspard for fruitful discussions.
Appendix A
A.0.1 Kinetic equations of the deterministic model for circadian oscillations
The time evolution of the 10 concentration variables in the model of Fig. 1 is governed by the following system of ordinary differential equations: (A.1)
A.0.2 Stochastic version of the model for circadian oscillations
For stochastic simulations, the model schematised in Fig. 1 and governed by the deterministic system of equations (A.1) is presented as a sequence of reaction steps in Table 1. Steps (21) to (30) refer to non-specific degradation reactions, of relatively reduced importance, which are not indicated in Fig. 1. The second column in Table 1 lists the sequence of reactions (see section 2.2). The probability of each reaction is given in the third column. Parameter values are the same as those listed above for the deterministic model, but units of nM have to be replaced by numbers of molecules, after multiplication of the parameter by , if need be, as indicated in the last column of Table 1.