INTRODUCTION

TOP ABSTRACT INTRODUCTION MECHANISM AND MATHEMATICAL... PHASE PLANE PORTRAITS COMPARISON TO EXPERIMENTS DISCUSSION REFERENCES

Wild-type fruit flies, Drosophila melanogaster, exhibit endogenous activity rhythms with a period of 24 h over a broad temperature range (18-33°C). The first mutation to interfere with this circadian rhythm was discovered by Konopka and Benzer (1971), who called the gene period (per, for short). Three mutant alleles of per have been studied: per^L and per^S, with endogenous activity rhythms of 27 and 19 h, respectively (at 18°C), and per⁰, a null allele with no overt rhythm (Huang et al., 1995). Remarkably, the per^L mutant has lost temperature compensation; the period of its endogenous rhythm increases from 25 h at 15°C to 33 h at 30°C (Huang et al., 1995).

A second important gene, timeless or tim, encodes a protein, TIM, that binds to PER (Gekakis et al., 1995; Myers et al., 1995; Sehgal et al., 1994, 1995; Vosshall et al., 1994; Zeng et al., 1996). Mutation of tim abolishes the circadian rhythm (Sehgal et al., 1994). During endogenous cycling in constant darkness, a brief light pulse causes a phase shift of the circadian rhythm (Myers et al., 1996; Pittendrigh, 1967). This phase shift has recently been attributed to rapid degradation of TIM upon exposure to light (Hunter-Ensor et al., 1996; Lee et al., 1996; Myers et al., 1996; Zeng et al., 1996).

PER protein and per mRNA fluctuate with a 24-h period, with protein lagging behind mRNA by 4-6 h (Hardin et al., 1990; Zeng et al., 1994). When PER protein is overexpressed from a constitutive promoter, expression of endogenous per mRNA is repressed (Zeng et al., 1994), suggesting that PER inhibits its own transcription (Hardin et al., 1990). Binding to TIM seems to be necessary for translocation of PER to the nucleus (Vosshall et al., 1994) to exert its inhibitory effect. PER forms homo- and heterodimers through its "PAS" domain (Gekakis et al., 1995; Huang et al., 1995; Lee et al., 1996; Zeng et al., 1996), which it shares with many transcription factors but not with TIM. The per^L mutation, which lies in the PAS domain, disrupts PER/PER (Huang et al., 1995) and PER/TIM binding (Gekakis et al., 1995). Expression of the per and tim genes is regulated by a pair of transcription factors, dCLOCK (also called JRK) and CYC, that appear to be inactivated by PER (Allada et al., 1998; Darlington et al., 1998; Rutila et al., 1998). This evidence for negative feedback of PER on transcription of its own mRNA is the basis for most current theoretical models of circadian rhythms (Goldbeter, 1995; Ruoff and Rensing, 1996; Leloup and Goldbeter, 1998; Scheper et al., 1999). However, we propose that a positive feedback loop, based on stabilization of PER by dimerization with TIM, may play an equally important role in generating oscillations. This proposal is supported by recent discoveries on PER phosphorylation and proteolysis.

PER is phosphorylated by a casein-like kinase called DBT (encoded by the double-time gene), which is present at roughly constant levels during the rhythm (Kloss et al., 1998; Price et al., 1998). PER phosphorylation seems to be a prelude to its degradation, as suggested by the phenotypes of dbt mutants. In dbt^P, which codes for a nonfunctional kinase and has no rhythm, PER accumulates in a hypophosphorylated form. dbt^S codes for a more active kinase, accumulates less PER than wild type, and has shorter cycles (18 h in homozygote). dbt^L codes for a less active kinase, accumulates more PER than wild type, and has longer cycles (26.8 h in homozygote). Experimental results suggest that PER is stabilized on association with TIM (Kloss et al., 1998; Price et al., 1998).

Other avenues of positive feedback are also possible in Drosophila. For instance, Suri et al. (1999) present evidence that PER/TIM dimers stabilize per mRNA, and the experiments of Bae et al. (1998) suggest that PER and TIM are transcriptional activators of dCLOCK, which in turn stimulates transcription of the per and tim genes.

	MECHANISM AND MATHEMATICAL MODEL

TOP ABSTRACT INTRODUCTION MECHANISM AND MATHEMATICAL... PHASE PLANE PORTRAITS COMPARISON TO EXPERIMENTS DISCUSSION REFERENCES

Following the lead of Kloss et al. (1998), we assume that PER monomers are rapidly phosphorylated and degraded, whereas PER/TIM dimers are less susceptible to proteolysis (poorer substrates for either DBT or the proteolytic machinery). Our model, summarized in Fig. 1, is similar in structure to that of Leloup and Goldbeter (1998), but with a crucial difference. In the Leloup-Goldbeter model, the role of PER phosphorylation is to introduce a time delay into the negative feedback loop. In our model, the role of PER phosphorylation is to introduce positive feedback in PER accumulation. As PER concentration increases, an ever greater proportion of protein is dimerized and protected from DBT. Therefore, as the total concentration of PER (monomer + dimer) increases, the rate of total PER degradation does not increase proportionally. This nonlinearity is a key factor in the following mathematical model of circadian rhythms.

View larger version (21K): [in this window] [in a new window]   FIGURE 1   A simple molecular mechanism for the circadian clock in Drosophila, adapted from Price et al. (1998), Kloss et al. (1998), Wilsbacher and Takahashi (1998), and Leloup and Goldbeter (1998). PER and TIM proteins (rectangle and oval, respectively) are synthesized in the cytoplasm, where they may be destroyed by proteolysis or they may combine to form relatively stable heterodimers. Heteromeric complexes are transported into the nucleus, where they inhibit transcription of per and tim mRNA. We assume that PER monomers are rapidly phosphorylated by DBT and then degraded. Dimers, we assume, are poorer substrates for DBT.

The mechanism in Fig. 1 could be translated into a set of six differential equations, for per and tim mRNAs, PER and TIM monomers, and PER/TIM dimers in the cytoplasm and nucleus. Such a complicated set of equations would not effectively illustrate the importance of positive feedback in the reaction mechanism. Noticing that PER and TIM messages and proteins follow roughly similar time courses in vivo, we lump them together into a single pool of clock proteins. In addition, we assume that the cytoplasmic and nuclear pools of dimeric protein are in rapid equilibrium. With these simplifying assumptions, our model reduces to three differential equations for [mRNA] = M, [monomer] = P₁, and [dimer] = P₂:

(1)

(2)

(3)

In Eq. 1, we assume some cooperativity (represented by a Hill coefficient of 2) for the inhibition of mRNA transcription by PER/TIM dimers. This assumption is not essential: the equations exhibit limit cycle oscillations even if the Hill coefficient = 1. However, we find it easier to fit some properties of mutant fruit flies with a Hill coefficient of 2. Of more importance is our assumption, in Eqs. 2 and 3, that both monomers and dimers bind to DBT, but P₁ is phosphorylated more rapidly (k_p1' k_p2). It is essential for oscillations that the DBT-catalyzed reaction shows saturable kinetics (e.g., Michaelis-Menten) and that the dimer is a competitive inhibitor of monomer phosphorylation. The extent of competitive inhibition is determined by r, the ratio of enzyme-substrate dissociation constants for the monomer and dimer. Oscillations are observed for r as small as 0.2, but not for r = 0. These properties of the DBT-catalyzed reaction have not yet been determined experimentally, so they constitute testable predictions of our theory. Finally, the terms involving k_p3 in Eqs. 2 and 3, which represent slow, first-order degradation of monomers and dimers, are not essential for oscillations, but they allow a better fit to some of the data.

If the dimerization reactions are fast (k_a and k_d large), then P₁ and P₂ are always in equilibrium with each other. Let P_t = P₁ + 2P₂ = [total protein]. Then, from the equilibrium condition, P₂ = K_eqP₁², K_eq = k_a/k_d, we can write P₁ = qP_t and P₂ = 1/2(1  q)P_t, where

(4)

Now our model reduces to a pair of nonlinear ordinary differential equations:

(5)

(6)

supplemented by the algebraic equation (Eq. 4), which specifies that q = q(P_t). In Eq. 6, k_p1 = k_p1'  k_p2  k_p1'. Also, in Eq. 6, we have assumed that r = 2. If r  2, the denominator of the Michaelis-Menten expression should be written as J_p + qP_t + (r/2)(1  q)P_t. A typical solution of Eqs. 4-6 is illustrated in Fig. 2.

View larger version (21K): [in this window] [in a new window]   FIGURE 2   Numerical solution of model of Eqs. 5 and 6, given the parameter values in Table 1.

	PHASE PLANE PORTRAITS

TOP ABSTRACT INTRODUCTION MECHANISM AND MATHEMATICAL... PHASE PLANE PORTRAITS COMPARISON TO EXPERIMENTS DISCUSSION REFERENCES

To reach the pair of differential equations (Eqs. 5 and 6), we have made a number of assumptions about interactions between PER and TIM. In a later publication, we will relax these assumptions and study the full set of kinetic equations implied by Fig. 1. But for our present purpose of emphasizing the role played by positive feedback in PER dynamics, the two-equation model has a great advantage in being amenable to the powerful tools of phase plane analysis (Edelstein-Keshet, 1988). The phase plane for Eqs. 5 and 6 is the Cartesian coordinate system representing our two variables, mRNA and total protein. In the phase plane (Fig. 3 A), we plot two nullclines: 1) where mRNA synthesis is exactly balanced by mRNA degradation,

(M-nullcline)

and 2) where protein synthesis and degradation are balanced,

(P-nullcline)

The M-nullcline is sigmoidally shaped because we assume cooperativity of the inhibition of transcription by PER/TIM dimers. The P-nullcline is N-shaped because 1) dimer is more stable than monomer and 2) dimer competitively inhibits binding of monomer to DBT.

View larger version (18K): [in this window] [in a new window]   FIGURE 3   Phase plane portraits. (A) M- and P-nullclines (see text) for the parameter values in Table 1. The closed orbit is a stable limit cycle oscillation. (B) Shift of the P nullcline as Keq is changed from 200 to 1. When Keq = 1, the system exhibits a stable steady state with abundant mRNA but a low protein level.

A sketch of the nullclines in the phase plane is called a "portrait" of the dynamical system, and it tells us much about the system's temporal behavior. For instance, the portrait in Fig. 3 A shows the system oscillating around a limit cycle in the phase plane. We obtain this portrait by adjusting the parameters of the model so that the M-nullcline intersects the P-nullcline on its intermediate branch. Then the rate constants are scaled to give an oscillation of (nearly) 24 h (see Table 1). If we were to change some of the parameters (e.g., by mutation), the nullclines would move across the phase plane, and the portrait would change (Fig. 3 B). In particular, the period of oscillation may change or the limit cycle may disappear altogether and be replaced by a stable steady state.

In Fig. 4 we show how the temporal behavior of the model depends on K_eq and k_p1. Within the U-shaped region, the model exhibits stable limit cycle oscillations. Clearly, oscillatory behavior requires that K_eq be large enough (i.e., protein subunits tend to dimerize) and k_p1 be large enough (i.e., protein monomers are sufficiently unstable). For K_eq > 100, the period of oscillation is close to 24 h and is quite insensitive to changes in K_eq or k_p1, suggesting that the rhythm of wild-type flies may be insensitive to temperature-induced changes in parameters. For K_eq < 50, the period of oscillation abruptly increases and becomes quite sensitive to both K_eq and k_p1. It is known that the defect in per^L-encoded protein reduces its tendency to form dimers (Gekakis et al., 1995; Huang et al., 1995), which leads naturally to longer periods of per^L mutants and the temperature sensitivity of their rhythm (Fig. 4). Of course, variations of the other parameters with temperature also affect the periods of wild-type and mutant flies. By choosing appropriately (Table 1) the activation energies for each parameter in the model (Ruoff et al., 1997; Ruoff, 1992), we readily account for temperature compensation of the wild-type rhythm and temperature dependence of the oscillator in per^L flies (Table 2).

View larger version (22K): [in this window] [in a new window]   FIGURE 4   Two-parameter bifurcation diagram for the model, calculated with the software tool AUTO (Doedel, 1986). As Keq and kp1 are allowed to vary, with all other parameter values fixed as in Table 1, we find a U-shaped region of two-parameter space where limit cycle oscillations are possible, bounded by a locus (H-H) of Hopf bifurcations. We foliate this region with curves of constant period, from 24 to 40 h. Notice that the period of oscillation varies little from 24 h in a large region of parameter space, but then becomes quite sensitive to parameter values when Keq is small.

	COMPARISON TO EXPERIMENTS

TOP ABSTRACT INTRODUCTION MECHANISM AND MATHEMATICAL... PHASE PLANE PORTRAITS COMPARISON TO EXPERIMENTS DISCUSSION REFERENCES

Our simple model can be tested against the phenotypes of several other circadian rhythm mutants. We account for the properties of dbt^S and dbt^L mutants (Table 2) by assuming that k_p1 is not much affected by these mutations (otherwise the rhythm would likely be destroyed; see Fig. 4), but k_p2 is increased or decreased 10-fold. On the other hand, period is not much affected by multiple doses of wild-type dbt (Table 2). The phenotypes of per⁰ (null mutant), per^OP (constitutive promoter), and dbt^p (null mutant) are easily explained (not shown). Regarding the dosage dependence of wild-type per, the model predicts a slight increase in period with increasing v_m (the period at v_m = 1 is 0.6 h longer than the period at v_m = 0.5), but genetic manipulations (Smith and Konopka, 1982; Cote and Brody, 1986) show a slight decrease in period (the period of per⁺/per⁺ is 0.5 h shorter than the period of per⁺/per⁰). In light of the simplicity of the model, this discrepancy does not seem too serious.

Any model of circadian rhythms should also be tested for its response to light pulses and its propensity to be entrained by light-dark cycles. To simulate typical phase-response curves (PRCs), we assume that the effect of light is to reduce K_eq. (The literature reports that light pulses increase the degradation of TIM (Hunter-Ensor et al., 1996; Lee et al., 1996; Myers et al., 1996; Zeng et al., 1996), suggesting that we should model a light pulse by increasing k_p3; however, we found that this assumption does not produce correctly shaped PRCs in our model, whereas decreasing K_eq does. Because monomeric protein is unstable in our model, if light absorption were to destabilize dimers, then the clock protein would rapidly be lost.) The pattern of delays and advances (Fig. 5 A) is qualitatively similar to PRCs observed experimentally (Myers et al., 1996; Pittendrigh, 1967), although close comparison will show many quantitative discrepancies. The delay or advance shows up in the first cycle after treatment, after which the oscillator is back on its 24-h track.

View larger version (15K): [in this window] [in a new window]   FIGURE 5   Interaction of the endogenous oscillator with light. (A) Phase response curve. At a chosen phase of the endogenous rhythm we perturb the equations by setting Keq = 15 for a definite period of time (10 or 30 min) and then return Keq to 200. We follow the perturbed rhythm for 100 h and compare the time when Pt = 1 to the time of this event in the unperturbed rhythm. We plot the phase difference (advance or delay) as a function of the phase at the onset of perturbation. Phase = 0 (Pt = 0.023, M = 0.82) corresponds roughly to subjective dusk (lights out). (B) Principal entrainment band. Equation 4 is modified so that Keq = 200 for one-half of a Zeitgeber period (lights off), and Keq = 200(1  a) for the other half (lights on). Presumably, a is proportional to the intensity of illumination. As a function of Zeitgeber period (abscissa), we look for the minimum value of a (ordinate) required for 1:1 entrainment of the circadian oscillator to the Zeitgeber.

In addition, the model is very rapidly entrained to an external Zeitgeber with a period different from 24 h (Fig. 5 B).

	DISCUSSION

TOP ABSTRACT INTRODUCTION MECHANISM AND MATHEMATICAL... PHASE PLANE PORTRAITS COMPARISON TO EXPERIMENTS DISCUSSION REFERENCES

The earliest attempts to model circadian rhythms were not based on molecular mechanisms, because nothing of the sort was known. Rather, they emphasized the generic properties of limit cycle solutions to nonlinear dynamical systems (Kalmus and Wigglesworth, 1960; Kronauer et al., 1982; Pavlidis, 1967; Tyson et al., 1976; Wever, 1960; Winfree, 1970). To their credit, they revealed the sort of behavior that can be expected of the circadian clock regardless of the actual make-up of its springs, gears, and levers. As soon as it became clear that repression of gene transcription plays an important role in circadian timekeeping, models based on Goodwin's (1965) negative-feedback paradigm appeared (Goldbeter, 1995; Ruoff and Rensing, 1996). Goodwin's equations, first used to model periodic enzyme synthesis in bacteria, describe a "pure" negative feedback loop (no autocatalytic terms):

View larger version (6K): [in this window] [in a new window]

(7)

Although the equations are simple and elegant, their analysis is complex and unintuitive. The conditions for sustained oscillations can be severe: n  3 (i.e., phase-plane portraiture is impossible), p > 8 for n = 3 (i.e., extreme levels of cooperativity are often required), and k_i  k_j for all i and j (i.e., all components of the loop must be comparably unstable) (Griffith, 1968; Thron, 1991; Tyson and Othmer, 1978). Furthermore, in the more than 30 years since the mechanism was first proposed, only one biochemical example of a pure negative-feedback oscillator has been carefully described (Bliss et al., 1982), even though end-product repression is a common theme of biochemical pathways.

Instead of concentrating on the negative feedback loop implied by PER's inhibition of per transcription, we propose that the crucial molecular interaction generating oscillatory behavior of the PER network is the positive feedback loop implied by the stabilization of PER protein upon dimer formation. To emphasize the role that positive feedback may play in the circadian system, we have drastically simplified the molecular machinery (lumping together PER and TIM) so that the model can be described by two differential equations. Our simple two-component model has many advantages over recent models based on pure negative feedback: phase-plane analysis now gives useful insight into the mechanism of oscillation, highly cooperative transitions are no longer required, and reasonable rate constants can be assigned to the elementary steps. Our model exhibits a remarkable insensitivity of the oscillatory period to certain crucial parameters (which preadapts the mechanism for temperature-compensated rhythms), and it is consistent with the phenotypes of many distinctive mutants at the per, tim, and dbt loci. We can account for certain qualitative features of phase resetting and synchronization in response to light by assuming that K_eq is light sensitive, but not by making the more realistic assumption that light exposure increases k_p3. These strengths and weaknesses of the model suggest that positive feedback plays a heretofore unrecognized role in the dynamics of circadian rhythms, but that more comprehensive molecular mechanisms and mathematical models will be required for accurate representation of the detailed properties of circadian rhythms in Drosophila.

The change in emphasis, from negative feedback to positive feedback, provides a clue to the evolution of the endogenous circadian clock. Consider a primitive mechanism that lacks negative feedback on transcription. In this case, dM/dt = v_m  k_mM replaces Eq. 5, and the phase-plane portrait of the system changes dramatically (Fig. 6 A). An organism with this primitive mechanism would not exhibit endogenous oscillations, but it could still be entrained by external light/dark cycles. With positive feedback in effect, the N-shaped P-nullcline creates a hysteresis loop (Fig. 6 B) that can be traversed by changes in the equilibrium constant for dimerization. When the lights are on, K_eq is small, and the protein level is small because it is mostly monomeric and rapidly degraded. When the lights are off, K_eq is large, and the protein level is large because it is mostly dimeric and more stable. Such an organism would still know the time of day, as long as it was subjected to external light/dark cycles. By adding negative feedback later, the process of natural selection could convert a "switch" into a "clock." An unexpected benefit is that the period of the endogenous oscillation is quite insensitive to the mechanism's crucial parameter values (Fig. 4), so it has a built-in preadaptation for temperature compensation.

View larger version (17K): [in this window] [in a new window]   FIGURE 6   A circadian "switch" based on positive feedback alone. (A) The P-nullcline is plotted for Keq = 15 (day) and Keq = 200 (night). The M-nullcline is a horizontal line (M = vm/km) because P does not inhibit synthesis of M in this model. During the day, the system sits on a stable steady state with Pt  0.008 (few dimers, rapid degradation of monomers). At night the system sits on a stable steady state with Pt  10 (mostly dimers, slow degradation). (B) Hysteresis loop. We plot the steady-state level of Pt as a function of Keq. Assuming, as before, that Keq is inversely proportional to the intensity of illumination, we predict that expression of PER protein will show distinct thresholds: turning on at low intensity (Keq = 155), but not turning off until a much higher intensity is reached (Keq = 21.6).

If by mutating the transcription factors, clock and cyc, geneticists can isolate flies that synthesize PER and TIM constitutively (i.e., no negative feedback); these mutant flies, although they have no endogenous rhythm, would still respond perfectly well to light, synthesizing PER at night and degrading it during the day. Any model in which TIM is rapidly degraded by exposure to light would predict this effect, because TIM is necessary for PER accumulation. What sets our model apart is the prediction that the switch between synthesis and degradation should show hysteresis as a function of intensity of illumination (Fig. 6 B). Furthermore, by knocking out the phosphorylation sites on PER, the hysteresis loop should be eliminated.

Our model of the circadian clock in Drosophila is surely incomplete, because the molecular basis of circadian rhythms is still vigorously studied and hotly debated. We think of it as a skeletal or minimal model to be elaborated and improved. We cannot expect such a simple model to explain correctly all features of circadian rhythms, but it does rest firmly on the current knowledge of the molecules involved, and it gives new insight into the central roles played by proteolysis, dimerization, and competitive inhibition in generating bistability and oscillations. In addition, our representation of the central control system by two differential equations means that the intuitively pleasing, geometric ideas of phase plane analysis can now be applied to circadian rhythms.