INTRODUCTION

TOP ABSTRACT INTRODUCTION 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 transcriptional repression, have been well characterized for a few organisms, notably Drosophila melanogaster. In Drosophila, the transcriptional activators dCLOCK and CYCLE form a heterodimer that activates per and tim transcription (Bae et al., 2000; Darlington et al., 1998). PER appears to bind the dCLOCK-CYCLE heterodimer so as to mask its DNA-binding activity (Bae et al., 2000; Lee et al., 1999) and thereby repress per and tim transcription. Positive feedback is also an important element of the Drosophila oscillator. The level of the core gene product dCLOCK oscillates (Lee et al., 1998) and represses dclock transcription (Glossop et al., 1999). PER appears to activate dclock by binding dCLOCK and blocking repression (Glossop et al., 1999; Bae et al., 1998). The positive and negative feedback loops are interlocked, because transcriptional regulation by dCLOCK is common to both loops.

Mathematical modeling has emerged as an important tool for gaining understanding of the dynamics of gene networks incorporating feedback loops and delays (Smolen et al., 2000; Keller, 1995). Detailed models of the Drosophila oscillator have been published (Smolen et al., 2001; Leloup and Goldbeter, 1998). These models consider multiply phosphorylated forms of PER, and the most recent model (Smolen et al., 2001) also represents complexes of PER with dCLOCK. These models have demonstrated that the negative and positive feedback loops, as currently understood, could sustain circadian oscillations of gene expression. Several simpler models have also been proposed to describe circadian rhythm generation (Lema et al., 2000; Ruoff et al., 1999; Scheper et al., 1999). However, these models do not represent both positive and negative feedback loops.

The goal of the present study was to construct a simplified model that represents the dynamics of the positive and negative feedback loops of the Drosophila oscillator, but that is implemented with as few differential equations as possible. Our earlier, detailed model (Smolen et al., 2001) was reduced to obtain a minimal representation of the feedback loops and their interactions. It is well established that such reduced models can greatly aid intuitive understanding of the dynamics of biophysical, biochemical, or genetic systems (Ermentrout, 2001; Smolen et al., 2000). The reduced model consists of two differential equations, each with a time delay. These delay differential equations describe the evolution of PER and dCLOCK concentrations. The delay differential equation for the evolution of [PER] is similar to that in the model of Lema et al. (2000). However, in the model of Lema et al. (2000), PER represses per transcription directly rather than by binding dCLOCK.

The reduced model succeeded in simulating circadian oscillations of [PER] and [dCLOCK]. The oscillation amplitude and period were robust to parameter variation. The oscillations entrained to simulated light pulses or light-dark cycles. For light pulses, a phase-response curve was constructed. The shape shared qualitative similarities with experimental curves.

It has recently been emphasized that simulated circadian oscillations should also be robust to random fluctuations in the molecule numbers of key gene products (Barkai and Leibler, 2000). Previous reduced models of circadian rhythmicity have not incorporated such stochastic fluctuations. Therefore, a stochastic version of the Drosophila model was constructed. This simulated robust circadian oscillations, with little variability in period. The average numbers of PER and dCLOCK molecules were both less than 80.

The role of the positive feedback loop in Drosophila is an issue of interest (Hastings et al., 2000). To examine possible roles of positive feedback, it was removed from the model by fixing the total amount of the transcriptional activator dCLOCK. Robust circadian oscillations in PER were preserved. We are able to provide an intuitive understanding of this result, because when the activators are fixed, both models reduce to a single delay differential equation with negative feedback. The oscillations simulated by this single equation are robust to variations in parameter values and to stochastic fluctuations in molecule numbers. These results suggest that the primary role of positive feedback may be to drive oscillations in the level of dCLOCK and in the expression of clock-controlled genes regulated by dCLOCK.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Model development and numerical methods

For simplicity, separate nuclear and cytoplasmic compartments are not considered below. Rather, concentrations are referenced to the total cell volume. Absolute concentrations are not well known for circadian proteins. We assume that the maximum concentration of PER during an oscillation is ~0.5 nM, which corresponds to ~250 PER molecules in a Drosophila lateral neuron with a radius of ~6 µm (Ewer et al., 1992).

Model of the Drosophila circadian oscillator

This model is schematized in Fig. 1 A. dCLOCK activates per transcription and thus PER synthesis. PER represses per transcription (and thus PER synthesis) by binding dCLOCK. PER also activates dCLOCK synthesis by binding dCLOCK and relieving dCLOCK's repression of dclock transcription. Thus, the model contains both a negative feedback loop, in which PER binds dCLOCK and thereby de-activates per transcription, and a positive feedback loop, in which activation of per transcription by dCLOCK results in binding of dCLOCK by PER and de-repression of dclock.

View larger version (24K): [in this window] [in a new window]   FIGURE 1   Simulation of circadian oscillations by a reduced model. (A) Schematic of model. The dCLOCK protein activates the synthesis of PER. PER represses its own synthesis indirectly, by binding and inactivating dCLOCK. dCLOCK also represses its own synthesis. A time delay 1 is included between changes in dCLOCK concentration and in PER synthesis, and a delay 2 is included between changes in dCLOCK concentration and in dCLOCK synthesis. (B) Simulated oscillations. Blue, dark green, and light green traces are, respectively, for [PER], [dCLOCK], and [dCLOCK]free. The standard set of parameter values, given in Methods, is used.

The model of Fig. 1 A is similar to our previous, more detailed model (Smolen et al., 2001) in that it neglects the dynamics of two other gene products known to be involved in the generation of circadian oscillations in Drosophila-TIM and CYCLE. The model of Fig. 1 A is reduced from the previous model in two ways. Multiple phosphorylations of PER are neglected, and complexes of PER and dCLOCK are not modeled explicitly. These simplifications, discussed further below, allow for a model with only two dependent variables, [PER] and [dCLOCK]. These variables refer to total concentrations of these proteins, whether free or bound together.

The differential equations for [PER] and [dCLOCK] have two terms, one for synthesis and one for degradation. Because regulation of degradation was not included, simple first-order degradation rate constants were assumed. The differential equation for [PER] is:

(1)

PER synthesis is assumed to be activated by free dCLOCK. With PER-dCLOCK complexes not modeled explicitly, the reduced representation of free dCLOCK is the difference [dCLOCK]  [PER]. This representation assumes that PER immediately binds any available dCLOCK. This assumption is not likely to be quantitatively correct, but it appears reasonable for a simplified and qualitative model. In the synthesis term of Eq. 1, the function R_sP is assumed to depend hyperbolically on free dCLOCK as follows:

(2)

with [dCLOCK_free] = ([dCLOCK]  [PER]) or zero, whichever is greater. Equations 1 and 2 are similar to a recent model of circadian oscillations based on a single delay differential equation (Lema et al., 2000). The model of Lema et al. (2000) consists of Eq. 1 and an expression similar to Eq. 2. In that expression, PER represses its own synthesis directly.

The parameter ₁ in Eq. 2 denotes the time delay between per transcription and the synthesis of new PER protein. This discrete time delay is implemented as follows. The fraction within the angle brackets represents activation of per transcription by dCLOCK. At each time step of a simulation, the value of this fraction is stored. The stored value is used ₁ h later to compute the rate of PER synthesis. ₁ includes the long (~5 h) delay between the time courses of per mRNA and PER protein (Glossop et al., 1999; Vosshall et al., 1994). ₁ also includes an ~2 h offset between the time course of per transcription and the time course of per mRNA (So and Rosbash, 1997).

The differential equation for [dCLOCK] is:

(3)

In the synthesis term of Eq. 3, the function R_sC represents repression of dclock transcription (and thus dCLOCK synthesis) by free dCLOCK:

(4)

₂ in Eq. 4 denotes the time delay between dclock transcription and the synthesis of new dCLOCK protein. ₂ has not been experimentally determined.

As mentioned above, the Drosophila model neglects a slow process that contributes to generating the circadian oscillation period. This process consists of multiple phosphorylations of PER protein before degradation (Edery et al., 1994). Our previous, detailed model (Smolen et al., 2001) included these phosphorylations, but that model had in excess of 20 dependent variables. One might think slow PER phosphorylation could be incorporated in Eq. 1 by adding another time delay. The degradation rate of PER might be made a delayed function of [PER]. However, in this case [PER] is no longer constrained to remain nonnegative over time. Simulations with this model variant were carried out, and [PER] was often observed to become negative. Thus this approach was rejected. Instead, the effective delay contributed by PER phosphorylations was incorporated into the delays ₁ and ₂ in Eqs. 2 and 4. Thus, these delays were assumed to be longer than in our previous model with PER phosphorylation (Smolen et al., 2001). In that model, values of ~7 h were used, whereas in the present model, ₁ and ₂ were set to 10 h. The present model also does not incorporate the posttranscriptional regulation of [dCLOCK], which has recently been demonstrated (Kim et al., 2002). Although the stability of dCLOCK is evidently regulated, not enough appears known about the mechanism to justify modeling it.

For most simulations with the Drosophila model (Eqs. 1-4), a standard set of parameter values was used, as follows: ₁ = 10 h, ₂ = 10 h, v_sP = 0.5 nM h¹, v_sC = 0.25 nM h¹, k_dP = 0.5 h¹, k_dC = 0.5 h¹, K₁ = 0.3 nM, and K₂ = 0.1 nM.

Experimental data to constrain values of kinetic parameters are generally lacking. As discussed above, some data exist to constrain ₁. 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.

Simulation of stochastic fluctuations in molecule numbers

Simulating fluctuations in molecule numbers requires, at each time step, a probabilistic determination of whether each type of chemical reaction takes place or not. For the Drosophila model, the procedure was as follows. The standard parameter value set was used as the starting point. Then, enzyme reaction velocities and Michaelis constants were rescaled so that the units of concentration variables were no longer nanomolar, but rather absolute numbers of molecules. To accomplish this, the parameters v_sP, v_sC, K₁, and K₂ were all multiplied by a common factor. Because the numbers of each type of molecule present per nucleus are not known accurately in Drosophila, the value of the common factor is arbitrary. As the factor is increased, the average molecule numbers increase. A factor of 250 was determined by trial and error to yield molecule numbers that simulated oscillations not overly degraded by large fluctuations.

After rescaling, at each time step of a simulation, the reaction probabilities were computed by multiplying the time step with the terms in Eqs. 1 and 3 that give the rates of the specific reactions. Two reaction terms create PER and dCLOCK. From Eqs. 1 and 3, these terms are v_sP × R_sP and v_sC × R_sC. The other two terms remove PER and dCLOCK; these terms are k_dP[PER] and k_dC[dCLOCK]. Multiplying these terms by the time step t gives the reaction probabilities per time step. The time step was fixed at a small value (5 × 10⁶ h) chosen so that the probability of each biochemical reaction was never larger than 2%. Therefore, in any time step the chance that the copy number of any given molecule will change by 1 is never larger than 4% (2% multiplied by 2 reactions because each protein is subject to two independent processes of synthesis and degradation). By using these small time steps, the probability of more than one reaction occurring in a time step can be considered negligible. Finally, at each time step a separate random number was generated for each reaction. Each random number was picked from a uniform distribution over (0,1). If the random number for a reaction was less than the probability of that reaction, then the reaction was assumed to occur. If the reaction synthesized or degraded a particular molecule, the copy number was changed by 1 or 1. For small time steps, this fixed time step algorithm is an explicit simulation of the master equation. To verify accuracy, simulations were repeated with the time step halved.

For some stochastic simulations, another algorithm was used. The Gillespie algorithm (Gillespie, 1977) takes time steps of varying length. It uses the following result from statistical physics. Suppose a given biochemical reaction occurs at a time arbitrarily taken as 0. Then, the probability P(t) that the next reaction of that type will occur within a small time interval (t, t + t) beginning at any specific time t > 0 is:

(5)

In Eq. 5, T_avg is defined as the reciprocal of the average, deterministic reaction rate. For example, after conversion to units of molecule numbers, T_avg for PER synthesis is given as (v_SP × R_SP)¹, with R_sP as in Eq. 2. At the beginning of each simulation time step, the set of average reaction rates is calculated, along with the sum of these rates. Denote the set of reaction rates by a_i, with i = 1, ... , m. For the Drosophila model, m = 4 because there is one reaction for synthesis of each protein and one for its degradation. Denote the sum of the reaction rates by a_TOT. Two random numbers (r₁, r₂) are picked from a uniform distribution on (0,1). The index i of the reaction that occurs during the time step is the value of i that satisfies the inequality:
The length of the time step is given by using the second random number:
To apply this algorithm to models with time delays, the appropriate delayed quantities were used to calculate the average reaction rates. For example, the average rate of PER synthesis was calculated by Eq. 2. To convert units to molecule numbers, parameters were rescaled by the same factor (250) as used for the fixed time step algorithm. For the simulations presented in this paper, both the Gillespie algorithm and the fixed time step algorithm ran at similar speeds (less than 30% difference in computer time).

Numerical methods

For integration of delay differential equations, the forward-Euler method was used, with storage of delayed quantities for later calculations. Integration time steps were reduced until no significant difference was seen upon further reduction. Final step sizes were ~2 × 10⁵ h. In deterministic simulations without fluctuations, and in stochastic simulations with the fixed time step algorithm, delayed quantities were updated by recall from storage every 0.05 h. In stochastic simulations with the Gillespie algorithm, time steps are of variable size, so delayed quantities were stored along with the time at which they were computed. Every 0.05 h, the entries in the storage arrays that were computed closest to ₁ and ₂ h previously were recalled. 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 METHODS RESULTS DISCUSSION REFERENCES

A reduced model with feedback loops and time delays can simulate Drosophila circadian oscillations

The Drosophila model (Fig. 1 A; Eqs. 1-4) readily simulated large-amplitude circadian oscillations in the levels of PER and dCLOCK (Fig. 1 B). Effects of light were not simulated, so these oscillations correspond to a free-running rhythm in constant darkness, with a period of 23.2 h. The oscillatory pattern was stable over time, and the oscillations were stable to modest changes in parameters (Fig. 2 A, discussed below). The oscillations in dCLOCK reflect the indirect activation of dCLOCK synthesis by PER, through sequestration of dCLOCK and de-repression of dclock transcription.

View larger version (14K): [in this window] [in a new window]   FIGURE 2   Effect of parameter variations on oscillations. (A) Robustness of oscillations. A scatter plot displays the periods and amplitudes of circadian oscillations in [PER]. To generate these oscillations, each value in the set of standard parameter values was increased or decreased by 20%. There are eight parameters and, therefore, 17 data points including the control with all values standard. The control point is at the intersection of the black lines. (B) Simulated and experimental photic PRCs. In constructing the model PRC (), each light pulse was simulated as an increase, for 5 h, in the first-order rate constant for PER degradation. The increase was 5.0 h1. Parameter values not affected by light exposure are as in Fig. 1 B. An experimental Drosophila PRC (Fig. 5 of Konopka et al., 1991, is also illustrated (). The means of the experimental phase shifts are displayed.

Fig. 1 B also illustrates the time course of free dCLOCK (not bound to PER), which exists only when the level of PER is below that of dCLOCK. The mechanism of oscillation is as follows. When [PER] rises at the beginning of an oscillation (at t  12 h), free dCLOCK is rapidly eliminated. This decreases the right-hand side of Eq. 2 and, after the delay ₁, terminates PER synthesis. The loss of free dCLOCK also increases the right-hand side of Eq. 4. As a result, after the delay ₂, dCLOCK synthesis is initiated (at t  22 h). Degradation of PER along with new dCLOCK synthesis rapidly regenerates free dCLOCK (at t  24 h). After another delay ₁, the free dCLOCK activates PER synthesis, beginning the next oscillation (at t ~34 h).

The delay ₁ is critical for oscillations. Decreasing ₁ decreased the oscillation period. For example, a ₁ of 5 h corresponds to a period of 17.7 h with other parameters as in Fig. 1 B. Eliminating ₁ always abolished oscillations, which could not be restored by varying other parameters. In contrast, when ₂ was eliminated oscillations were preserved, although v_sP had to be increased slightly (to 0.8 nM h¹). The oscillation period decreased to 21.5 h, and the lags between upstrokes of [PER] and succeeding upstrokes of [dCLOCK] were eliminated.

Simulation of mutations affecting PER phosphorylation

In Drosophila, the doubletime-S mutation accelerates PER phosphorylation and subsequent degradation, shortening the period to ~18 h (Price et al., 1998). In the Drosophila model the time required for phosphorylations of PER was incorporated in the delays ₁ and ₂, as discussed above. Therefore, to simulate mutations that accelerate or retard phosphorylation, ₁ and ₂ were, respectively, diminished or enhanced. Decreasing ₁ and ₂ to 7 h decreases the period to 17.1 h, similar to doubletime-S. Another mutation, perS, also shortens the oscillation period to ~19 h, and this shortened period is associated with an increased instability of PER, which may be due to accelerated phosphorylation of PER (Marrus et al., 1996). Thus, the simulation of the doubletime-S mutation may also represent the perS mutation.

The reduced model is robust to parameter variation and can simulate light responses

Sensitivity of oscillations to parameter variation.

Simulation of light responses

1

1

Simulated oscillations are robust to fluctuations in molecule numbers

Recently, the importance of testing models of circadian rhythmicity for robustness to stochastic noise has been emphasized (Barkai and Leibler, 2000; Smolen et al., 2000). This noise is due to stochastic fluctuations in the numbers of molecules, because of the random timing of individual biochemical reaction events. Barkai and Leibler (2000) point out that current models of circadian rhythmicity usually do not incorporate stochastic fluctuations. For any given model, inclusion of fluctuations might be found to produce unacceptably large random variation in oscillation period or amplitude. To test the model of Fig. 1 A, stochastic simulations were performed to determine the minimal numbers of protein molecules necessary to sustain oscillations without large random variations in period or amplitude. For smaller average molecule numbers, random fluctuations will always be relatively larger, and for sufficiently small numbers, such fluctuations would destroy periodic oscillations.

For the model of Fig. 1 A, stochastic fluctuations were incorporated by an algorithm that uses fixed time steps to simulate the master equation governing transitions between all possible sets of molecule numbers (see Methods for details). Fig. 3 A illustrates that despite fluctuations, a robust oscillation pattern was preserved. Parameters were as in Fig. 1 B except that concentration units were converted to numbers of molecules (see Methods). When averaged over 50 oscillations, the mean period was circadian and the standard deviation was modest (23.0 ± 1.1 h). These oscillations were sustained with mean molecule numbers of less than 80. Over 50 oscillations, the mean numbers were 70 for PER and 76 for dCLOCK.

View larger version (29K): [in this window] [in a new window]   FIGURE 3   Simulation of oscillations when fluctuations in molecule numbers are incorporated in the model of Fig. 1 A. (A) The fixed time step algorithm (Methods) was used to simulate fluctuations. Blue and green traces are, respectively, for PER and dCLOCK. (B) Comparison of oscillations simulated by the fixed time step algorithm (blue PER time course) and the Gillespie algorithm (brown PER time course). (C) Entrainment of oscillations by simulated light pulses. Parameters were as in A. Light pulses were modeled as an increase from 0.5 h1 to 5.0 h1 in the first-order rate constant for PER degradation. The duration of each increase was 3 h.

Relatively few simulations of stochastic models with time delays have been reported. Also, analytic understanding of the dynamics of such models is limited. For models based on a single linear stochastic differential equation, methods exist for determining steady states and analyzing their stability to perturbations (Mackey and Nechaeva, 1995; Ohira and Yamane, 2000). However, few analytical results exist concerning the dynamics of stochastic differential equations with delay, such as those in our models. As a consequence, it is an open question which algorithms are best for simulating such models. Because of this situation, we compared simulations with the fixed time step algorithm to simulations using the Gillespie algorithm (Gillespie, 1977). The Gillespie algorithm uses variable time steps to simulate stochastic fluctuations (see Methods). Fig. 3 B compares PER oscillations simulated by the fixed time step algorithm and the Gillespie algorithm. The figure illustrates that oscillation amplitude and variability as well as the size of fluctuations are qualitatively the same with both algorithms.

The robustness of the oscillations of Fig. 3 A to variations of model parameters was assessed in the same manner as for oscillations without fluctuations (Fig. 2 A). Changes of ±20% were made in each parameter, and the resulting oscillation patterns were examined. Oscillations were preserved in all cases, and their periods and amplitudes were not very different from those of Fig. 3 A. Mean periods with standard deviations varied from 21.0 ± 0.85 h to 25.7 ± 1.0 h. The average of PER varied between 57 and 84 molecules.

Entrainment by light

Oscillations are preserved when positive feedback is eliminated

Because regulation of the rate of synthesis of the transcriptional activator dCLOCK is central to the positive feedback loop in the Drosophila model, simulations were carried out with the total concentration of dCLOCK fixed, eliminating the regulation responsible for the positive feedback. With [dCLOCK] fixed, the negative feedback loop still operates. PER still inhibits its own synthesis by binding to dCLOCK and blocking transcriptional activation.

With [dCLOCK] fixed, the Drosophila model reduces to a single delay differential equation for the rate of change of [PER]. This equation is simply Eq. 1, with the time delay ₁. [dCLOCK] is a parameter in this equation. PER synthesis is described by Eq. 2. For clarity these equations are repeated here:

(6)

(7)

with [dCLOCK_free] = [dCLOCK]  [PER] or zero, whichever is greater. The model of Eqs. 6 and 7 differs from earlier models of the Drosophila circadian oscillator that incorporated only negative feedback (Gonze et al., 2000; Leloup and Goldbeter, 1998; Goldbeter, 1995). Those models assumed PER directly repressed its own synthesis, and they did not incorporate time delays. The model of Eqs. 6 and 7 is, however, similar to that of Lema et al. (2000).

A standard parameter value set was chosen for Eqs. 6 and 7 such that all parameter values except [dCLOCK] are the same as the corresponding values in the model with positive feedback. [dCLOCK] was given a value close to the mean value of [dCLOCK] in the oscillations of Fig. 1 B. The resulting standard set is ₁ = 10 h, v_sP = 0.5 nM h¹, k_dP = 0.5 h¹, K₁ = 0.3 nM, and [dCLOCK] = 0.25 nM.

With these parameter values, circadian oscillations in the concentration of PER are obtained as illustrated in Fig. 4 A. The mechanism of oscillation is as follows. A rise in [PER] leads to a fall in free dCLOCK. After the delay ₁, the loss of free dCLOCK terminates the rise in [PER]. PER degrades, and free dCLOCK is regenerated. After a second delay ₁, this free dCLOCK leads to the next rise in [PER]. Thus, the oscillation period is somewhat more than twice ₁, with the excess depending on the speed of changes in [PER] (as determined by the parameters v_sP and k_dP).

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Simulation of circadian oscillations after positive feedback has been removed by fixing [dCLOCK]. The model of Eqs. 6 and 7 results. (A) Deterministic simulation. Blue, dark green, and light green time courses are, respectively, for [PER], [dCLOCK], and [dCLOCKfree] not bound to PER. (B) Stochastic simulation. (C) Robustness of oscillations to parameter variation. The scatter plot displays the periods and amplitudes of circadian oscillations in [PER]. The control point is at the intersection of the red lines.

Stochastic simulations without positive feedback, but including fluctuations in the numbers of PER and dCLOCK, were also carried out with Eqs. 6 and 7. The fixed time step algorithm was used. To convert concentration units to molecule numbers, the parameters v_sP and K₁ were multiplied by a factor of 250. To model fluctuations in dCLOCK, synthesis and degradation of dCLOCK needed to be represented. The average synthesis rate of dCLOCK was set to 37.5 h¹. The deterministic rate constant for degradation of dCLOCK was set to 0.5 h¹, the same as in the simulation of Fig. 3 A. The average number of dCLOCK molecules was thereby set to 75 (the ratio of the synthesis rate to the degradation rate constant). This is close to the average value in Fig. 3 A. Circadian oscillations in PER resulted (Fig. 4 B). The amplitude, period, and standard deviation of the period were similar to those of the PER oscillations in Fig. 3 A. Over 50 oscillations, the mean PER copy number was 60 and the period was 23.3 ± 1.0 h.

The simulation of Fig. 4 B illustrates that when positive feedback is eliminated, the amplitude of oscillations is only slightly diminished, and the variability of oscillation period or amplitude is not increased (compare Fig. 4 B and Fig. 3 A). With the model of Eqs. 6 and 7, entrainment to light pulses could also be simulated (not shown). Therefore, because positive feedback does not appear essential for simulating large-amplitude circadian oscillations or for simulating entrainment to light, what function can be attributed to the positive feedback loops? Perhaps positive feedback increases the robustness of oscillation amplitude and period to modest variations in parameters. To test this hypothesis, scatter plots of amplitude and period for different parameter variations were constructed, in the same manner as was done for Fig. 2 A. Fig. 4 C displays the plot obtained when the parameters in the model without positive feedback (Eqs. 6 and 7) were varied by ±20%. The only significant effect on oscillation period occurs for variations in ₁. The variability in period and amplitude is not significantly larger in the absence versus the presence of positive feedback (compare Fig. 4 C with Fig. 2 A). Therefore, these simulations do not support the hypothesis that positive feedback significantly enhances the robustness of oscillations to parameter variations.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Circadian oscillations in Drosophila can be simulated by a reduced model incorporating feedback and time delays

A model (Fig. 1 A) with only two delay differential equations is able to represent important biochemical elements of the circadian rhythm generator in Drosophila. These elements are 1) time delays to represent the intervals between changes in the concentrations of proteins that regulate transcription and changes in the rates of appearance of gene products and 2) positive and negative feedback loops underlying the regulation of gene expression. The Drosophila model simulates circadian oscillations of core gene product concentrations (Fig. 1 B). The oscillation period and amplitude do not undergo large changes given modest (20%) variations in parameter values (Fig. 2 A). The oscillations can be entrained to simulated light pulses. The photic PRC simulated by the present model shares qualitative similarities with experimental curves (Fig. 2 B), which was not the case for the PRC simulated by our previous, detailed model (Smolen et al., 2001). These results suggest that despite their simplicity, the models capture important features of the processes underlying circadian rhythmicity.

In Fig. 1 B, the oscillations of dCLOCK and PER are approximately antiphase. However, experimental dCLOCK oscillations lag PER oscillations by only 4-6 h (Lee et al., 1998). We found that the simulated lag between PER and dCLOCK could be reduced to ~6 h by increasing ₁ to 12 h and decreasing ₂ to 5 h. However, these changes markedly degraded the simulated Drosophila PRC. Only a narrow region of phase advance remained, centered on CT 0. For most CT, the phase shift was near zero. Therefore, the reduced Drosophila model fails to fully represent the molecular processes responsible for generating both the PRC and the phase relationship between dCLOCK and PER oscillations.

Comparison with an alternative model with a positive feedback loop

One process not represented in our model is an apparent stabilization of PER upon dimerization with TIM (Price et al., 1998). Such stabilization could constitute a posttranscriptional positive feedback loop, in which a rise in PER favors dimerization and PER stabilization. In an alternative model (Tyson et al., 1999), this positive feedback loop is essential for sustaining circadian oscillations. Positive feedback steepens each rise in PER. The negative feedback loop whereby PER represses its own transcription is also present in that model, to terminate each rise in PER. A decline in PER follows due to degradation of phosphorylated PER. In the model of Tyson et al. (1999), if positive feedback is removed, the negative feedback loop is not capable on its own of sustaining oscillations. Modeling suggests that certain conditions must be met for a biochemical negative feedback loop to sustain oscillations. Either the loop must contain a time delay (MacDonald, 1989) or the loop must contain a combination of many sequential reaction steps and highly cooperative feedback repression (Griffith, 1968). The negative feedback loop in the model of Tyson et al. (1999) does not meet either of these conditions, so it cannot sustain oscillations in the absence of the positive feedback loop. By contrast, in the model of Fig. 1 A, positive feedback plays an entirely different role: transcriptional regulation of dCLOCK synthesis rather than posttranslational regulation of PER degradation. In this model, positive feedback is not essential to oscillations (Fig. 4 A) because the time delay ₁ lies within the negative feedback loop.

Although important positive and negative feedback loops in Drosophila appear to be based on transcriptional regulation, additional experiments may characterize important feedback based on posttranscriptional regulation, such as the positive feedback loop postulated in the model of Tyson et al. (1999). Expressions such as Eqs. 2 or 4 might usefully represent such regulation, because these expressions incorporate saturation and delays, which are common elements in biochemical regulation. Therefore, the model of Fig. 1 A may be generic, in that it may be able to include additional regulation with only minor changes (additional terms similar to Eq. 4). Circadian regulation of per mRNA stability has been reported (So and Rosbash, 1997), but it is not known whether this regulation lies within a feedback loop.

Robust oscillations can be simulated with low (<100) average molecule numbers

It was important to develop stochastic variants of the models of Figs. 1 A and 4 A, because simplified stochastic models have proven useful in understanding the qualitative dynamics of biochemical systems with small average molecule numbers. For example, such models have illustrated that stochastic fluctuations will usually lessen, sometimes greatly, the steepness of zero-order ultrasensitivity functions (Berg et al., 2000). In some model systems, fluctuations amplify the sensitivity of enzyme activity or gene expression to changes in effector levels (Paulsson et al., 2000). We carried out simulations (Fig. 3) to determine how large the average numbers of PER and dCLOCK molecules needed to be to sustain circadian oscillations without large random variations in period or amplitude. Simulations using the Gillespie algorithm or a fixed time step algorithm (see Methods) gave qualitatively the same results (Fig. 3 B). Robust oscillations (Fig. 3 A) and entrainment to light pulses (Fig. 3 C) were simulated with relatively low molecule numbers (less than 80 for each species averaged over 50 oscillations).

Simulations suggest that time delays, but not positive feedback, are essential to generate circadian oscillations

With the Drosophila model, a time delay of hours (₁ in Eq. 2) is needed to sustain oscillations. A second delay (₂ in Eq. 4) is not necessary for oscillations but is needed to create a lag of several hours between each rise of [PER] and the following rise of [dCLOCK]. Such lags are observed in Drosophila (Lee et al., 1998).

It has been suggested that both positive and negative feedback are necessary for sustaining circadian oscillations (Hastings, 2000; Crosthwaite et al., 1997). Hastings (2000) suggested that an oscillator based on a single negative feedback loop would progressively dampen. Without positive feedback, the Drosophila model reduces to a model with a single differential equation, in which the total concentration of dCLOCK is fixed (Eqs. 6 and 7). Stable circadian oscillations of PER were simulated without positive feedback (Fig. 4, A and B). These simulations suggest positive feedback is not required to sustain circadian oscillations. An experimental prediction follows. A transgenic Drosophila line, based on dclock-null mutant animals, might be constructed in which dclock expression in neurons responsible for circadian rhythm generation is driven by a transgenic promoter expressed constitutively in those neurons. Then circadian oscillations of PER should still be evident in these neurons, as long as the level of dCLOCK is similar to the average level during oscillations in normal animals.

We also found that with or without positive feedback, oscillation period and amplitude do not vary by large amounts when modest variations are made in all parameter values (Fig. 4 C). Therefore, modeling fails to support the hypothesis that positive feedback is necessary to sustain circadian oscillations that are robust to modest variations in parameters. Our simulations also fail to provide support for the suggestion that positive feedback increases robustness to stochastic fluctuations (Hastings, 2000). This failure is evident from Fig. 4 B, in which a robust oscillation in PER is sustained with a low (<70) average molecule number. The fluctuations do not appear larger than in the analogous simulation with positive feedback (Fig. 3 A). Finally, eliminating positive feedback does not significantly reduce the amplitude of PER oscillations (compare Fig. 4, A and B, with Figs. 1 B and 3 A). Thus, our simulations do not support the recent suggestion (Cheng et al., 2001) that a role of positive feedback is to increase the amplitude of oscillations.

In what way might positive feedback be important? One possibility is that positive feedback is required to regulate output, or clock-controlled, genes (CCGs). CCGs are not part of the core feedback loops, but they are responsible for circadian variation in behaviors such as locomotion. In Drosophila, the positive feedback loop appears essential to drive circadian oscillations in the level of total dCLOCK. Positive feedback may therefore be essential to drive circadian oscillations in the expression of CCGs regulated by dCLOCK. Microarray assays at time points from CT 0 to CT 20 have recently identified more than 100 Drosophila CCGs, the majority of which are regulated (directly or indirectly) by dCLOCK (McDonald and Rosbash, 2001).

It is likely that future study of the Drosophila oscillator will reveal additional components. For example, an essential Drosophila core gene, vrille, encodes a transcription factor of as yet unknown function (Blau and Young, 1999). Therefore, the detailed descriptions of the positive and negative feedback loops within these oscillators may change. However, it is likely that reduced models, based on two or three differential equations, will remain important for representing essential mechanistic elements to aid intuitive understanding. Positive or negative feedback involving transcriptional activators and repressors also characterizes circadian rhythms in other organisms. In mammals, there are positive and negative feedback loops based on interactions of the transcriptional activator CLOCK with isoforms of PER and/or cryptochrome proteins (Shearman et al., 2000). In the cyanobacterium Synechococcus, the transcriptional activator KaiA appears to enhance the expression of the transcriptional repressors KaiC and/or KaiB (Iwasaki and Dunlap, 2000). Therefore, reduced models similar to ours, with minimal representations of the essential feedback interactions, might be useful for gaining intuitive understanding of circadian rhythm generation in Synechococcus, mammals or other organisms.