This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal 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 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.

Simulation of Drosophila Circadian Oscillations, Mutations, and Light Responses by a Model with VRI, PDP-1, and CLK

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; and ^* Department of Biology and Biochemistry, University of Houston, Houston, Texas

Correspondence: Address reprint requests to John H. Byrne, Dept. of Neurobiology and Anatomy, W. M. Keck Center for the Neurobiology of Learning and Memory, The University of Texas-Houston Medical School, PO Box 20708, Houston, TX 77225. Tel.: 713-500-5602; Fax: 713-500-0623; E-mail: john.h.byrne{at}uth.tmc.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT AND NUMERICAL... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
A model of Drosophila circadian rhythm generation was developed to represent feedback loops based on transcriptional regulation of per, Clk (dclock), Pdp-1, and vri (vrille). The model postulates that histone acetylation kinetics make transcriptional activation a nonlinear function of [CLK]. Such a nonlinearity is essential to simulate robust circadian oscillations of transcription in our model and in previous models. Simulations suggest that two positive feedback loops involving Clk are not essential for oscillations, because oscillations of [PER] were preserved when Clk, vri, or Pdp-1 expression was fixed. However, eliminating positive feedback by fixing vri expression altered the oscillation period. Eliminating the negative feedback loop in which PER represses per expression abolished oscillations. Simulations of per or Clk null mutations, of per overexpression, and of vri, Clk, or Pdp-1 heterozygous null mutations altered model behavior in ways similar to experimental data. The model simulated a photic phase-response curve resembling experimental curves, and oscillations entrained to simulated light-dark cycles. Temperature compensation of oscillation period could be simulated if temperature elevation slowed PER nuclear entry or PER phosphorylation. The model makes experimental predictions, some of which could be tested in transgenic Drosophila.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT AND NUMERICAL... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Circadian rhythms in physiology and behavior depend on the oscillating expression of genes, a few of which act as core clock components. The first core clock components identified in Drosophila were per and tim. Per and tim are activated by a heterodimer of the transcription factors CLK and CYC (Bae et al., 2000; Lee et al., 1999; Allada et al., 1998; Darlington et al., 1998). PER and TIM repress per and tim transcription, forming a negative feedback loop. This repression arises from binding of PER and TIM or of PER alone to the CLK-CYC heterodimer, preventing activation of per and tim (Bae et al., 2000; Rothenfluh et al., 2000; Lee et al., 1999). Positive feedback is also present. Elevated CLK represses Clk (Glossop et al., 1999). PER indirectly activates Clk by binding CLK and blocking this repression (Glossop et al., 1999; Bae et al., 1998). CLK's activation of per then closes a positive feedback loop. Another core clock component, vri (or vrille), is activated by CLK-CYC (Blau and Young, 1999) and in turn represses Clk (Glossop et al., 2003). Par Domain Protein 1 (Pdp-1) is another core clock component that activates Clk and is activated by CLK-CYC (Cyran et al., 2003).

Modeling is important for gaining understanding of the dynamics of gene systems containing complex regulatory motifs, such as multiple feedback loops (Hasty et al., 2001; Smolen et al., 2000). Previous models have considered the negative feedback loop based on per and tim repression (Leloup and Goldbeter, 1998; Goldbeter, 1995) and a few models have considered regulation of Clk as well as Clk's activation of per (Smolen et al., 2002, 2001; Ueda et al., 2001). However, no model has yet considered regulation of vri expression, VRI's repression of Clk, and the role of Pdp-1. To encapsulate current understanding and to make predictions that can guide further experiments, we believe it is timely to develop a model that represents the mechanisms of circadian transcriptional regulation as currently understood in Drosophila.

Fig. 1 A schematizes our model's representation of the transcriptional feedback loops. On the left is the per negative feedback loop, in which CLK activates per expression and PER then represses Clk. The TIM gene product is not represented in the model. The rationale for this simplification is discussed below (Model Development and Numerical Methods). On the right is the vri negative feedback loop, in which CLK activates vri and VRI then represses Clk. The positive feedback loop discussed above, involving Clk repression, subsumes both these negative feedback loops as follows. When Clk is activated, CLK levels increase. Activation of per by CLK results in PER synthesis and the binding of CLK by PER. Thus, activation of vri by CLK is diminished, and VRI levels fall. Finally, repression of Clk by vri is relieved, and CLK increases further. On the top of Fig. 1 A, the reciprocal activation of Clk by the Pdp-1 gene product and vice versa forms a second positive feedback loop. Fig. 1 A emphasizes the ways in which all the feedback loops include regulation of Clk expression or CLK function.

View larger version (19K): [in this window] [in a new window]   FIGURE 1  Schematic of the model for Drosophila circadian rhythm generation. (A) Feedback loops of transcriptional regulation. PER interacts with CLK, inhibiting CLK's activation of per and reducing the level of PER, forming a negative feedback loop. Vri is activated by CLK, and VRI in turn represses Clk. A second negative feedback loop is thus formed. A positive feedback loop also exists in which Pdp-1 is activated by CLK, and PDP-1 in turn activates Clk. A second positive feedback loop is formed from CLK, PER, and VRI as follows: Increased PER binds with CLK and suppresses CLK's activation of vri. The level of VRI then falls and Clk is thereby de-repressed. Then, more CLK is synthesized and the expression of per is further increased. (B) The role of PER in more detail. PER undergoes two cytosolic phosphorylations and then enters the nucleus. PER then interacts with CLK, suppressing CLK's activation of per. Nuclear PER undergoes further phosphorylations before degradation. All phosphorylation states of nuclear PER are assumed competent to interact with CLK (as illustrated by dashed box).

Etchegaray et al. (2003) reported that mammalian per genes exhibit circadian rhythms in acetylation of histone proteins, synchronous with rhythms of per mRNA. CRY proteins inhibit CLK-activated transcription, suggesting CRY inhibits histone acetylation and transcriptional activation by interacting with CLK. Our model assumes Drosophila CLK similarly induces histone acetylation, activating transcription of per, vri, and Pdp-1. PER, rather than CRY, is assumed to inhibit acetylation and repress transcription. Previous models (Ueda et al., 2001; Leloup and Goldbeter, 1998; Goldbeter, 1995) have found it essential to describe activation of per expression with nonlinear Hill functions of [CLK] or [PER]. Hill coefficients of 3–5 were used to create steep, sigmoidal functions of transcription rate versus [CLK] or [PER]. Requirements for CLK or PER to bind to multiple E-box enhancer elements in the per promoter region could justify nonlinear Hill functions. However, a promoter fragment containing a single E-box, coupled to a per transgene, is sufficient to drive robust per mRNA cycling and rescue behavioral rhythmicity in per⁰¹ mutants (Hao et al., 1999). Therefore, we suggest that the kinetics of multiple histone acetylation events are more likely to generate the steep, sigmoidal relationships of per expression rate versus levels of CLK and PER that modeling suggests is necessary for robust circadian oscillations. We developed equations to describe these sigmoidal relationships. In the resulting model, as in previous models (Ueda et al., 2001; Leloup and Goldbeter, 1998; Goldbeter, 1995), we find that if sigmoidicity is removed, circadian oscillations can no longer be simulated.

Our model simulates oscillations in constant darkness, photic entrainment of oscillations, the photic phase-response curve, and temperature compensation of oscillation period. The model also simulates the effects of null mutations of per and Clk, and the effects of vri, Clk, and Pdp-1 heterozygous null mutations on oscillation period. Simulations suggest the negative feedback loop in which PER interacts with CLK to inhibit per expression is essential for circadian oscillations of gene expression. Simulations also suggest the positive feedback loop involving Clk repression is not essential for oscillations in the expression of per, vri, or Pdp-1. However, this feedback loop is likely to play an essential role in driving oscillations in Clk expression and in regulating many clock-controlled genes regulated by CLK (McDonald and Rosbash, 2001). Therefore, this feedback loop is likely to mediate behavioral aspects of rhythmicity. Simulations also suggest the second positive feedback loop of reciprocal Pdp-1 and Clk activation may not be essential for circadian oscillations of the expression of genes other than Pdp-1.

The model makes experimental predictions. Expression of a reporter gene with a promoter including CLK binding sites is predicted to depend steeply on levels of PER and/or CLK. Circadian oscillations of per expression are predicted to be preserved with constitutive Pdp-1 or Clk expression. Either the rate of PER nuclear entry or that of PER phosphorylation is predicted to decrease with temperature.

	MODEL DEVELOPMENT AND NUMERICAL METHODS

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT AND NUMERICAL... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
For simplicity, all concentrations are referenced to the total cell volume. Absolute concentrations of circadian proteins in Drosophila neurons have not been determined, although relative abundances of PER, CLK, CYC, and TIM have been quantified (Bae et al., 2000). We chose parameters such that the average concentration of PER during an oscillation has a plausible value, 2 nM. In a Drosophila lateral neuron with a radius of 5–6 µm (Ewer et al., 1992), 2 nM would correspond to roughly 1000 PER molecules.

Circadian expression of per, vri, and Pdp-1 is activated by binding of CLK-CYC heterodimers to E-box sequences (CACGTG) in their promoter regions (Cyran et al., 2003; Blau and Young, 1999; Hao et al., 1999). The concentration of CYC in head extract is relatively constant and more than two orders-of-magnitude above the concentration of CLK, and most CLK in head extract interacts stably with CYC (Bae et al., 2000). Therefore, it is reasonable to assume that the level of CLK-CYC is limited by [CLK], and to model transcriptional activation as a function of [CLK].

Etchegaray et al. (2003) provide evidence that mammalian CLK regulates transcription of circadian genes by promoting histone acetylation. CLK/BMAL1 heterodimers activate the expression of three per genes (Per1–3) and two cryptochrome genes (Cry1–2). CRY proteins inhibit CLK-activated transcription. The promoter regions of Per1, Per2, and Cry1 exhibit circadian rhythms in acetylation of H3 histone proteins and RNA polymerase II binding. These rhythms are synchronous with those of Per1, Per2, and Cry1 mRNA. The histone acetyltransferase p300 precipitates with CLK in a circadian manner. CRY proteins inhibit the p300-induced increase in CLK-activated transcription. These results suggest CLK interacts with p300 to promote histone acetylation, chromatin structure alteration, and transcriptional activation, whereas CRY interacts with CLK-p300 to inhibit acetylation and reduce transcription. It is plausible that regulation of circadian transcription by Drosophila CLK and PER also involves histone acetylation by p300 histone acetyltransferase. In Drosophila, CRY proteins do not appear important for repressing CLK-activated transcription. Instead, the interaction of PER with CLK may reduce histone acetyltransferase activity, allowing histone deacetylation. We developed a phenomenological representation of multiple histone acetylations that generates steep, sigmoidal relationships of per, vri, and Pdp-1 expression rates versus [CLK] and [PER]. As discussed in the Introduction, previous models have suggested that such nonlinear relationships are essential for robust circadian oscillations.

We introduce first-order rate constants to describe acetylation of amino acid residues in histone proteins near E-box elements in the vri, per, and Pdp-1 promoters. These rate constants are respectively denoted k_Vacet, k_Pacet, and k_PDacet. Acetylation is assumed to require the binding of CLK to an E-box, and to be repressed by binding of PER. The acetylation rate constants are therefore represented as a product of hyperbolic Hill functions of [CLK] and of [PER]. The latter function describes repression by nuclear PER with concentration [PER_nuc]. The following expressions show the result:

(1)

(2)

(3)

In Eqs. 1–3, F_V, F_P, and F_PD are arbitrary scaling factors. Corresponding rate constants for deacetylation at the vri, per, and Pdp-1 promoters are introduced. These are denoted k_Vdeac, k_Pdeac, and k_PDdeac, and are assumed to be fixed. This assumption is consistent with data of Etchegaray et al. (2003) illustrating nonrhythmic histone deacetylase activity at the mammalian per1 and per2 promoters.

Variables AC_Vri, AC_Per, and AC_Pdp describe acetylation of amino acid residues in a representative histone protein at the vri, per, or Pdp-1 promoters. These variables denote the fraction of eligible residues that are acetylated. Other rhythmic histone modifications, such as phosphorylations (Etchegaray et al., 2003), may also affect transcription and could be subsumed into AC_Vri, AC_Per, and AC_Pdp. Because the number of residues eligible for acetylation or for other rhythmic modification is not known, we approximate AC_Vri, AC_Per, and AC_Pdp as continuous variables ranging from 0 to 1. First-order differential equations describe their dynamics, using acetylation rate constants defined by Eqs. 1–3 and the corresponding deacetylation rate constants:

(4)

(5)

(6)

These equations do not yet imply any sigmoidal relationships between levels of CLK or PER and transcriptional activation. To introduce sigmoidicity, we make the plausible assumption that multiple histone proteins need to be acetylated at the vri, per, and Pdp-1 promoters to bring about an "open" DNA configuration readily accessible to RNA polymerase. We describe accessibility to RNA polymerase by variables OP_Vri, OP_Per, and OP_Pdp. A requirement for acetylation of multiple histones can be represented as a nonlinear dependence of OP_Vri, OP_Per, and OP_Pdp on the single-histone acetylation fractions AC_Vri, AC_Per, and AC_Pdp. We chose a phenomenological, sigmoidal nonlinearity. Accessibilities were assumed to approach steady-state values given by a power N of the acetylation fractions (usually N = 5). These steady-state values were denoted OP_{Vri, ss}, OP_{Per, ss}, and OP_{Pdp, ss}:

(7)

(8)

(9)

OP_Vri, OP_Per, and OP_Pdp relax to these steady-state values with time constants _{Vri, op}, _{Per, op}, and _{Pdp, op}:

(10)

(11)

(12)

The expression rates of per, vri, and Pdp-1 are denoted R_Per, R_Vri, and R_Pdp. These rates are assumed proportional to the accessibility of promoters to RNA polymerase, i.e., to OP_Per, OP_Vri, and OP_Pdp. Small basal rates in the absence of CLK are denoted R_Pbas, R_Vbas, and R_PDbas. The following equations show the result:

(13)

(14)

(15)

Fig. 2 illustrates the sigmoidicity of the relationships that the above equations imply for the steady-state per expression rate versus the levels of CLK and PER. The expression rate varies from a minimum of R_Pbas (in the absence of CLK, or with high [PER_nuc]) to a maximum of V_Per + R_Pbas (for saturating [CLK] with PER absent). Graphs of the steady-state vri and Pdp-1 expression rates are very similar (not shown).

View larger version (14K): [in this window] [in a new window]   FIGURE 2  Steady-state per expression rates implied by the model's description of histone acetylation kinetics and of the effect of acetylation on transcription (Eqs. 1–15). Equations 2, 5, 8, and 13 were used to determine the expression rate as a function of [CLK] and [PERnuc], with model parameters at standard values (Table 1). (A) Log-log plot of per expression rate versus [CLK] in the absence of repression by PER ([PERnuc] = 0). (B) Log-log plot of per expression rate versus [PERnuc] with [CLK] fixed at a saturating level (2.5 nM).

(16)

The PDP-1 and VRI DNA binding domains are similar. Competition between PDP-1 and VRI has been demonstrated at one PDP-1 binding site and is thought to take place at the other sites (Cyran et al., 2003). Eq. 16 does not represent the details of this competition. Instead, Eq. 16 was constructed to qualitatively represent nonlinear activation by PDP-1 and repression by VRI. We note that activators of Clk other than PDP-1 may also play a significant role, because Clk^Jrk mutant Drosophila have low levels of PDP-1, but high levels of Clk mRNA (near the wild-type peak; Cyran et al., 2003; Glossop et al., 1999).

View larger version (16K): [in this window] [in a new window]   FIGURE 3  Simulation of circadian oscillations in constant darkness (DD). Parameters are set to standard values (Table 1). (A) Time courses of [VRI], [PERtot], [PDP-1], and [CLK]. (B) Time courses of PER cytosolic and nuclear species concentrations [P0cyt], [P2cyt], and [PERnuc]. For clarity, time courses of [P0cyt] and [P2cyt] are scaled vertically by a factor of 10. (C) Oscillations in [VRI], [PERtot], [PDP-1], and [CLK] when third powers of ACVri, ACPer, and ACPdp are used in Eqs. 7–9. Values of [CLK] are multiplied by 2 for ease of visualization. (D) Periods and amplitudes of simulated circadian oscillations in [PERtot]. To generate these oscillations, each parameter in the set of standard parameter values (Table 1, excluding N and klight) was increased or decreased by 20%. There are 77 data points including the control with all parameter values standard. The control is at the intersection of the horizontal and vertical red lines.

View larger version (19K): [in this window] [in a new window]   FIGURE 7  Simulated effects of light exposure and simulated temperature compensation. All parameters are at standard values except as specified. (A) Entrainment of oscillations to a 22-h light-dark cycle. The overbar indicates the dark phase. Light is assumed to increase the degradation of cytosolic PER. During the light phase, a first-order degradation term was added to the right-hand side of the differential equations for [P0cyt], [P1cyt], and [P2cyt]. The first-order rate constant was 0.685 h–1. (B) Photic phase-response curves (PRCs). When constructing the model PRC (solid curve), each light pulse was simulated by adding, for 2 h, a first-order degradation term to the differential equations for each cytosolic form of PER. The first-order rate constant was 0.685 h–1. An experimental Drosophila PRC (Fig. 5 of Konopka et al., 1991) is also illustrated (circles). The means of the experimental phase shifts are displayed. (C) Temperature-compensated oscillations. The standard parameter value set (Table 1) was used as the starting point. Parameters containing units of time (h) or inverse time (h–1) were respectively divided or multiplied by a Q10 of 2 to simulate a temperature increase of 10°C. However, the maximal velocity for PER nuclear transport (vtrans) was assumed to decrease with increasing temperature, and was multiplied by a Q10 of 0.62. A period of 24.2 h results. Values of [CLK] are multiplied by a factor of 3 for ease of visualization.

(17)

(18)

The differential equation for [PDP-1] is similar, but includes a delay of several hours between regulation of Pdp-1 as given by Eq. 15 and regulation of PDP-1 synthesis. Pdp-1 and vri are both activated by CLK, but during a circadian oscillation, the rise in [PDP-1] lags the rise in [VRI] by 3–4 h (Cyran et al., 2003). The mechanism underlying this delay is not known. To implement this delay, the rate of PDP-1 synthesis, R_Pdp in Eq. 15, was continuously computed and stored. The stored values were used to calculate the rate of PDP-1 synthesis 3 h later (_delay in Table 1). The following equation describes these dynamics, with brackets denoting the discrete time delay:

(19)

PER is progressively phosphorylated over a period of hours after its synthesis (Edery et al., 1994). The model assumes sequential phosphorylations, as have previous models (Smolen et al., 2001; Leloup and Goldbeter, 1998; Goldbeter, 1995). The number of phosphorylations and their distribution between cytosol and nucleus is not known. The model represents four phosphorylations of PER, 2 of which are cytosolic. Nonphosphorylated cytosolic PER is denoted P0_cyt, and its synthesis rate is given by Eq. 13. One phosphorylation gives P1_cyt, and a second yields P2_cyt. For simplicity, phosphorylations are assumed to be irreversible and described by Michaelis-Menten rate expressions. The Michaelis constant is denoted K_pcyt, and the maximal velocity is denoted v_pcyt. P2_cyt is removed by transport into the nucleus, which is described by a Michaelis constant of K_pcyt and a maximal velocity of v_trans. These assumptions yield the following differential equations for [P0_cyt], [P1_cyt], and [P2_cyt]:

(20)

(21)

(22)

Nuclear PER phosphorylation states are denoted P0_nuc, P1_nuc, and P2_nuc. P0_nuc is formed from P2_cyt by transport into the nucleus. The resulting differential equations for [P0_nuc] and [P1_nuc] are

(23)

(24)

P2_nuc is formed by phosphorylation of P1_nuc. Experimentally, PER is seen to degrade relatively rapidly after its phosphorylation (Edery et al., 1994). In the model, degradation of P2_nuc is described by a Michaelis-Menten term with maximal velocity v_degp and Michaelis constant K_degp. The differential equation for [P2_nuc] is therefore

(25)

The relative peak-to-trough amplitude (% difference) of simulated circadian oscillations was found to be significantly enhanced if phosphorylation and degradation of nuclear PER were modeled as approximately zero-order processes. Therefore, the corresponding Michaelis constants (K_pnuc and K_degp in Eqs. 23–25) are assigned very low values (0.001 nM and 0.01 nM, respectively).

Certain combinations of molecular species play a significant role in the dynamics of the model, or display time courses that can be compared with experimental data. The total concentration of nuclear PER is given as

(26)

The model dynamics are fully determined by Eqs. 1–26.

The concentration of cytosolic PER is given by an expression analogous to Eq. 26,

(27)

Because concentrations are referenced to the total cell volume, the total concentration of PER, [PER_tot], is the sum of the nuclear and cytosolic concentrations

(28)

The above model does not represent expression of tim or the dynamics of the TIM gene product. Recent data suggests monomeric nuclear PER is the species responsible for inhibiting CLK's activation of per, tim, and other core circadian genes (Weber and Kay, 2003; Rothenfluh et al., 2000). Thus, TIM does not appear to directly regulate transcription. However, TIM indirectly regulates transcription by regulating the level and localization of PER. Tim null homozygotes are arrhythmic, and [PER] is decreased (Price et al., 1995). PER nuclear localization is blocked in homozygotes (Vosshall et al., 1994). PER and TIM heterodimerize, and it has been suggested that heterodimer formation suffices for PER nuclear entry (Saez and Young, 1996). However, Shafer et al. (2002) report PER accumulates in the nucleus of ventrolateral neurons during early night, whereas most TIM remains cytosolic. These data suggest a mechanism in which TIM interacts with an unidentified cytosolic factor and blocks its ability to retain PER in the cytosol, as first noted by Vosshall et al. (1994). Alternatively, PER-TIM heterodimerization might be required for nuclear entry, but during the early night, TIM may be re-exported (Shafer et al., 2002). Because the mechanism of TIM's facilitation of PER nuclear entry is uncertain, and because TIM appears not to directly regulate transcription, it appears reasonable not to explicitly model TIM dynamics. The number of equations is thereby reduced considerably. However, in future work we plan to consider the effect of cyclical modulation of PER nuclear entry by TIM.

For most simulations, a standard set of model parameter values was used. Table 1 divides the parameters into groups that describe similar processes and displays the standard values. As noted in the footnote to Table 1, these values include a scaling factor to adjust the free-running oscillation period to 24 h.

Data to constrain values of kinetic parameters are generally lacking. To obtain standard parameter values, it was necessary to rely on trial-and-error variation. Values were found that generated simulations similar to experimental data for the following phenomena: stable circadian oscillations robust to small parameter changes, simulation of entrainment to light pulses and simulation of a photic phase-response curve, and oscillations and steady states of circadian gene mutants. For example, in Table 1, the maximal rates of gene expression (V_Clk, V_Vri, V_Per, and V_Pdp in Eqs. 13–16) cover a broad range. These values were chosen via extensive efforts to simulate both normal circadian oscillations and steady states of mutants (Fig. 3 A and Fig. 4, A–C).

View larger version (17K): [in this window] [in a new window]   FIGURE 4  Simulations of clock gene mutants. Parameters are at standard values (Table 1) except as specified for each simulation. In Fig. 4 and subsequently, time courses of [CLK], [VRI], [PERtot], and [PDP-1] are respectively solid black, dashed black, solid shaded, and dashed shaded. (A) A per01 homozygous mutation. PER is assumed unable to interact with CLK or regulate transcription. This is implemented by setting the Hill functions of [PERnuc] to 1 in Eqs. 1–3. (B) A ClkJrk homozygous mutation. CLK is assumed unable to activate transcription. This is implemented by setting the Hill functions of [CLK] to zero in Eqs. 1–3. (C) Decrease in period with increased PER phosphorylation. The maximal phosphorylation velocities and the nuclear transfer velocity (vpcyt, vpnuc, and vtrans) were increased by 65% from their standard values.

Numerical methods
The forward-Euler method was used for integration of differential equations. Integration time steps were reduced until no significant difference was seen upon further reduction (2 x 10^–4 h). The model was programmed in Java and simulated on a Pentium 3 microcomputer. Simulation programs are available from the authors upon request.

	RESULTS

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT AND NUMERICAL... RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Simulated oscillations of Drosophila core gene product levels resemble experimental oscillations
Fig. 3 A illustrates that the model of Fig. 1 simulates large-amplitude circadian oscillations in the concentrations of CLK, VRI, PDP-1, and total PER (cytosolic plus nuclear, [PER_tot] in Eq. 28) with a period of 24 h. Effects of light are not simulated, so these oscillations correspond to a free-running rhythm in constant darkness (DD). For [PER_tot], the peak and average are approximately twofold higher than for [CLK]. The peak of [VRI] is also higher than the peak of [PER_tot], and the [VRI] peak leads the [PER_tot] peak by 5 h. In these respects, the simulation resembles experimental data (Glossop et al., 2003; Bae et al., 2000; Lee et al., 1998). The delay between the accumulation of [VRI] and that of [PER_tot] is due to the time required for the multiple phosphorylation states of PER to accumulate. Fig. 3 A displays a lag of 3 h between the rise of [VRI] and that of [PDP-1]. A similar lag is seen experimentally (Cyran et al., 2003). Because CLK activates both vri and Pdp-1 expression, this lag could only be simulated by introducing a time delay (_delay in Eq. 19) between CLK's regulation of Pdp-1 expression and PDP-1 protein accumulation.

Fig. 3 B illustrates oscillations in nuclear [PER] concentration ([PER_nuc]), which is assumed to mediate repression of CLK-driven transcription (Eqs. 1–3). Fig. 3 B also illustrates oscillations of newly synthesized cytosolic PER ([P0_cyt]), and doubly phosphorylated cytosolic PER ([P2_cyt]). There is a substantial delay (6 h) between the peak of [P0_cyt] and that of [PER_nuc], due to the time required for cytosolic PER to undergo two phosphorylation events and then to be transported into the nucleus (Eqs. 20–23). After PER enters the nucleus, it is assumed to undergo two more phosphorylations and then to be degraded (Eqs. 23–25). Similar dynamics are seen experimentally; multiple phosphorylations of PER over hours are followed by rapid degradation of phosphorylated PER (Edery et al., 1994).

Although per mRNA is not an explicit model variable, comparison of the time courses of [PER_tot] (Fig. 3 A) and of [P0_cyt] (Fig. 3 B) illustrates that the model represents much of the 5–6 h delay observed between the time courses of per mRNA and [PER_tot] (Glossop et al., 1999; Vosshall et al., 1994; Hardin et al., 1990). The time required for PER phosphorylation and nuclear entry results in a delay of 5 h between the peak of the simulated time course of the initial form of PER, [P0_cyt] (at t = 11 h), and the peak of [PER_tot] (at t = 16 h). A similar delay is seen between the initial rise of [P0_cyt] and the rise of [PER_tot].

The mechanism underlying the oscillations in Fig. 3, A and B, can be summarized as follows. At time t = 4 h, [CLK] is rising close to its peak level, and [PER_tot] and [PER_nuc] are falling. These changes activate histone acetylation at the vri, per, and Pdp-1 promoters (Eqs. 1–3) and transcription of vri, per, and Pdp-1 (Eqs. 13–15). [VRI] begins to rise at t = 7 h. The continuing decay of nuclear PER from the previous cycle of PER accumulation prevents [PER_tot] from rising until t = 10 h, and the delay _delay between Pdp-1 transcription and PDP-1 synthesis prevents [PDP-1] from rising until t = 11 h. [PER_nuc] rises slightly after [PER_tot] (Fig. 3 B). From t = 9 h to t = 17 h, [VRI] is high, Clk is repressed, and [CLK] is falling. Elevated [PER_nuc] promotes histone deacetylation by decreasing the acetylation rate constants (Eqs. 1–3). Therefore vri, per, and Pdp-1 transcription are inhibited after t = 13 h, when [PER_nuc] is high. Because of this inhibition, [VRI], [PER_tot], [PER_nuc], and [PDP-1] pass in sequence through their peaks during t = 12–16 h and then decline. By t = 24 h [VRI], [PER_tot], and [PDP-1] have dropped well below their peak levels. The fall in [VRI] has removed its inhibition of Clk, and [CLK] is therefore rising again, to begin another cycle.

There is considerable variation in standard values for the basal synthesis rates of CLK, PER, VRI, and PDP-1 (R_Cbas, R_Pbas, R_Vbas, and R_Pdbas, respectively, in Table 1). The values vary between R_Cbas = 0.001 nM h^–1 and R_PDbas = 0.36 nM h^–1. We found empirically that, for simulating circadian oscillations (Fig. 3 A), a small R_Cbas was required to make the trough of [CLK] oscillations near zero. This, in turn, was required to sustain large peak-to-trough ratios of oscillations in [PER_tot], [VRI], and [PDP-1]. In contrast, the peak of [CLK] occurs just before the trough of [PDP-1] (Fig. 3 A). To maintain substantial CLK synthesis and a significant [CLK] peak, it was necessary that the [PDP-1] trough not be too low. Therefore R_PDbas had to be relatively high.

As discussed in Model Development and Numerical Methods, activation of vri, per, and Pdp-1 expression is proportional to the variables OP_Vri, OP_Per, and OP_Pdp, which represent promoter accessibilities to RNA polymerase (Eqs. 7–15). A putative requirement for acetylation of multiple histones is represented by a sigmoidal, "threshold" dependence of these promoter accessibilities on the variables AC_Vri, AC_Per, and AC_Pdp, which represent acetylation of single histones. The promoter accessibilities relax to steady-state values given by fifth powers of AC_Vri, AC_Per, and AC_Pdp (Eqs. 7– 9). Simulations suggested this sigmoidicity is essential. If only first powers of AC_Vri, AC_Per, and AC_Pdp were used in Eqs. 7–9, circadian oscillations could not be simulated irrespective of the values for the parameters in Table 1. Oscillations could be sustained if third or fourth powers were used. However, Fig. 3 C illustrates that with third powers, the peak-to-trough ratio of oscillations is considerably reduced. Experimental time courses (Cyran et al., 2003; Glossop et al., 2003; Bae et al., 2000; Lee et al., 1998) exhibit much larger peak-to-trough ratios, more similar to Fig. 3 A.

Circadian rhythms of histone acetylation at mammalian per1 and per2 promoters appear synchronous with the rhythms of per1 and per2 mRNA (Etchegaray et al., 2003, their Fig. 1). Do our simulated oscillations display similar synchrony? Our model does not include mRNA concentrations. However, time courses of the per, vri, and Pdp-1 transcription rates (R_Per, R_Vri, and R_Pdp, Eqs. 13–15) can be compared with time courses of acetylation at the per, vri, and Pdp-1 promoters (AC_Vri, AC_Per, and AC_Pdp, Eqs. 4–6). For the oscillations of Fig. 3, the comparison reveals very little lag (1.5 h) between the peaks of acetylation and the following peaks of transcription rates (not shown). Therefore, this simulation appears qualitatively consistent with the data of Etchegaray et al. (2003).

The model describes Clk activation by PDP-1, and Clk repression by nuclear PER, with nonlinear functions. In Eq. 16, a Hill coefficient of 2 is assumed to describe both these processes. We repeated the simulation of Fig. 3 A with these nonlinearities diminished (Hill coefficients of 1 in Eq. 16) or enhanced (Hill coefficients of 3 in Eq. 16). The effects of these changes were minor. Oscillations were preserved, with circadian periods (23.8 and 24.3 h, respectively). The relative phases and amplitudes of [PER], [CLK], [PDP-1], and [VRI] oscillations remained similar to those in Fig. 3 A. Peak-to-trough ratios were slightly reduced with Hill coefficients of 1 (not shown).

Simulated oscillations are robust to parameter variations
Despite heterogeneity between individuals in the values of parameters describing gene expression regulation, individual Drosophila maintain very similar circadian periods in DD. For example, Bao et al. (2001) reported periods of 24.3 ± 0.06 h for 43 wild-type flies. A model of circadian rhythm generation should therefore be robust in the sense that small parameter variations should not result in large period alterations. Robustness was tested with a method used in previous investigations (Smolen et al., 2001; Lema et al., 2000; Goldbeter, 1995). Each parameter was increased or decreased by a moderate amount (20%) from its standard value in Table 1. The effects on oscillation amplitude and period were determined, with all other parameters kept at standard values. There are 38 model parameters (Table 1, not including k_light or N). Therefore, 77 simulations were carried out including the control with standard parameter values. Fig. 3 D plots the period and amplitude of the resulting oscillations. The amplitude was measured as the peak-to-trough difference in [PER_tot].

Oscillations were preserved in all simulations. The mechanism of oscillation remained as described for control oscillations (Fig. 3, A and B). Most oscillations have periods and amplitudes close to control values (24 h, 3.30 nM). Three parameter changes yielded periods differing >3 h from control. The largest period and amplitude decreases (20.3 h, 1.54 nM) occurred for the 20% decrease in v_pnuc. A 20% increase in v_pnuc produced the largest amplitude increase (5.50 nM) and a significant period increase (27.5 h). Increasing v_pnuc acts to inhibit and delay accumulation of PER_nuc by promoting the formation of the most highly phosphorylated species of PER, P2_nuc, which is rapidly degraded. The inhibition of PER_nuc accumulation enhances CLK's activation of per, vri, and Pdp-1 expression. Therefore, [PER_tot] reaches a higher peak and takes longer to degrade. The largest period increase (28.3 h) occurred for a 20% decrease in K_PP (Eq. 2). Perhaps uncharacterized mechanisms exist in Drosophila to reduce the sensitivity of oscillations to v_pnuc and K_PP. Overall, the lack of large period or amplitude changes suggests the model is sufficiently robust to be regarded as a reasonable representation of biochemical mechanisms responsible for circadian oscillations in Drosophila.

Simulations of mutations in core clock genes, and of enhanced PER phosphorylation, are qualitatively similar to experimental data
Starting from standard parameter values (Table 1, Fig. 3, A and B), the model was perturbed to simulate the effects of homozygous null mutations in per and Clk. Following established nomenclature, these mutations are denoted per⁰¹ and Clk^Jrk. In these mutations, nonfunctional proteins appear to be produced. Therefore, to simulate Clk^Jrk, the transcriptional activation terms for per, Pdp-1, and vri in Eqs. 1–3 (the Hill functions of [CLK]) were set to 0. To simulate per⁰¹, the transcriptional repression terms in Eqs. 1–3 (the Hill functions of [PER_nuc]) were set to one. The per⁰¹:Clk^Jrk double mutant was also simulated, by applying both the Clk^Jrk and per⁰¹ perturbations.

Fig. 4, A and B, illustrate simulations of per⁰¹ and Clk^Jrk. Levels of [PER_tot], [VRI], [PDP-1], and [CLK] can be compared with simulated wild-type (WT) oscillations (Fig. 3 A). For per⁰¹, the PER_tot level is approximately one-third of the WT peak, qualitatively consistent with observations that per mRNA levels are 50% of the WT peak (So and Rosbash, 1997; Hardin et al., 1990). [CLK] is very low, qualitatively consistent with the observation of a low Clk mRNA level in per⁰¹ (20% of the WT peak, Glossop et al., 1999). However, [VRI] is nevertheless high, near its peak in Fig. 3 A. This apparent contradiction results because the low [PER_tot] and low [PER_nuc] result in little repression of vri histone acetylation and vri expression (Eqs. 1 and 16). Under this condition, even low [CLK] suffices to drive substantial VRI synthesis. This simulation does show some discrepancy with data, in that experimentally, vri mRNA is at intermediate, not peak, levels in per⁰¹ flies (Blau and Young, 1999).

In the Clk^Jrk simulation (Fig. 4 B), PER_tot is very low. This is qualitatively consistent with the observation that in Clk^Jrk flies, per mRNA levels are very low (<15% of the WT peak, Allada et al., 1998). In contrast, the level of [CLK] is similar to the peak of [CLK] in Fig. 3 A. This is consistent with the observation that in Clk^Jrk flies, Clk mRNA is steady and near the WT peak (Glossop et al., 1999). Because the nonfunctional CLK does not activate vri, [VRI] is low (near the trough of [VRI] in Fig. 3 A). Because VRI is not repressing Clk, a relatively high [CLK] is maintained even though the level of Clk's activator, PDP-1, is relatively low. In per⁰¹:Clk^Jrk double mutants, Clk mRNA is also observed to be near the WT peak (Glossop et al., 1999). The per⁰¹:Clk^Jrk simulation yielded a high [CLK] level, identical to Clk^Jrk. This result was expected, because after the Hill functions of [CLK] in Eqs. 1–3 are set to 0 to simulate Clk^Jrk, no further change occurs if the Hill functions of [PER_nuc] are set to 1 to simulate per⁰¹.

If any of the multiple phosphorylations of PER are obligatory for nuclear entry of PER, then these phosphorylations constitute kinetic steps within the negative feedback loop in which PER acts to repress per expression. In the model, both cytosolic phosphorylations of PER lie within this feedback loop. Therefore, increasing the PER phosphorylation rate is predicted to speed up the dynamics of this loop and shorten the oscillation period. Fig. 4 C confirms this prediction. A 65% increase was applied to the maximal phosphorylation velocities of PER and to the velocity of PER nuclear transfer (v_pcyt, v_pnuc, and v_trans, Eqs. 20–25). The period decreased to 18.6 h. The increase in v_trans was necessary to reduce the period below 19 h. It is plausible that PER nuclear transfer involves a phosphorylation, so that its rate would be affected by altering PER phosphorylation kinetics. In Fig. 4 C, [PDP-1] and [VRI] oscillations are larger than control oscillations (Fig. 3 A). This occurs because the increased v_pnuc allows more rapid conversion of nuclear PER into fully phosphorylated P2_nuc. P2_nuc rapidly degrades, so PER_nuc accumulation is inhibited. Therefore, the expression of vri and Pdp-1 is not as repressed by PER_nuc, yielding higher peaks of [PDP-1] and [VRI].

The per^S mutation shortens the free-running circadian period to 19 h (Konopka and Benzer, 1971) and results in more rapid PER phosphorylation, with a greater proportion of PER^S extensively phosphorylated at earlier Zeitgeber time (Edery et al., 1994). The simulation of Fig. 4 C may qualitatively represent these aspects of the per^S phenotype. Another mutation, dbt^S, also shortens the free-running circadian period (Price et al., 1998). The dbt gene product (DBT) is a casein kinase I homolog that phosphorylates PER (Kloss et al., 1998). Its activity may be increased in dbt^S, but this has not been biochemically confirmed. The period in Fig. 4 C is similar to that of dbt^S (18 h, Price et al., 1998). Two other mutations, dbt^g and dbt^h, increase the free-running period and may decrease DBT activity (Suri et al., 2000). However, if these mutations are simulated by a parameter change opposite to that in Fig. 4 C (decreasing v_pcyt, v_pnuc, and v_trans), the period decreases. Therefore, the model may not satisfactorily simulate dbt^g and dbt^h mutations, although a 32% decrease in v_pcyt alone does yield a period of 29 h, similar to a dbt^h homozygote or dbt^g heterozygote.

Blau and Young (1999) reported that circadian rhythms are preserved in flies having one functional copy of vri (i.e., vri^null/+), but that the period was decreased by 0.4–0.8 h. In Fig. 5 A, a vri^null/+ mutation was modeled with a 50% decrease in the maximal and basal velocities of vri transcription (V_Vri and R_Vbas, Eq. 14). The [VRI] and [PER_tot] oscillations are similar to WT oscillations (Fig. 3 A) except that the amplitude of the [VRI] oscillation is diminished. Repression of CLK synthesis is thereby relieved and [CLK] levels are above the WT control (Fig. 3 A). The period is decreased to 23.3 h, 0.7 h below WT. We also simulated heterozygous null mutations of Clk, Pdp-1, and per (Clk^Jrk/+, Pdp^null/+, per⁰¹/+) with 50% decreases in the corresponding maximal and basal transcription rates in Eqs. 16, 15, and 13 (V_Clk, R_Cbas, V_Pdp, R_PDbas, V_Per, R_Pbas). Experimental periods of these heterozygotes in DD are similar to WT: 24.8 h for Clk^Jrk/+ (Allada et al., 1998), 23.6 h for Pdp^null/+ (Cyran et al., 2003), and 25.2 h for per⁰¹/+ (Konopka and Benzer, 1971). Simulated Clk^Jrk/+ and Pdp^null/+ periods were close to the experimental periods: 25.0 h for Clk^Jrk/+, 23.9 h for Pdp^null/+. The appearance of the oscillations was similar to WT (Fig. 3 A), except the [CLK] and [PDP-1] amplitudes were respectively decreased by 50% in Clk^Jrk/+ and Pdp^null/+. In contrast, the per⁰¹/+ simulation failed to sustain oscillations, contradicting experimental data. A steady state was obtained, with [VRI] high (7.0 nM) and [CLK] and [PER_tot] very low.

View larger version (24K): [in this window] [in a new window]   FIGURE 5  Simulations with altered vri and Clk expression. (A) A vrinull/+ mutant. The VRI synthesis rate Rvri (Eq. 14) was scaled by a factor of 0.5 throughout the simulation. All species continue to show large-amplitude oscillations. The period is 23.3 h. (B) Clk overexpression. A constant value of 0.3 nM h–1 was added to the CLK synthesis rate Rclk (Eq. 16). The oscillation period is 23.4 h.

Hao et al. (1999) rescued circadian rhythmicity in per⁰¹ mutants by introducing a per transgene with a promoter more active than that of WT per. However, the rhythm had a short period (22.5 h). The model simulates a similar period decrease with per overexpression. If the maximal and basal rates of per expression (v_per and R_Pbas) are increased by 60%, the period decreases to 22.4 h.

Removal of positive feedback eliminates circadian expression of some clock genes, but not of per; negative feedback is essential for oscillations of all clock genes
The positive feedback loop relying on repression of Clk was elucidated by Glossop et al. (1999) and is summarized in the legend of Fig. 1. Constitutive vri or Clk expression can be simulated by fixing [VRI] or [CLK]. Either manipulation eliminates the positive feedback loop by eliminating dynamic repression of Clk. In previous models, fixation of [CLK] did not prevent circadian oscillations of per expression (Smolen et al., 2002, 2001). Similarly, Fig. 6 A illustrates that, when [CLK] is fixed at half of its peak (as in Fig. 3 A),oscillations in [PER_tot] persist, with period near circadian (26 h). The oscillations in [PER_tot] and [PER_nuc] result in peaks of PER_nuc that cyclically inhibit CLK's ability to activate Pdp-1, vri, and per. Therefore, [PDP-1] and [VRI] also oscillate.

View larger version (25K): [in this window] [in a new window]   FIGURE 6  Simulations with fixed Clk and vri expression. All parameters are at standard values except as specified. (A) Fixed Clk expression. [CLK] is fixed at one-half of its maximal concentration in Fig. 3 A (at 0.79 nM). Oscillations in [PERtot], [VRI], and [PDP-1] persist, with period near circadian (26.0 h). (B) Constitutive vri expression is simulated by fixing [VRI] at 3.0 nM. [PERtot] and [PDP-1] continue to oscillate with a period of 19.7 h, and oscillations in [PDP-1] drive low-amplitude oscillations in [CLK]. The plotted [CLK] values are increased 10-fold for visualization.

Experimentally, vri overexpression often leads to arrhythmicity (Blau and Young, 1999). In the model, maintaining [VRI] at a high level strongly inhibits circadian oscillations. Fixing [VRI] 60% above its level in Fig. 6 B yields an amplitude of only 0.85 nM for [PER_tot]. The period is decreased to 14.4 h. This decrease is discrepant with data, in that Blau and Young (1999) found that the free-running period of rhythmic Drosophila with vri overexpressed was lengthened by 1.5–3 h. In the model, higher [VRI] abolishes oscillations. In contrast, the dynamics in Fig. 6 A are relatively insensitive to variations in [CLK]. [PDP-1] and [PER_tot] oscillations of nearly circadian period are preserved when [CLK] is fixed at any value between 0.1 nM and infinity, with "infinite [CLK]" implemented in Eqs. 1– 3 by setting the functions of [CLK] on the right-hand sides to 1. Eqs. 1–3 therefore qualitatively model the following situation: Even with very high [CLK], only a limited and saturated amount of [CLK] is bound at E-box sites to promote histone acetylation. PER_nuc still binds cyclically to CLK–E-box complexes, repressing CLK's action and sustaining oscillations of acetylation and gene expression. Circadian oscillations of [PDP-1] and [PER_tot] are also simulated if [VRI] is fixed at a low value or at 0 (not shown).

The model contains a second positive feedback loop in which PDP-1 activates Clk and CLK reciprocally activates Pdp-1 (Fig. 1 A). To remove only this feedback loop, Pdp-1 expression was fixed by holding [PDP-1] constant. Oscillations of [VRI], [CLK], and [PER_tot] were preserved. The shapes and relative phases were similar to the control oscillations (Fig. 3 A) for a wide range of fixed [PDP-1] (from <0.3 nM to arbitrarily high levels). The period remained circadian.

The above simulations suggest neither of the positive feedback loops in Fig. 1 A may be required to sustain circadian oscillations in the levels of PER, VRI, or free CLK not complexed with PER. What about the core negative feedback loop in which PER represses per transcription by binding CLK? Simulations were carried out with this loop removed by fixing the PER synthesis rate (R_Per, Eqs. 13 and 20). Under this condition, oscillations were abolished (not shown). Only steady states of all concentrations were obtained, irrespective of the values chosen for the parameters in Table 1.

It can be hypothesized that, although the positive feedback loops in Fig. 1 A are not required for circadian oscillations per se, they may increase the robustness of oscillation amplitude and period to modest variations in parameters. To determine whether our model supports this hypothesis, we constructed scatter plots of oscillation period versus [PER_tot] oscillation amplitude, analogous to Fig. 3 C. Plots were constructed for three conditions: 1), [VRI] fixed at 3.0 nM as in Fig. 6 B, thereby removing the PER–CLK positive feedback loop and the VRI–CLK negative feedback loop; 2), [CLK] fixed at 0.79 nM as in Fig. 6 A, thereby eliminating all feedback loops except the PER–CLK negative feedback loop; and 3), [PDP-1] fixed at 3.0 nM, thereby eliminating only the CLK–PDP-1 positive feedback loop. To construct each scatter plot, each model parameter was increased or decreased by 20% from its standard value. There are 38 model parameters. Therefore, 76 simulations were carried out for each scatter plot. We found that none of the plots showed significantly less robustness to parameter variation than did the control plot with all feedback loops present (Fig. 3 C). The percent difference between minimal and maximal amplitudes or periods in each plot was not greater than in Fig. 3 C, with the single exception that decreasing v_pnuc by 20% abolished oscillations for fixed [CLK]. Therefore, these simulations fail to support the hypothesis that either the positive feedback loops in Fig. 1 A, or the VRI–CLK negative feedback loop, significantly increase the robustness of oscillations to moderate parameter variations.

The model simulates photic entrainment and a phase-response curve similar to experimental curves
In Drosophila, light enhances degradation of phosphorylated TIM (Myers et al., 1996; Zeng et al., 1996). When TIM is removed from the complex of PER and TIM, phosphorylation of PER is strongly enhanced (Kloss et al., 2001). It is plausible that accelerated degradation of highly phosphorylated PER would result, and experiments have demonstrated light-induced enhancement of PER phosphorylation and degradation (Lee et al., 1996). Therefore, we modeled light exposure by adding a first-order degradation term to the right-hand side of each differential equation for a form of cytosolic PER. The degradation rate constant was denoted k_light. For example, Eq. 21 for [P1_cyt] has a term –k_light[P1_cyt] added.

Fig. 7 A illustrates entrainment of the oscillations of Fig. 3 A to a circadian light-dark (LD) cycle. The LD cycle was divided into equal light and dark portions and the period was set to 22 h. During the light phase, k_light was set to 0.685 h^–1 for P0_cyt, P1_cyt, and P2_cyt. The light phase occurs during the falling portion of the PER time course, as it does experimentally (Lee et al., 1998). Oscillations entrained to a 24-h LD cycle look very similar to Fig. 7 A (not shown). The amplitudes of the oscillations in PER and the other species are not significantly greater than in simulated DD conditions (Fig. 3 A). The entrainment range appears reasonably broad. For parameters as in Fig. 7 A, entrainment occurs for an LD cycle length of 21–25 h.

A photic phase-response curve (PRC) was also constructed. Simulated light pulses were applied at 40 evenly spaced intervals during the circadian cycle in DD (Fig. 3, A and B). During the pulse duration, k_light was set to 0.685 h^–1 for P0_cyt, P1_cyt, and P2_cyt. The pulse duration chosen was 2 h. In previous models, durations of 1–3 h have been used for degradation of TIM in response to light pulses (Leloup and Goldbeter, 1998) or for degradation of an unspecified clock protein (Lema et al., 2000). Data illustrating PER disappearance subsequent to a light pulse lack sufficient temporal resolution to clearly determine the duration of increased PER degradation (Lee et al., 1996). These data also illustrate that under LD conditions, light pulses during the dark phase do not always accelerate PER disappearance. At Zeitgeber time 15, light pulses delay PER disappearance (Lee et al., 1996). Because the mechanism of this delay is not understood, and because PER regulation during LD cycles may differ significantly from that in DD, we have not modeled this additional complexity.

The phase advance or delay was determined seven cycles after each simulated light pulse (sufficient time for any transient dynamics to disappear). Fig. 7 B illustrates the resulting PRC (solid curve). Circadian time (CT) zero was chosen as coinciding with the peak of PER_tot during the unperturbed oscillation (Shafer et al., 2002). By the classification of PRCs given in Winfree (1987), this PRC is Type 1 (average slope of 0 when plotted as in Fig. 7 B). Fig. 7 B also illustrates an experimental PRC for Drosophila locomotor activity (data from Fig. 5 of Konopka et al., 1991). The agreement between the simulated and experimental PRCs appears quite good. The PRCs are similar in the magnitude of advances and delays, the number of hours of CT that correspond to advance versus delay, and the steepness of the crossover from delay to advance. The PRCs both have a dead zone of 0 phase shift at CT 5–9. The crossover from delay to advance occurs at CT 17 (simulated) and CT 18 (experimental). The experimental PRC of Matsumoto et al. (1994) is similar, with a dead zone at CT 5–9 and a steep crossover from delay to advance at CT 19.

Increases in the strength of the simulated light pulse (i.e., in k_light) continued to yield a Type 1 PRC. However, Saunders et al. (1994) illustrate Type 0 PRCs for Drosophila melanogaster. These PRCs were obtained with light pulses of very long duration (6 h). We therefore simulated very long-lasting increases in k_light and examined whether a Type 0 PRC could result. A duration of 8 h was used, along with values of k_light as high as 2.0 h^–1. The PRC remained Type I. In contrast, if light was assumed to degrade nuclear as well as cytosolic PER, Type 0 PRCs could be obtained with strong stimuli (e.g., k_light = 2.0 h^–1 for 3 h). With weaker stimuli, Type I PRCs were also obtained, but the shape of these PRCs failed to resemble experimental PRCs.

Temperature compensation of oscillation period can be simulated
A well-known characteristic of circadian oscillators is temperature compensation. At different ambient temperatures, the constant-darkness period of circadian oscillations remains virtually constant in Drosophila (Pittendrigh et al., 1973; Pittendrigh, 1954). This finding is a priori unexpected, because a temperature increase tends to increase the rate of each biochemical reaction. The temperature dependence of a reaction is commonly described by the Q₁₀ factor, which equals the ratio of the rate 10°C above a reference temperature to the rate at the reference temperature. Q₁₀ values of 2–3 are common for biochemical reactions (Segel, 1975). Since the period of a circadian oscillation is a function of the rates of many individual reactions, the period would be expected a priori to decrease as temperature increases.

One possible mechanism for temperature compensation is as follows. Processes such as nuclear transport of PER are composed of several elementary reactions (e.g., binding of protein to nuclear pore complexes, conformation changes of nuclear pore complexes). The rate of elementary reactions that hinder PER transport, such as dissociation of PER from nuclear pore complexes, might increase more rapidly with temperature than do the rates of reactions that favor PER transport. In this case a temperature increase could decrease the rate of PER nuclear transport, and a Q₁₀ < 1 would describe this decrease. Inhibiting PER nuclear transport would delay PER's repression of CLK-mediated transcription and lengthen the oscillation period (Hong and Tyson, 1997; Leloup and Goldbeter, 1997). This mechanism for temperature compensation is sometimes considered unsatisfying because it requires fine-tuning Q₁₀ values for separate kinetic processes. However, molecular evolution might have selected for phenotypes with a period resistant to moderate temperature changes, thereby resulting in fine-tuning of kinetic parameters.

Simulations suggest the model can qualitatively represent temperature compensation resulting from such fine-tuning. The control simulation of free-running oscillations (Fig. 3 A) was used as the starting point. Initially, a Q₁₀ of 2 was applied to every kinetic process. This corresponds to multiplying the scaling factor in Table 1 by the Q₁₀, doubling every rate constant and maximal velocity and halving each time consta