New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Stochastic simulation of the mammalian circadian clock

Contributed by Charles S. Peskin, November 14, 2004
Abstract
Circadian (nearly 24h) clocks are remarkably accurate at timing biological events despite the randomness of their biochemical reactions. Here we examine the causes of their immunity to molecular noise in the context of a detailed stochastic mathematical model of the mammalian circadian clock. This stochastic model is a direct generalization of the deterministic mammalian circadian clock model previously developed. A feature of that model is that it completely specifies all molecular reactions, leaving no ambiguity in the formulation of a stochastic version of the model. With parameters based on experimental data concerning clock protein concentrations within a cell, we find accurate circadian rhythms in our model only when promoter interaction occurs on the time scale of seconds. As the model is scaled up by proportionally increasing the numbers of molecules of all species and the reaction rates with the promoter, the observed variability scales as 1/n ^{0.5}, where n is the number of molecules of any species. Our results show that gene duplication increases robustness by providing more promoters with which the transcription factors of the model can interact. Although PER2 mutants were not rhythmic in the deterministic version of this model, they are rhythmic in the stochastic version.
The unicellular organism's timing of daily (circadian) biological events like luminescence (1) and O_{2} consumption (2) can be attributed to an intracellular clock consisting of oscillating protein feedback loops. Higher organisms can time sleep, hormone release, and other biological processes by means of a population of cells, each of which seems to contain a biochemical clock similar in basic design to that found in unicellular organisms (3). Isolated cellular circadian clocks can time events with an error of less than ±10% of the 24h circadian day (2, 4) and might be able to time events with an error of <1% of the circadian day (1). Because chemical reactions in cells involve finite (and often low) numbers of molecules (5), individual molecular interactions can be important. One cause of circadian mistiming is molecular noise (the fact that interactions between individual molecules are stochastic). The accuracy of circadian clocks can be essential for survival, and there is likely a strong evolutionary selection to overcome molecular noise.
This article is concerned with the mechanisms by which biochemical circadian clocks maintain high accuracy with low numbers of molecules despite the stochasticity of individual molecular interactions (molecular noise). Much recent attention has focused on model predictions of the accuracy of circadian rhythms in the presence of molecular noise, and the predictions of these models have been somewhat conflicting (6, 7). Although there is a widely accepted scheme for stochastic simulation of molecular processes, the Gillespie method (8), circadian clock models often are not detailed enough to describe individual molecular interactions. Instead, several reaction steps are typically summarized by Michaelis–Menton kinetics or by terms that involve Hill coefficients, and there is ambiguity in formulating a stochastic version of such reaction schemes, particularly because certain reaction rates required in the stochastic simulation are not specified in the corresponding deterministic models (9). In addition, there have not been experimental data on the numbers of molecules of key clock proteins.
Previously, we developed a detailed model of the circadian clock within the mammalian pacemaker (suprachiasmatic nucleus) that is based on individual molecular reactions (11). This model accounts for many more biochemical species than other models and has features (transcription, translation, protein feedback, dimerization, and the sequestration of proteins, i.e., localization) that appear to be used in the circadian clocks of many organisms. Mammalian data, which we use to estimate the number of molecules of key chemical species of our model in the cell, are now available. This model and these data put us in a unique situation to study the effects of molecular noise on circadian clocks and to study the first stochastic mammalian circadian clock model. In particular, we are interested in (i) predictions on the accuracy of the clock within individual cells, (ii) the differences in behavior between our deterministic and stochastic models, and (iii) what design principles contribute to immunity to molecular noise.
Simulation Methods
Our mammalian model consists of five key proteins, PER1, PER2, CRY1, CRY2, and CKI. CKI phosphorylates PER1 and PER2 to regulate their nuclear entry and stability. CRY1 or CRY2 are transported by PER1 and PER2 into the nucleus of the cell, where they can bind to the control regions of the PER1, PER2, CRY1, and CRY2 genes to affect transcription. Reverbα is not included, because we previously found that its effect is negligible (10). Thus, our model consists of a biochemical feedback loop where PER1, PER2, CRY1, and CRY2 oscillate with an ≈24h period. A more detailed description of the model can be found in refs. 10 and 11. In particular, a complete list of reactions and parameters, suitable for a Gillespie simulation, can be found in ref. 11.
Stochastic simulations of our model used the Gillespie method (8), but in an efficient way. Because the number of reaction classes (a reaction that can occur in essentially the same way in many different molecular species; see ref. 9) in our model is much less than the total number of individual reactions, we first determined in which reaction class the next individual reaction would be, and then we chose randomly which specific reaction within that class would occur. On one hand, this method could increase computation cost, because it requires the generation of two random numbers (rather than one in the original Gillespie method) to choose which reaction occurs (one random number is needed to determine the reaction class, and another is needed to distinguish between reactions in this class). On the other hand, this method causes fewer comparisons to be made when determining which reaction will occur (whole classes of reactions can be ruled out with a single comparison), and fewer summations once the reactions occur (we only need to resum the rates in reaction classes that were changed). An informal comparison of this method and the original Gillespie method showed a significant decrease in computation time. Our method is similar to that of Gibson and Bruck (12), which is somewhat harder to implement but may be more computationally efficient.
Because many circadian gene mutations show semidominant behavior (13), we report simulations where both copies of each gene are active, unless otherwise stated.
Our deterministic model already simulated correctly the relative values of the maximum and minimum concentrations of the rhythms of key proteins in the circadian clock in liver cells (relative to the peak of the CRY1 rhythm) found by Lee et al. (14). From the number of cells used in these studies, one can estimate that there are ≈5,000 molecules of CRY1 at its rhythm peak (Chogoon Lee and David Weaver, personal communication). It is important to note that CRY1 is the most prevalent protein in the clock, and many other protein species exist in much smaller numbers. To convert from a deterministic model (which describes concentrations) to a stochastic model (which describes numbers of molecules), we needed to know the volume of a cell. This volume was chosen so that there are ≈5,000 molecules of CRY1 at its rhythm peak. We adjusted this volume when we wanted to study the effect of the overall number of molecules on the variability of the circadian clock.
Initial Simulations and a Reassessment of the Promoter Binding/Unbinding Rates
In our initial stochastic simulations, we found rhythms that varied greatly in amplitude and period from cycle to cycle (Fig. 1 Top, blue curve). In searching for a cause of this irregularity, we found that the average time between the activations of each PER gene was ≈39 h. Indeed, although there can be hundreds or thousands of molecules of each of the other molecular species, there are typically only a few copies of any particular gene. Therefore, compensatory mechanisms are particularly needed to reduce the molecular noise associated with transcription. In our model, a molecule of either of the repressors CRY1 or CRY2 is typically bound to any of the five possible binding sites (specifically Eboxes with CLK:BMAL1) on the PER genes. For transcription to occur, there must be no bound repressor to any of the sites, and this is a somewhat rare event. Our deterministic model did not reveal this because, being macroscopic, it could only consider the average behavior of an ensemble of identical promoters. Thus, although the average activation of many promoters oscillates with a 24h time scale, each individual promoter may be turned on less frequently. Although the resulting irregular rhythms could become more regular when many cells are coupled, we here consider a mechanism that can reduce this irregularity within an isolated cell.
When we increased the promoter binding and unbinding rates by a factor of 100 in comparison to those used above, the average time between activations of a typical PER gene was reduced from 39 h to 1.22 h. This increase gave much more robust oscillations (Fig. 1 Top, black curve). Rapid promoter interactions have been seen experimentally (15, 16) and are also assumed in most circadian clock models through Hill terms (17). Changing these parameters made almost no difference in the behavior of our deterministic model (Fig. 1 Middle). The stochastic model with increased promoter binding and unbinding showed relatively good agreement with the deterministic model (Fig. 1 Bottom). From now on, we will use these faster promoter binding/unbinding rates.
It is remarkable that events involving the promoter must occur on such a fast time scale to achieve robust rhythms on a 24h time scale. Although no data are available on the binding rates of circadian transcription factors, other genes are turned on or off on the time scale of minutes or faster (15, 16). This short time scale is not at all unreasonable for promoter binding, especially because transcription factors can be localized to the nucleus of the cell or possibly even to subnuclear compartments directly surrounding a promoter, which allows for rapid binding by increasing the effective concentration. Recent experimental and theoretical work demonstrated that rapid interactions with a promoter can greatly reduce stochasticity in eukaryotic genetic networks (18–20). Rapid interaction might also be facilitated by CLK:BMAL1 or by other transcription factors that have a high binding affinity (21–23). The unbinding of protein from promoters often occurs on a much faster time scale than the binding of protein to promoters (16). Bacteria may not be as efficient in sequestering molecules as eukaryotes and therefore might be able to generate only noisy rhythms. This fact may partially explain why many bacteria (although not all; see ref. 3) do not contain circadian clocks (24).
The Stochastic Treatment Can Approach the Deterministic Treatment in a Characteristic Way
Next, we explored in what limits the stochastic treatment of our model would converge to its deterministic treatment. By increasing the cell volume while keeping concentrations the same, we created a situation in which more molecules would be present and more reactions would occur. One might then expect the stochastic treatment to approach the deterministic treatment, but this is not the case, for the following reason: There is a fixed number of copies of a particular gene within a cell (typically two), and this number will not increase with increasing cell volume. The fixed number of copies of a particular gene will result in residual stochasticity that cannot be diminished by increasing the cell volume. To see this, imagine the limit of large volume, where all reactions can be modeled deterministically except those involving the promoter. When a binding site on the promoter is unoccupied, there still will be a random time for one of the transcription factors to bind. Moreover, once a transcription factor is bound, its unbinding also will occur stochastically. Because it is not realistic to add more genes, we can instead increase the binding and unbinding rates with the promoters. Although the individual binding and unbinding events will still be stochastic, these events will occur on such a fast time scale that the randomness of the individual events will average out and will not affect the slower time scale of the circadian rhythm. Thus, by both increasing cell volume and promoter binding–unbinding rates, the overall stochasticity of the model can be reduced to arbitrarily small levels.
The actual stochasticity of the model is difficult to calculate analytically because of its many nonlinear interactions and feedback; however, we can estimate how the stochasticity diminishes with increasing cell volume and promoter binding/unbinding rates over a small time interval. Choose some time interval that is small enough that we can consider all reaction rates, with the exception of those involving the promoter binding and unbinding, as effectively constant. Next choose any reaction besides binding and unbinding to the promoter. Over this time interval, the number of times this reaction will occur, with the Gillespie method, will scale with increasing cell volume and obey Poisson statistics. Thus, the mean (over many trials) number of times the reaction will occur in this interval will scale with the cell volume, and the standard deviation divided by the mean will decrease in proportion to 1/(cell volume)^{0.5}. Next consider the binding and unbinding to the promoter. On a time scale comparable to this short time interval, the principle of averaging tells us that the overall system's dynamics depend on the amount of time the gene in “on” rather than the individual fluctuations. Because the promoter is continually being turned on and off, the total amount of time the promoter is “on” is the sum of the individual durations of time the promoter is “on” each time it is turned on. These “on” times are random variables that are independent and identically distributed on this time scale. Therefore, because the binding and unbinding rates are increased by a factor n, we will, on average, sum over n times as many individual “on” times (each of which, on average, has a duration that scales as 1/n), and the standard deviation of the sum (which is independent of n) of these “on” times should decrease as 1/n ^{0.5}. It is interesting to note that Ehret and Trucco (24) used a similar argument in the Chronon model of the molecular circadian clock 35 years ago. Accordingly, if we scale the cell volume and the promoter binding and unbinding rates by n, we should see stochasticity of the model decrease as 1/n ^{0.5}.
The periodtoperiod variability of our unscaled model had a standard deviation of 3.68 h. This standard deviation is larger than what has been found experimentally (4), but these experimental results were from cells that were likely coupled. Liu et al. (31) postulated that this coupling led to more accurate timekeeping. As suggested above, simulations showed decreased variability as more molecules and promoter events were present, with a 1/n ^{0.5} scaling (Fig. 2).
Simulation of Circadian Clock Mutations
Considering the relatively low numbers of molecules of certain chemical species, it is surprising that such robust rhythms are seen. Perhaps this has something to do with the structure of the mammalian circadian clock. To test this hypothesis, we studied the stochastic behavior of our model when each of the key proteins was removed (Fig. 3 Upper) or when only one copy of both of the PER or CRY genes was active (but producing mRNA at twice the normal rate) (Fig. 3 Upper).
One might initially think that removing any of the clock proteins would cause fewer molecules to be produced and more variable rhythms. However, because the proteins are within negative feedback loops, if one PER protein or one CRY protein were removed, more of the other PER or CRY protein respectively would be produced. This fact explains why CRY1 or CRY2 mutants (either CRY1 or CRY2 was removed) were not less precise than our wildtype model. Interesting, the CRY mutants were more precise than the wildtype model. When different nonlinear oscillations are coupled, they often have a lower amplitude than when similar oscillations are coupled (25). Our deterministic model shows larger amplitude rhythms with homogeneous CRY (CRY1 or CRY2 mutants) than with heterogeneous CRY (wild type). As discussed below, large amplitude rhythms can lead to more robust rhythms.
It is rather remarkable that the PER2 mutant (i.e., PER2 removed) was rhythmic (Fig. 3 Lower) in our stochastic model (although every so often the amplitude of the rhythms briefly decreased to undetectable levels) even though the deterministic PER2 mutant was not rhythmic (11). Although there may be an inherent drive to dampen oscillations in the stochastic model (as in the deterministic model), this damping appears to be offset by noise in the stochastic case. A similar result can be seen with a version of our model where just three Eboxes are active on PER promoters; although the deterministic model in this case is not rhythmic, the corresponding stochastic model does oscillate (data not shown).
When the number of copies of each PER, or CRY gene was reduced from two to one, the rhythms became more variable, because fewer gene activations occurred. A much stronger effect was seen in the PER onecopy mutant, because the CRY gene is activated so frequently that added benefit of having two copies is minimal. This fact explains why the PER1 and PER2 mutant models were more noisy than the wildtype models: four copies of the PER genes are better than two.
Discussion
The starting point of this study is a deterministic model whose stochastic counterpart is unambiguous. Instead of considering relatively simple models to illustrate general principles, we employed a biochemically detailed, validated, and distinctly mammalian model to reverse engineer the design principles of in vivo circadian clocks. We also made use of an experimental estimate of the number of molecules of each protein species (a key parameter in stochastic stimulation). Very recent experimental preparations now make it possible to assess the accuracy of individual oscillating mammalian clock cells (26). It is our hope that these experimental preparations will soon be used to test the predictions of our model and that, because of the model's detail, these predictions will hold.
In line with previous work, we find that there can be significant differences in behavior between the stochastic and deterministic versions of a model. Similar to stochastic resonance, and in line with studies of much simpler models (27, 28), we find that our stochastic model can oscillate even when the corresponding deterministic model does not.
Amazingly, individual events on the promoter must occur on a very fast time scale, in our model, to have accurate 24h rhythms. This finding may explain some of the discrepancy between Barkai and Leibler (6), who report very irregular rhythms in a stochastic version of the Leloup and Goldbeter model of the intracellular clock in Drosophila (29), and Gonze et al. (7), who use infinitely fast reaction rates with the promoter in the same model and find very regular rhythms. This interpretation of the discrepancy is in agreement with that of Gonze et al. (30).
Stochastic simulations using the Gillespie method (as well as circadian clocks in vivo) differ from deterministic simulations in both their randomness and their coarseness. (Deterministic solutions of circadian clock models are smooth, whereas the corresponding stochastic solutions must jump between integers because they represent numbers of molecules.) In some limit, a stochastic solution may approach a deterministic solution, but without taking such a limit, it is, at best, a step function approximation to the deterministic solution. The more steps (i.e., the more molecules), the better the approximation will be. However, the quality of the approximation also depends on the time scale of the solution. There must be many steps (molecular events) within the time scale of the oscillation for good correspondence.
Likewise, the smaller the amplitude (relative to its mean) of a deterministic circadian oscillation, the more molecular events per unit time that are needed to accurately resolve the oscillation in a stochastic simulation. Thus, mechanisms for increasing the amplitude of circadian oscillations [e.g., changes in promoter cooperativity (30) or dimerization (data not shown)] can increase immunity to molecular noise. However, such noise reduction methods change the deterministic properties of the system and therefore must be weighed against any detrimental effects they might have on the clock's ability to entrain to environmental signal or function in a variety of environmental conditions.
The numerical experiments reported here suggest that certain aspects of the mammalian clock's structure might play a key role in noise immunity. Multiple copies of genes and rapid interactions with promoters are features that reduce the variability of the period of the clock. The robustness of circadian clocks appears to increase as more molecules are present or more frequent promoter interactions occur. In certain cases, notably the study of mutants, we find that the deterministic and stochastic models exhibit qualitatively different behaviors. This finding underscores the importance of using a stochastic model when possible, because the stochastic model corresponds more closely to reality.
Acknowledgments
We thank David Weaver, Chogoon Lee, Justin Blau, Willard Larkin, and J. Woodland Hastings for helpful discussions. This work was supported by Defense Advanced Research Planning Agency–Air Force Research Laboratory Grant F30602020554 and funding from The Sloan Foundation (to D.B.F.).
Footnotes

↵ † To whom correspondence should be addressed at: Department of Biology, New York University, 100 Washington Square East, New York, NY 10003. Email: forger{at}post.harvard.edu.

Author contributions: D.B.F. designed research, performed research, analyzed data, and wrote the paper; and C.S.P. participated in design of research, in analysis of data, and in writing of the paper.
 Copyright © 2005, The National Academy of Sciences
References
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Forger, D. B. & Peskin, C. S. (2003) Proc. Natl. Acad. Sci. USA 100 , 1480614811. pmid:14657377
 ↵

↵
Vitaterna, M. H., King, D. P., Chang, A. M., Kornhauser, J. M., Lowrey, P. L., Mcdonald, J. D., Dove, W. F., Pinto, L. H., Turek, F. W. & Takahashi, J. S. (1994) Science 264 , 719725. pmid:8171325
 ↵

↵
Femino, A. M., Fay, F. S., Fogarty, K. & Singer, R. H. (1998) Science 280 , 585590. pmid:9554849

↵
McClure, W. R. (1980) Proc. Natl. Acad. Sci. USA 77 , 56345638. pmid:6160577
 ↵
 ↵

↵
Raser, J. & O'Shea, E. (2004) Science 304 , 18111814. pmid:15166317

↵
Reppert, S. M. & Weaver, D. R. (2000) J. Biol. Rhythms 15 , 357364. pmid:11039914

Alberts, B., Bray, D., Lewis, J., Raff, M., Roberts, K. & Watson, J. (1994) Molecular Biology of the Cell (Garland, New York).

↵
Crosthwaite, S. K., Dunlap, J. C. & Loros, J. J. (1997) Science 276 , 763769. pmid:9115195
 ↵

↵
Jordan, D. & Smith, P. (1995) Nonlinear Ordinary Differential Equations (Clarendon, Oxford).

↵
Yamaguchi, S., Isejima, H., Matsuo, T., Okura, R., Yagita, K., Kobayashi, M. & Okamura, H. (2003) Science 302 , 14081412. pmid:14631044
 ↵

↵
Vilar, J. M., Kueh, H. Y., Barkai, N. & Leibler, S. (2002) Proc. Natl. Acad. Sci. USA 99 , 59885992. pmid:11972055

↵
Leloup, J. C. & Goldbeter, A. (1998) J. Biol. Rhythms 13 , 7087. pmid:9486845

↵
Gonze, D., Halloy, J. & Goldbeter, A. (2002) Proc. Natl. Acad. Sci. USA 99 , 673678. pmid:11792856
 ↵
Citation Manager Formats
More Articles of This Classification
Biological Sciences
Biophysics
Related Content
 No related articles found.
Cited by...
 Computational and experimental insights into the circadian effects of SIRT1
 Design Principles of PhosphorylationDependent Timekeeping in Eukaryotic Circadian Clocks
 Circadian clocks optimally adapt to sunlight for reliable synchronization
 Inverse Gillespie for inferring stochastic reaction mechanisms from intermittent samples
 Stochastic properties of the plant circadian clock
 A mechanism for robust circadian timekeeping via stoichiometric balance
 Signal processing in cellular clocks
 Synchronization and entrainment of coupled circadian oscillators
 A model of the cellautonomous mammalian circadian clock
 Temporal switching and celltocell variability in Ca2+ release activity in mammalian cells
 Modelling reaction kinetics inside cells
 Biological switches and clocks
 Synchrony and entrainment properties of robust circadian oscillators
 Rate constants rather than biochemical mechanism determine behaviour of genetic clocks
 Oscillation patterns in negative feedback loops
 Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks
 Reversible Protein Phosphorylation Regulates Circadian Rhythms
 An opposite role for tau in circadian rhythms revealed by mathematical modeling
 Balance between DBT/CKI{varepsilon} kinase and protein phosphatase activities regulate phosphorylation and stability of Drosophila CLOCK protein
 The relationship between FRQprotein stability and temperature compensation in the Neurospora circadian clock