QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

HOME  |   SEARCH  |   ARCHIVE  |   SUBSCRIBE  |   CONTACT  |   HELP

This Article Abstract Full Text (PDF) Submit an eLetter Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (65) Citing Articles via Google Scholar Google Scholar Articles by Smolen, P. Articles by Byrne, J. H. Search for Related Content PubMed PubMed Citation Articles by Smolen, P. Articles by Byrne, J. H.

 Previous Article  |  Next Article 

The Journal of Neuroscience, September 1, 2001, 21(17):6644-6656

Modeling Circadian Oscillations with Interlocking Positive and Negative Feedback Loops

Paul Smolen, Douglas A. Baxter, and John H. Byrne

Department of Neurobiology and Anatomy, W. M. Keck Center for the Neurobiology of Learning and Memory, The University of Texas-Houston Medical School, Houston, Texas

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Both positive and negative feedback loops of transcriptional regulation have been proposed to be important for the generation of circadian rhythms. To test the sufficiency of the proposed mechanisms, two differential equation-based models were constructed to describe the Neurospora crassa and Drosophila melanogaster circadian oscillators. In the model of the Neurospora oscillator, FRQ suppresses frq transcription by binding to a complex of the transcriptional activators WC-1 and WC-2, thus yielding negative feedback. FRQ also activates synthesis of WC-1, which in turn activates frq transcription, yielding positive feedback. In the model of the Drosophila oscillator, PER and TIM are represented by a "lumped" variable, "PER." PER suppresses its own transcription by binding to the transcriptional regulator dCLOCK, thus yielding negative feedback. PER also binds to dCLOCK to de-repress dclock, and dCLOCK in turn activates per transcription, yielding positive feedback. Both models displayed circadian oscillations that were robust to parameter variations and to noise and that entrained to simulated light/dark cycles. Circadian oscillations were only obtained if time delays were included to represent processes not modeled in detail (e.g., transcription and translation). In both models, oscillations were preserved when positive feedback was removed.

Key words: circadian; Drosophila; Neurospora; PER; CLOCK; FRQ; positive feedback; simulation

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Circadian rhythms reflect oscillating expression of genes, one or a few of which act as clock components, or core genes. The mechanisms by which core genes generate oscillations have been the subject of extensive experimental investigation. Negative feedback loops, involving indirect mechanisms of transcriptional repression, are now known for a few organismsnotably Drosophila melanogaster and Neurospora crassa. In Drosophila the transcriptional activators dCLOCK and CYCLE form a heterodimer that activates per and tim transcription (Darlington et al., 1998; Rutila et al., 1998; Lee et al., 1999; Bae et al., 2000). This heterodimer appears to be bound by PER and TIM to mask its DNA binding activity (Lee et al., 1999; Bae et al., 2000) and thereby repress per and tim transcription. In Neurospora, the white-collar proteins WC-1 and/or WC-2 activate frq transcription (Crosthwaite et al., 1997). The WC proteins may be bound by FRQ, and this may be the mechanism whereby FRQ represses frq transcription (Dunlap, 1999).

Recent results indicate that positive feedback also characterizes these circadian oscillators. In Drosophila, dCLOCK protein represses dclock transcription. Repression is mediated by dCLOCK-CYCLE heterodimers (Bae et al., 2000). PER and TIM activate dclock (Bae et al., 1998). A positive feedback loop can be envisioned as follows. If dclock expression is slightly activated, then the resulting activation of per and tim transcription by dCLOCK-CYCLE results in binding of dCLOCK-CYCLE by PER and TIM. Thus, repression of dclock is relieved, and the level of dCLOCK increases further. In Neurospora, increases in FRQ levels are followed by increases in WC-1 levels (Lee et al., 2000). Therefore, a positive feedback loop can be envisioned as follows. If wc-1 expression is initially activated, activation of frq transcription by WC-1 will increase the level of FRQ, which will in turn further increase the level of WC-1. A major question is the functional role of positive feedback in these systems. Intuitively, it appears that positive feedback might "cancel out" negative feedback. However, time delays between transcription and appearance of functional gene product (Garceau et al., 1997; Luo et al., 1998; Glossop et al., 1999) might minimize the conflict. Clearly, the interacting feedback loops with multiple time delays make it difficult to understand intuitively the mechanism of oscillation of even a simple configuration of core genes. Mathematical modeling is emerging as a powerful tool to supplement experimental work and provide insights into the operation of such gene networks (Leloup and Goldbeter, 1998, 2000; Leloup et al., 1999; Ruoff et al., 1999; Scheper et al., 1999; Gonze et al., 2000; Smolen et al., 2000a, 2000b).

The goal of the present investigation was to construct qualitative, differential equation-based models incorporating the above positive and negative feedback loops. Two models were constructed (Fig. 1) that represent, for Neurospora and Drosophila, positive and negative feedback of core genes on their own expression. Both models succeeded in simulating circadian oscillations that were robust to parameter variation and to stochastic noise in biochemical reaction terms. The oscillations would entrain to simulated light pulses or light/dark cycles. Positive feedback was not essential for simulation of circadian oscillations, although it may be important for driving "output" circadian genes regulated by core gene products. However, time delays between transcription and appearance of functional protein were required to be of considerable length (~7 hr).

View larger version (22K): [in this window] [in a new window]   Figure 1.   Models for the circadian oscillators in Neurospora and Drosophila. A, Neurospora. The FRQ gene product is multiply phosphorylated before degradation. WCC represents a complex of WC-1 and WC-2. All forms of FRQ are assumed equally competent at repressing frq transcription (dashed box) via binding to WCC. Degradation of WCC is included. The time delay 1 lies between changes in the level of WCC and resultant changes in the level of FRQ. 2 lies between changes in the level of FRQ and resultant changes in the level of WCC. B, Drosophila. "PER" represents a combination of PER and TIM levels. PER is multiply phosphorylated before degradation. All forms of PER are assumed equally competent at repressing the transcription of per (dashed box) via binding to dCLOCK. Degradation of dCLOCK is included. The time delay 1 lies between changes in the level of dCLOCK and resultant changes in the level of PER. 2 lies between a change in the level of free dCLOCK and the subsequent change in the level of total dCLOCK because of regulation of clock transcription.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

In the models of the Neurospora and Drosophila oscillators presented below, we do not consider separate concentration variables for the nuclear and cytoplasmic compartments. Rather, for simplicity, we model only the concentrations of macromolecular species averaged over a whole cell. Therefore, all concentrations are referenced to the total cell volume. Absolute in vivo concentrations are not well known for circadian proteins. We assumed standard nanomolar units.

A model of the Neurospora circadian oscillator. The model of the circadian oscillator in Neurospora that incorporates interactions between FRQ and the WC proteins is schematized in Figure 1A. The WC-1 and WC-2 proteins are found together in a complex in vivo (Talora et al., 1999; Denault et al., 2001). Recent evidence indicates that FRQ interacts with this complex and may thereby block its activation of frq expression (Denault et al., 2001). We use WCC (white-collar complex) as an abbreviation for the WC-1/WC-2 heterodimer, and we assume WCC activates frq expression. Our model contains a negative feedback loop in which FRQ represses frq transcription by binding and sequestering WCC, and a positive feedback loop in which activation of FRQ synthesis by WCC increases the level of FRQ, which in turn further increases the level of WCC. Experimentally, increases in FRQ are followed by increases in WC-1 (Lee et al., 2000) and thus plausibly in WCC.

The model incorporates a cascade of sequential phosphorylation events for FRQ protein, with only the most highly phosphorylated form being degraded (Fig. 1A). The number of FRQ phosphorylations that occur before degradation appears to be rather high (Garceau et al., 1997) but is not known definitely. Ten phosphorylations were assumed, as indicated by the notation "×10" in Figure 1A. F_0-10 denotes FRQ monomers with 0, 1, ..., 10 phosphorylations, respectively. Phosphorylations of FRQ were assumed to be sequential and irreversible. In the absence of detailed characterization of the enzymes involved, it was assumed for simplicity that the phosphorylations were all performed by a single kinase. Phosphorylation rates are given by Michaelis-Menten expressions in which all species of FRQ that are not fully phosphorylated can act to saturate the kinase. Michaelis-Menten terms represent degradation of fully phosphorylated FRQ and of WCC. All forms of FRQ are assumed able to form a complex with WCC. We use F-WCC to denote an FRQ-WCC complex.

The first differential equation represents synthesis and removal of unphosphorylated FRQ, F₀:

(1)

In this equation, the first term represents synthesis of F₀. The rate of synthesis of F₀ is proportional to a function, R_sF, of the concentration of free WCC (Eq. 2). The second term in Equation 1 represents removal of F₀ via phosphorylation, with a maximal velocity v_ph and a Michaelis constant K_ph. In the denominator, the quantity TOT_UNPHOS is a sum over the concentrations of all forms of FRQ that have <10 phosphorylations, whether bound to WCC or not. The use of this quantity reflects the assumption that a single kinase catalyzes all phosphorylations of FRQ. This kinase can therefore bind, and be saturated by, all species of FRQ that are not fully phosphorylated. The third term in Equation 1, R, represents removal of F₀ by association with WCC to form a complex F-WCC. The association process is assumed second-order, with rate constant k_f. Thus, R = k_f[F₀][WCC]. Dissociation of F-WCC is neglected for simplicity.

In the first term of Equation 1, the function R_sF takes the following form:

(2)

This form makes the rate of FRQ synthesis, v_sF * R_sF, a saturable and delayed function of the level of free WCC. A justification for Equation 2 is as follows. The rate of transcription of frq is assumed proportional to the probability that WCC is bound to an upstream regulatory sequence. This probability is assumed to be a saturable function of [WCC], yielding:

(3)

Because our model does not consider separate cytoplasmic and nuclear compartments, the concentrations in Equation 3, as in all other equations, are averaged over the total cell volume. To go to Equation 2, it is assumed that the rate of FRQ synthesis is proportional to the rate of frq transcription, but with a time delay ₁, taken as 7 hr. Much of ₁ is needed to account for the ~4 hr delay between the time courses of frq mRNA and FRQ protein during circadian oscillations (Garceau et al., 1997; Luo et al., 1998). ₁ would also include time needed for WC-1 and WC-2 to move to the nucleus and regulate frq.

The use of time delays such as ₁ to represent slow biochemical processes not modeled in detail is an effective way of encapsulating slow processes whose details are too complex or uncertain to model (Smolen et al., 1999). We used distributed, rather than discrete, time delays. With a distributed delay, the derivative of a variable depends on the averaged values of one or more variables over a specified range of previous time (MacDonald, 1989). For example, if the range of previous time is denoted by (t₁, t₂), a distributed delay for a variable X(t) whose rate of change depends on a function of itself is given by:

(4)

Distributed delays are general enough to encapsulate transcription, translation, or transport. For example, if movement of mRNA from a transcription site to translation sites relies on active transport with a range of transport times for individual mRNAs, a distributed delay is a proper modeling framework.

In equations defining delayed functions, angled brackets indicate that the quantity within is averaged over a time window centered at a mean delay  (₁ in Eq. 2). For convenience, the mean delay may be termed "the delay." The width of the time window is denoted . In Equation 4,  = (t₂  t₁).  is stated as a fraction of . In Equation 2, the standard value of  was taken to be 100% of ₁. Such a large value is consistent with experimental results. Garceau et al. (1997) and Luo et al. (1998) found that frq mRNA exhibits a rather narrow oscillation, with a width at half-maximum on the order of 6-8 hr, whereas the rising phase alone of oscillation in FRQ protein appears somewhat longer than this, lasting ~8-10 hr. These findings suggest a large variability in the times during which individual frq mRNAs are translated. This corresponds to a large .

The differential equation for the rate of change of the level of free WCC, the transcriptional activator complex, is as follows:

(5)

In this equation, the first two terms represent synthesis and degradation of WCC. In the denominator of the second term, WCC_tot denotes the total concentration of WCC, both free and bound to dCLOCK. The third term, R, represents removal of free WCC by association with any of the 11 species of FRQ (unphosphorylated, or with 1-10 phosphorylations). The sum runs over 11 association rates, one for each species of FRQ. Each association process is assumed second-order, with the same rate constant k_f in all cases. For example, R = k_f[F₁₀][WCC].

FRQ appears to activate the synthesis of one of the WC proteins, WC-1, post-transcriptionally (Lee et al., 2000). The expression of the second white-collar gene, wc-2, does not display circadian oscillations (Crosthwaite et al., 1997). Because WC-2 is in excess over WC-1 (Denault et al., 2001), activation of WC-1 synthesis should correspond to increasing the level of WCC. Therefore, to reflect the activation of WC-1 synthesis by FRQ, the rate of synthesis of WCC in Equation 5 is proportional to a delayed function R_sW of the total concentration of FRQ not complexed to WCC, [F_tot]. [F_tot] = [F₀] + [F₁] + ··· + [F₁₀]. Analogously to Equation 2, a saturable function was used:

(6)

The delay ₂ appears to be large (~7-8 hr; Lee et al., 2000). We used a value of 7 hr. The window width  was taken as 30% of the mean delay.

Phosphorylated species of FRQ can be formed and removed by phosphorylation reactions, except that removal of F₁₀ is by degradation. Also, these species can be removed by association with WCC to form FRQ-WC complexes. Each differential equation for these species therefore has three terms:

(7)

(8)

Within complexes of FRQ and WCC, FRQ can also be phosphorylated. The variables C are used to denote F-WCC complexes with 0, 1, ..., 10 FRQ phosphorylations. In Equations 1, 7, and 8, the variable TOT_UNPHOS is equal to [C] + [C] + ... + [C] + [F₀] + [F₁] + ... + [F₉]. There is no evidence to suggest the kinetics of FRQ phosphorylation depend on whether FRQ is complexed with WCC or not. Therefore, for simplicity, the same velocities and Michaelis constants are assumed to govern phosphorylation of FRQ in both cases. For the variables [C] through [C], differential equations analogous to Equation 7 can therefore be written. One additional term is added to account for degradation of these complexes if the WCC component is degraded. This term has the same form as that in Equation 5 for the degradation of free WCC:

(9)

(10)

The differential equation for the rate of change of [C] needs another term, to account for degradation of the complex upon degradation of fully phosphorylated FRQ. This term is analogous to that in Equation 8 for the degradation of fully phosphorylated FRQ:

(11)

For most simulations with the Neurospora model (Eqs. 1, 2, and 5-11), a standard set of parameter values was used, with exceptions noted in the text or figure legends. The standard set is:

Experimental data to estimate parameter values is lacking, except in the case of the delays ₁ and ₂, which were chosen as discussed above. Therefore, to obtain standard values for the other parameters, it was necessary to rely on trial-and-error variation. Values were found that allowed simulation of stable circadian oscillations robust to small parameter changes and simulation of entrainment to light pulses. Parameter values for the model of the Drosophila oscillator were obtained in the same manner.

A model of the Drosophila circadian oscillator. The PER-TIM heterodimer inhibits per and tim expression by binding to a heterodimer of the transcriptional activators dCLOCK and CYCLE and hindering binding of dCLOCK-CYCLE to the per and tim promoter regions (Lee et al., 1999; Bae et al., 2000). Because CYCLE is present in excess (Bae et al., 2000) it is likely that the concentration of dCLOCK-CYCLE is proportional to that of dCLOCK. Thus, the interaction of dCLOCK with CYCLE was not modeled explicitly. Rather, it was assumed that activation of per and tim is governed by the level of dCLOCK not complexed with PER-TIM. dclock transcription is activated by PER and TIM acting in concert. The activation is actually a de-repression. The PER-TIM heterodimer binds to and sequesters dCLOCK, which otherwise would repress its own gene (Glossop et al., 1999; Bae et al., 2000). Our model assumes this mechanism.

TIM alone does not appear to regulate transcription. Also, the time courses of PER and TIM proteins are similar in shape and largely overlap (Lee et al., 1998; Bae et al., 2000). Because of these considerations, the dynamics of PER and TIM may be represented by a single "lumped" variable, "PER," which can be thought of as an "average" of PER and TIM levels. Combining PER and TIM into a single variable is also economical. If PER and TIM were modeled separately the number of differential equations would increase greatly. Therefore, in our model, only the variable "PER" was used. Binding of dCLOCK with PER forms an inactive complex, denoted by P-C.

The model of the circadian oscillator in Drosophila, which includes binding of PER to dCLOCK, is schematized in Figure 1B. It embodies the negative feedback loop in which PER binds dCLOCK and therefore deactivates per transcription. In addition, it embodies the positive feedback loop that operates as follows. Activation of per transcription by dCLOCK results in sequestration of dCLOCK by binding to PER. Thus, repression of dclock by dCLOCK is relieved, and the level of dCLOCK increases further.

The model assumes sequential phosphorylation events for PER, with only the most highly phosphorylated form being degraded (Fig. 1B). The number of PER phosphorylations that occur is considerable (Edery et al., 1994; Zeng et al., 1996). Ten was assumed as was done above for FRQ. Michaelis-Menten terms represent degradation of PER and dCLOCK. All forms of PER, irrespective of phosphorylation, are assumed able to bind dCLOCK. The resulting model has 23 dependent variablesthe same as in our model of the Neurospora oscillator. The equations are almost identical to that model.

The first differential equation represents synthesis and removal of unphosphorylated PER, P₀. The rate of synthesis of P₀ is given as a function, R_sP, of the concentration of free dCLOCK:

(12)

In this equation, the first two terms represent synthesis and phosphorylation of P₀. In the first term, the function R_sP takes the form:

(13)

This form makes the rate of PER synthesis a saturable, delayed function of free dCLOCK. The mean delay ₁ was taken as 8 hr. Much of ₁ can be accounted for by the ~4 hr delay between the time courses of per mRNA and PER protein (Zerr et al., 1990; Vosshall et al., 1994; Glossop et al., 1999). An additional component of ₁ is evident as an experimental offset between the time course of the per transcription rate and the level of per mRNA. Surprisingly, this delay is on the order of 2 hr (So and Rosbash, 1997). Similar delays exist between the time courses of tim transcription, tim mRNA, and TIM protein (So and Rosbash, 1997). Movement of dCLOCK into the nucleus and binding to DNA would also add to ₁. In the denominator of the second term in Equation 12, TOT_UNPHOS is a sum of the concentrations of all forms of PER with <10 phosphorylations, whether bound to dCLOCK or not. The third term in Equation 12, R, represents removal of P₀ by association with dCLOCK to form a complex P-C. The association process is assumed second-order, thus R = k_f[P₀][dCLOCK]. Dissociation of P-C is neglected for simplicity.

The differential equation governing [dCLOCK] is:

(14)

In this equation the first two terms represent synthesis and degradation. The third term, R, represents removal of free dCLOCK by association with the species of PER with 0-10 phosphorylations. For example, R = k_f[P₁₀][dCLOCK]. In the first term, the function R_sC was chosen to reflect the repression of dCLOCK synthesis:

(15)

In Equation 15, ₂ represents the delay between changes in free dCLOCK and subsequent changes in the rate of appearance of new dCLOCK protein. Thus, ₂ encompasses dCLOCK transcription and translation. A value of 5 hr was used for ₂. A caveat is that ₂ has not been experimentally determined. Equation 15 embodies the key difference between the Drosophila and Neurospora models. In Equation 15, dCLOCK synthesis is repressed by dCLOCK protein. The analogous Equation 6 in the Neurospora model embodies activation of WCC synthesis by FRQ protein.

Phosphorylated species of PER can be formed and removed by phosphorylation reactions, except that for P₁₀, removal is by degradation. Also, these species can be removed by association with dCLOCK to form P-C complexes. Each differential equation for these species therefore has three terms:

(16)

(17)

Within complexes of PER and dCLOCK, PER can also be phosphorylated. The variables C are used for P-C complexes with 0, 1, ..., 10 PER phosphorylations. Differential equations analogous to Equations 9-11 in the Neurospora model govern the concentrations of these complexes:

(18)

(19)

(20)

For most simulations with the Drosophila model (Eqs. 12-20), a standard set of parameter values was used. The standard set is:



As in the model of the Neurospora oscillator, the standard value of the window width  was taken to be 100% of the mean delay for ₁ and 30% of the mean delay for ₂.

A list of all the parameters in the models of Figure 1, together with summaries of their biochemical significance, is presented in Table 1.

View this table: [in this window] [in a new window]   Table 1.   Parameters of the models

For numerical integration of the equations, the forward-Euler method was used with storage of averages for use in calculations of delayed quantities. Integration time steps were reduced until no significant difference was seen after further reduction. Final step sizes were 4-5 × 10⁵ hr. All models were programmed in FORTRAN 77 and simulated on a Compaq XP1000 workstation. Programs are available from the authors upon request.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Models incorporating feedback loops and time delays can represent the Drosophila and Neurospora circadian oscillators

Simulation of oscillations in Drosophila

The Drosophila model (Fig. 1B) readily simulated large-amplitude circadian oscillations in the level of total PER and total dCLOCK (Fig. 2, top panel). Here and subsequently, "total" protein concentration is the sum over all phosphorylation states and over free and bound states (e.g., for PER, the sum of the model variables [P₀] through [P₁₀] and [C] through [C]). The period of the simulated oscillations is 23.6 hr. Effects of light were not simulated in Figure 2, so these oscillations correspond to a rhythm in constant darkness. The oscillatory pattern in Figure 2 is stable over time and to modest changes in parameters as discussed further below.

View larger version (21K): [in this window] [in a new window]   Figure 2.   Simulation of circadian oscillations in Drosophila. Subsequent to an initial transient (data not shown), the model converged to oscillations that are independent of initial conditions. The standard set of parameter values (Materials and Methods) was used. Top panel, Time courses are displayed for the level of total PER (top, gray time course) and for the level of total dCLOCK (bottom, black time course). Middle panel, Expanded view displaying in more detail the dynamics of dCLOCK (top, gray time course) and of the level of the PER-dCLOCK complex (black time courseoften overlying the dCLOCK time course). Bottom panel, Time course of the level of free dCLOCK not complexed with PER.

The oscillations of PER concentration in Figure 2 are quite symmetric in appearance, with a gradual rise and fall. The time course of the concentration of total dCLOCK and that of the complex of PER and dCLOCK (the sum of [C] through [C]), are displayed in the middle panel of Figure 2. The oscillations of Figure 2 qualitatively agree with experiment in that the average level of dCLOCK is considerably less than that of PER (Bae et al., 2000).

The bottom panel of Fig. 2 illustrates the time course of free dCLOCK (not bound to PER), which is only at significant levels during the intervals when the level of total PER is below total dCLOCK. When [PER] rises at the beginning of an oscillation, the formation of PER-dCLOCK complex rapidly eliminates free dCLOCK. The loss of free dCLOCK acts to terminate per transcription. Phosphorylation and degradation of PER continue, and after several hours a decline in [PER] leads to dissociation of the PER-dCLOCK complex and the regeneration of free dCLOCK. Free dCLOCK can then activate per transcription and [PER] can again rise. It is evident from comparing the top and bottom panels of Fig. 2 that a considerable delay, on the order of the parameter ₁, separates changes in the level of free dCLOCK (such as the increase beginning at t = 6 hr) from the resulting changes in the rate of synthesis of PER protein (the increase beginning at t = 11 hr).

Oscillations of significant amplitude in the concentration of total dCLOCK are seen (Fig. 2, top panel). These oscillations reflect the indirect regulation of dCLOCK synthesis by PER, through sequestration of dCLOCK and de-repression of dclock transcription. In the middle panel of Figure 2, a circadian cycle of the level of PER-dCLOCK complex displays two humps. The first hump, peaking at t = 13 hr, is caused by the binding of newly synthesized PER to dCLOCK. [PER-dCLOCK] then declines because dCLOCK, both free and in complex, is being degraded. At t = 17 hr dCLOCK synthesis begins, and dCLOCK can bind with PER, so [PER-dCLOCK] rises again, to the second maximum at t  28 hr. It is not yet known experimentally whether such dynamics characterize the level of a PER-TIM-dCLOCK complex.

Experimentally, oscillations in PER during 24 hr light/dark cycles tend to peak at approximately zeitgeber time (ZT) 20 (Edery et al., 1994; Lee et al., 1998), where ZT0 is lights on and ZT12 is lights off. The dCLOCK time course is clearly later than that of PER throughout the day, peaking at approximately ZT0 (Bae et al., 2000). In the oscillations of Figure 2, a delay on the order of the parameter ₂, ~7 hr, separates the start of each increase in PER from the start of the resulting increase in dCLOCK. Thus, both in experiment and simulation, there is a significant lag between the PER and dCLOCK time courses.

One shortcoming of the simulated oscillations in Figure 2 is that as the level of PER rises at the beginning of an oscillation, PER-dCLOCK complex begins to form immediately, rather than after an interval of a few hours (Saez and Young, 1996). In the model, newly formed PER can always complex readily with dCLOCK. To more adequately represent the kinetics of PER-dCLOCK (or PER-TIM-dCLOCK) complex formation, a more detailed model would be necessary and would have to include kinetic expressions for PER transport into the nucleus.

The delay parameter ₁ is critical. Decreasing ₁ decreased the oscillation period. For example, a ₁ of 3 hr corresponds to a period of 12.8 hr with other parameters as in Figure 2. The period is more sensitive to ₁ than to any other parameter, and eliminating ₁ abolished oscillations. ₁ appeared essential, because oscillations could not be restored by varying other parameters. In contrast, if ₂ was eliminated, oscillations were preserved. The period is relatively insensitive to ₂. Decreasing ₂ acts to decrease the lag between increases or decreases in [PER] and the resulting increases or decreases in [dCLOCK].

Simulation of oscillations in Neurospora

The Neurospora model (Fig. 1A) readily simulated large-amplitude circadian oscillations in the levels of total FRQ and total WCC (Fig. 3, top panel). The period of these oscillations is 24.0 hr. Effects of light were not simulated in Figure 3, so these oscillations correspond to a free-running rhythm in constant darkness. The oscillatory pattern in Figure 3 is stable over time, and the oscillations are stable to modest changes in parameters (see Fig. 6). Decreasing ₁ decreased the oscillation period. Elimination of ₁ abolished oscillations. Again, ₁ appeared essential, because oscillations could not be restored by varying other parameters. For the standard parameter values, elimination of ₂ also abolished oscillations, but these could be restored by increasing v_sF from 8.0 to 16.0 nM/hr.

View larger version (21K): [in this window] [in a new window]   Figure 3.   Simulation of circadian oscillations in Neurospora. Subsequent to an initial transient, the model converged to oscillations that are independent of initial conditions. The standard set of parameter values (Materials and Methods) was used. Top panel, Time courses are displayed for the level of total FRQ (top, gray time course) and for the level of total WCC (black time course). Middle panel, Expanded view displaying in more detail the dynamics of WCC (top, gray time course) and of the level of the FRQ-WCC complex (black time courseoften overlying the WCC time course). Bottom panel, Time course of the level of free WCC not complexed with FRQ.

The oscillations in Figure 3 display a gradual rise and fall in the concentration of total FRQ. The middle panel of Figure 3 displays in more detail the time courses of total WCC and of the complex of FRQ and WCC (the sum of [C] through [C]). As [FRQ] rises, WCC is rapidly bound so that soon essentially all WCC is complexed with FRQ. Over an oscillation, the average concentration of WCC is ~40% that of FRQ. Experiments have not yet established the relative levels of FRQ and WCC. In Fig. 3 (top), the oscillations of FRQ and WCC are approximately antiphase. A large value of ₂ is required for this relationship. Experimental oscillations of FRQ and WC-1 are also approximately antiphase (Lee et al., 2000).

Figure 3, bottom panel, illustrates the time course of the concentration of free WCC. Free WCC is only at significant levels during the intervals when [FRQ] is below [WCC]. When [FRQ] rises at the beginning of an oscillation, the formation of FRQ-WCC complex rapidly eliminates free WCC. The loss of free WCC acts to terminate frq transcription. Phosphorylation and degradation of FRQ continue, and after ~15 hr a decline in FRQ leads to dissociation of the FRQ-WCC complex and the regeneration of free WCC. Free WCC can then activate frq transcription, and [FRQ] can again rise. It is evident from comparing the top and bottom panels of Figure 3 that a considerable delay, on the order of the parameter ₁, separates changes in the level of free WCC (such as the increase beginning at t = 28 hr) from the resulting changes in the rate of synthesis of FRQ protein (the increase beginning at t = 33 hr).

Comparison of simulated and experimental oscillations

Figure 4A compares experimental time courses of Neurospora FRQ and WC-1 levels (Lee et al., 2000) with time courses of FRQ and WCC simulated by our model (Fig. 3). Absolute in vivo concentrations of these proteins are not known. Thus, the simulated time courses were scaled vertically to match the magnitude of the experimental time courses, and specific units were not used for amounts of proteins in Figure 4. The experimental WC-1 time course is used because no experimental time course for WCC, the complex between WC-1 and WC-2, is available. The overall shape of the simulated FRQ time course, including the relative timing of the peak versus the minimum, appears qualitatively consistent with experiment when the random variability between experimental oscillations is considered (Garceau et al., 1997; Lee et al., 2000). The shape of the simulated WCC time course is quite similar to that of the WC-1 time course. Experimentally, because WC-2 is in excess (Denault et al., 2001), WC-1 levels should determine the levels of WCC, so that one might expect the shapes of WCC and WC-1 time courses to be similar. Our model is qualitative and does not include separate cytoplasmic and nuclear compartments. The dynamics of FRQ and of the WC proteins are expected to differ between compartments, for example, WC-2 and WCC are almost entirely nuclear proteins (Denault et al., 2001). A more detailed model would be required to capture these effects. For this reason, and because of the considerable random variation between experimental FRQ oscillations, more detailed fitting of simulated time courses to experimental data via extensive parameter variation was not performed.

View larger version (20K): [in this window] [in a new window]   Figure 4.   Comparison of simulated and experimental time courses of circadian gene product levels. A, Neurospora. The simulated time course of FRQ (red) is compared with experimental time course (gray with squares). The simulated time course of WCC (green) is compared with the experimental WC-1 time course (black with squares). Data is from Lee et al. (2000), their Figure 1B. B, Drosophila. The simulated time course of PER (red) is compared with experimental time courses of PER and TIM (black and gray). Data is from Lee et al. (1998), their Figure 2B.

Figure 4B compares experimental time courses of Drosophila PER and TIM levels (Lee et al., 1998) with the simulated time course of PER. Both experimental time courses are shown because, as discussed in Materials and Methods, "PER" in our model represents an average of PER and TIM. Experimental time courses are during a 12 hr light/dark cycle. Therefore, as discussed in the following section, the simulated time course is also during a 12 hr light/dark cycle (parameters as in Fig. 5C). The locations of the peak and minimum, and the shape of the rising phase, are similar in the experimental and simulated time courses. The decline of the simulated time course is somewhat steeper than experiment, however.

View larger version (22K): [in this window] [in a new window]   Figure 5.   Entrainment of circadian oscillations simulated by the Neurospora model (A, B) and the Drosophila model (C). A, Entrainment of the oscillations of Figure 3 by simulated light pulses. The interstimulus interval was 20 hr. Short square bars mark the effect of each light pulse, which was assumed to induce an increase in FRQ synthesis. Each 90 min increase initiates an upstroke in total FRQ. The time course of total WCC is also shown. B, Limits of entrainment for the interstimulus interval of the simulated light pulses of A. The free-running oscillation period (24.0 hr) is marked. C, Entrainment of the oscillations of the Drosophila model (Fig. 2) and the Neurospora model (Fig. 3) by a 12 hr light/dark cycle. The light and dark phases are indicated by the overbar. Each light phase was assumed to induce a first-order degradation of phosphorylated PER. The degradation rate constant was 0.9 hr1. Each light phase also decreased the degradation velocity for dCLOCK, vdC, to 1.5 nM/hr. For Neurospora, each light phase induced a constant synthesis of unphosphorylated FRQ (3.0 nM/hr).

The Neurospora and Drosophila models can simulate entrainment by light

Models of circadian rhythms must be able to simulate responses of the rhythm to light pulses or light/dark cycles (Leloup and Goldbeter, 1998, 2000; Scheper et al., 1999; Gonze et al., 2000). In Neurospora, an increase in frq transcription is observed after a light pulse. Increased frq mRNA levels persist for somewhat longer than 60 min (Crosthwaite et al., 1995). Thus, each light pulse was modeled as an increase in the first term on the right-hand side of Equation 1, which denotes the rate of synthesis of unphosphorylated FRQ (the variable F₀). A delay of a few hours, not explicitly included in simulations, would separate light pulses from increases in FRQ synthesis, because of processing and transport of frq mRNA (Garceau et al., 1997).

In one set of simulations, light pulses were modeled by increasing the velocity of synthesis (v_sF) of unphosphorylated FRQ, F₀, to a constant rate of 10.0 nM/hr for 1.5 hr. Parameter values were otherwise as in Figure 3. Figure 5A demonstrates entrainment of oscillations by these simulated light pulses, with an interstimulus interval (ISI) of 22 hr. The brief increases in FRQ synthesis are each seen to initiate upstrokes of the level of total FRQ from a very low value. ISIs of longer than the free-running period (24 hr) can also entrain the oscillations of Figure 3. With the stimulus parameters of Figure 5A, the limits of entrainment for the ISI are 20-27 hr, as illustrated in Figure 5B.

In Drosophila, light enhances the degradation of phosphorylated TIM (Myers et al., 1996; Zeng et al., 1996). Because the Drosophila model does not have a separate variable for the level of TIM, the degradation of phosphorylated PER was enhanced to simulate the effect of light. The effect of such enhancement is to free dCLOCK from the PER-dCLOCK complex, allowing dCLOCK to regulate transcription. In vivo, the degradation of TIM may analogously release free dCLOCK. The effect of light on subsequent transcription of core circadian genes might therefore be similar in simulations with our model and in vivo. To simulate light exposure, first-order degradation rate constants (k_deg[P_i]) were assumed for the phosphorylated forms of PER (the variables [P₁], ..., [P₁₀] and [C], ..., [C]) during the light phase of light/dark cycles or during light pulses. When PER complexed with dCLOCK was degraded, it was assumed that free dCLOCK is released.

Figure 5C demonstrates entrainment of the simulated Drosophila and Neurospora circadian oscillations to a light/dark cycle. The cycle used was 12 hr. During the 12 hr of "light," a first-order degradation rate constant of 0.9 hr¹ was assumed for all phosphorylated forms of PER. Because the enhanced PER degradation acted to advance the phase of the succeeding peak in the concentration of free dCLOCK, and then the phase of the next peak in [PER], a compensating parameter change was necessary to maintain an oscillation period of 24 hr. It was therefore assumed that dCLOCK degradation was reduced during the light phase (v_dC was decreased to 1.5 nM/hr). This change enhanced and prolonged the succeeding peak of free dCLOCK and the next peak in [PER]. It appears plausible that during the light phase, enhanced PER degradation could use elements of the pathway mediating dCLOCK degradation, competitively reducing v_dC. For Neurospora, each light phase was modeled as addition of a constant synthesis rate of FRQ to Equation 1 (3.0 nM/hr). In Figure 5C, discontinuities in the slope of the time course of total PER are evident that coincide with the onset of "light" and the imposition of the degradation process. The switch from dark to light occurs on the falling portion of the PER time course, as it does experimentally (Lee et al., 1998; Bae et al., 2000). With the stimulus parameters of Figure 5C, the shape of the time course of total PER during a simulated oscillation is significantly different from the time course in constant darkness (Fig. 2). With the light/dark cycle, the falling portion of the time course is considerably steeper than the rising portion. Experimentally, during entrainment to light/dark cycles, the falling portions of the PER and TIM time courses are indeed steeper than the rising phases (Marrus et al., 1996; Lee et al., 1998). Comparing the FRQ and PER time courses in Figure 5C, it is seen that the light phase coincides with the rising portion of the FRQ time course but with the falling portion of the PER time course. This agrees with observation (Dunlap, 1999). Entrainment of the oscillations of the Drosophila model by light pulses could also be simulated (data not shown). Entrainment to pulses with an ISI considerably shorter than 24 hr was readily obtained. However, entrainment to pulses with an ISI longer than ~25 hr was not obtained.

Simulated circadian oscillations are robust to parameter variation

Biochemical parameters are expected to vary somewhat from cell to cell and from one member of a species to another. Nevertheless, individual Neurospora or Drosophila are generally observed, in constant darkness, to sustain circadian rhythms with similar period. Because circadian rhythmicity is well preserved from one individual to the next, it is important that in a model of circadian rhythmicity, small parameter variations such as might be expected among individuals should not cause large changes in the period or amplitude of simulated circadian oscillations.

For simulations to test robustness of oscillations to small parameter variations, the standard parameter value sets (Materials and Methods section) were used as the starting point. With the Neurospora and Drosophila models, each individual parameter was increased and decreased by 15% of its standard value. For each model, there are 13 parameters, including the delays ₁ and ₂. Therefore, 27 simulations were performed including the control with standard parameter values. Figure 6 plots the period and amplitude of these simulations for the Neurospora model. The amplitude was measured as the peak-to-minimum difference in the oscillation of the concentration of total FRQ. In four additional simulations, also plotted, the widths of the distributed delays _1-2 were also increased and decreased by 25% from their standard values.

View larger version (11K): [in this window]