A computational model clarifies the roles of positive and negative feedback loops in the Drosophila circadian clock

Physics Letters A 374 (2010) 2743–2749

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

A computational model clarifies the roles of positive and negative feedback loops in the Drosophila circadian clock Junwei Wang a,∗ , Tianshou Zhou b a b

Cisco School of Informatics, Guangdong University of Foreign Studies, Guangzhou 510006, China School of Mathematics and Computational Science, Sun Yat-Sen University, Guangzhou 510275, China

a r t i c l e

i n f o

Article history: Received 7 April 2009 Received in revised form 18 April 2010 Accepted 21 April 2010 Available online 24 April 2010 Communicated by C.R. Doering Keywords: Drosophila Circadian clock Transcriptional feedback loop Negative feedback loop Positive feedback loop

a b s t r a c t Previous studies showed that a single negative feedback structure should be sufficient for robust circadian oscillations. It is thus pertinent to ask why current cellular clock models almost universally have interlocked negative feedback loop (NFL) and positive feedback loop (PFL). Here, we propose a molecular model that reflects the essential features of the Drosophila circadian clock to clarify the different roles of negative and positive feedback loops. In agreement with experimental observations, the model can simulate circadian oscillations in constant darkness, entrainment by light–dark cycles, as well as phenotypes of per01 and clkJrk mutants. Moreover, sustained oscillations persist when the PFL is removed, implying the crucial role of NFL for rhythm generation. Through parameter sensitivity analysis, it is revealed that incorporation of PFL increases the robustness of the system to regulatory processes in PFL itself. Such reduced models can aid understanding of the design principles of circadian clocks in Drosophila and other organisms with complex transcriptional feedback structures. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The circadian clock is a daily time-keeping mechanism responsible for the ∼ 24-hour biochemical rhythm of almost all organisms ranging from unicellular cyanobacteria through filamentous fungi and plants to mammals. Its main function is to anticipate daily environmental changes to organize the physiology and behavior of organisms such that they occur at biologically advantageous times. In recent decades, many of the molecular aspects comprising the circadian clock have been elucidated for model organisms as diverse as Drosophila, Neurospora and mammals [1], thanks to advances in molecular biology experiments. At the molecular level, circadian oscillations are based on the rhythmic expression of clock genes. In Drosophila, a number of circadian clock genes have been identified and cloned since Konopka and Benzer [2] provided the first evidence that a single gene can have an impact on circadian rhythms. These genes include period (per), timeless (tim), clock (clk), cycle (cyc), doubletime (dbt), cryptochrome (cry), and vrille (vri) [3]. As is now commonly accepted, the Drosophila intracellular circadian oscillator is composed of interacting negative feedback loop (NFL) and positive feedback loop (PFL) that drive recurrent rhythms in the mRNA and protein levels of key clock components [4,5]. The NFL involves rhythmic

*

Corresponding author. Tel.: +86 20 39328577. E-mail address: [email protected] (J.W. Wang).

0375-9601/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2010.04.059

inhibition of the CLK–CYC driven by negative regulators. Specifically, CLK and CYC first form a heterodimer and then activate the rhythmic transcription of per and tim genes by binding to E-box sequences (CACGTG) in their promoters. The resultant PER and TIM proteins accumulate during the night and form a complex that enters the nucleus and binds to transcription factors CLK and CYC. The complex between PER (and perhaps TIM) and the CLK/CYC prevents CLK/CYC from binding to E-box sequences in the per and tim promoters and results in a cessation of transcriptional activity [6]. This process creates a NFL where PER and TIM inhibit their own transcription. In a PFL, clk mRNA oscillations are caused by the rhythmic inhibition of clk transcription through VRI protein. If clk expression is slightly activated, then resulting activation of per and tim transcription by CLK/CYC results in binding of CLK/CYC by PER/TIM. Thus, repression of clk expression is relieved, and the level of CLK increases further, indicating that clk transcription is essentially positively regulated by PER and TIM. Mathematical modeling has emerged as a powerful tool for gaining understanding of the dynamic of complex gene regulatory networks. The models can be used to identify the structure of the physiologic system, integrate experimental data, make predictions about properties of the system that are not apparent, conduct formal statistical analyses of parameters and, if well constructed, guide the design and direction of future investigations [7–9]. A number of models have been created on the basis of experimental observations to better characterize the molecular network underlying circadian rhythmicity in Drosophila [10–20].

2744

J.W. Wang, T.S. Zhou / Physics Letters A 374 (2010) 2743–2749

The models where PER represses its own gene expression [10] or where the complex PER/TIM inhibits transcription of per and tim genes [11,12] have demonstrated that a single negative feedback loop could produce circadian oscillations of mRNAs and proteins. At this stage, it is natural to ask why current description of Drosophila circadian clocks suggests that they have interlocked NFL and PFL [4,21]. Several authors have begun to speculate on the reasons for this complexity. For example, it has been previously suggested that the PFL would increase robustness either to stochastic fluctuations [22] or to gene mutations [16]. However, it has not been shown that the robustness results from the observed structure and there is no convincing explanation of why one could expect this [14,23,24]. These discussions suggest that the function of NFL and PFL in Drosophila rhythm generation is not well understood. The goal of the present Letter is to construct a molecular model for the Drosophila circadian clock to further clarify the different roles of NFL and PFL. The model with a reduced set of variables reflects the essential of the Drosophila circadian oscillator. The model can simulate circadian oscillations in constant darkness (DD) and entrainment by light–dark (LD) cycles with agreement to experimental observations. We also use the model to simulate phenotypes of per01 and clkJrk mutants. Sustained oscillations persist when the PFL is removed, implying the crucial role of NFL for rhythm generation. Through parameter sensitivity analysis, we found that elimination of PFL only decreases the robustness of oscillations to regulatory processes in PFL itself, whereas the robustness of oscillations toward regulatory processes in NFL remains stable with and without PFL. Such reduced models can aid understanding of the design principles underlying circadian clocks in Drosophila and other organisms with complex transcriptional feedback structures. 2. Description of the mathematical model

Fig. 1. Schematic diagram representing the model for the Drosophila circadian clock involving genes per and clk. CLK protein activates transcription of Per gene but inhibits transcription of clk gene. In the negative transcriptional feedback loop, binding of the PER to CLK prevents the expression of Per gene, whereas in the positive transcriptional feedback loop, PER indirectly enhances clk expression by binding to CLK to form the complex PER–CLK. Variable names used in the model are indicated in the parentheses.

PER (PERP ) and CLK (CLKP ) is assumed as was done by Leloup and Goldbeter [28]. We assume that these phosphorylated proteins are subject to degradation by enzymes. Such degradation is also considered for mRNAs and the PER–CLK complex. The dynamics of the model, schematized in Fig. 1, is described by the following system of 9 differential equations:

dx1 dt

The model of the Drosophila circadian oscillator is schematized in Fig. 1, whilst the variables representing the concentration of mRNAs and proteins are indicated. This model incorporates the essential features of molecular processes of Drosophila and describes the dynamics of both NFL and PFL. The following assumptions and rationale of the assumptions are invoked in developing the mathematical model given in Fig. 1:

dx2

(1) At the mRNA level, it is already clear that cyc expression in Drosophila is constitutive and produces constant level of CYC [4]. Our model only considers the action by CLK without explicitly modeling the interaction of CYC with CLK. (2) The negative elements PER and TIM are essential for maintenance of the Drosophila circadian rhythms. The expression of per and tim genes is coregulated by CLK/CYC complex. Also, the phase of nuclear accumulation PER and TIM proteins is similar [25]. Based on these considerations, we represent the dynamics of per and tim genes, i.e., per and tim mRNAs, by a single “lumped” variable, “per mRNA”, which can be thought of as an average level of per and tim mRNAs; similarly, PER and TIM proteins are represented by a single variable “PER”. This representation also reduces the number of ordinary differential equations describing the dynamics of our model. (3) Except for activating per and tim transcription, the CLK/CYC heterodimer also activates expression of gene vri. The VRI protein then represses clk transcription [21]. In our model, we assume that clk expression is directly, negatively autoregulated by CLK without taking dynamic of vri gene as explicit model variables. (4) Although the number of protein phosphorylations that occur may be considerable [26,27], a single phosphorylated state for

dt dx5

dt dx3 dt dx4

dt dx6 dt

= v per

xm 8 km per

+

xm 8

= v PER x1 − T 1P

− v mp

x1 kmp + x1

x2

= T 1P

L 1P + x2

− T 2P

(1)

x3

+ T 2P

L 1P + x2 − k1P x2 + k2P x4 − d2 x2 x2

− d1 x1 L 2P + x3

(2)

x3

− v dPC

L 2P + x3

x3 L d + x3

− d3 x3

= k1P x2 − k2P x4 − k1PC x4 x8 + k2PC x9 − d4 x4 = v clk

knclk knclk + xn8

= v CLK x5 − T 1C

− v mc

x5 kmc + x5

x6 L 1C + x6

(4)

− d5 x5

+ T 2C

(5)

x7 L 2C + x7

− k1C x6 + k2C x8 − d6 x6 dx7 dt dx8 dt dx9 dt

= T 1C

x6 L 1C + x6

− T 2C

(6)

x7 L 2C + x7

− v dCC

x7 L d + x7

= k1C x6 − k2C x8 − k1PC x4 x8 + k2PC x9 − d8 x8 = k1PC x4 x8 − k2PC x9 − v dn

(3)

x9 L d + x9

− d9 x9

− d7 x7

(7) (8) (9)

The variables of the system represent the following species: x1 , per mRNA; x2 , PER protein in the cytosol; x3 , phosphorylated PER protein in the cytosol (PERP ); x4 , PER protein in the nucleus; x5 , clk mRNA; x6 , CLK protein in the cytosol; x7 , phosphorylated CLK protein in the cytosol (CLKP ); x8 , CLK protein in the nucleus; x9 , PER– CLK complex in the nucleus. The biochemical meaning of the corresponding parameters is explained in Table 1. Concentrations of clock components are in units of nM.

J.W. Wang, T.S. Zhou / Physics Letters A 374 (2010) 2743–2749

2745

Table 1 Model parameters. Parameter

Value

Biochemical description

v per m kper v mp kmp v PER T 1P T 2P L 1P L 2P k1P k2P v dPC Ld k1PC k2PC v clk n kclk v mc kmc v CLK T 1C T 2C L 1C L 2C k1C k2C v dCC v dn dj

1.5 nM h−1 4 0.7 nM 1.1 nM h−1 0.31 nM 0.37 h−1 1.4 nM h−1 1.3 nM h−1 0.1 nM 0.2 nM 0.4 h−1 0.2 h−1 0.7 nM h−1 0.3 nM 0.5 nM−1 h−1 0.1 h−1 0.9 nM h−1 3 2.2 nM 0.9 nM h−1 0.4 nM 0.16 h−1 1.5 nM h−1 1.1 nM h−1 0.1 nM 0.1 nM 0.4 h−1 0.2 h−1 0.5 nM h−1 0.8 nM h−1 0.01 h−1

Maximum rate of per transcription Hill coefficient of activation of per transcription Activation constant for enhancement of per expression by CLK Maximum rate of degradation of per mRNA Michaelis constant for degradation of per mRNA Rate constant for synthesis of PER Maximum rate for phosphorylation of PER in cytosol Maximum rate for dephosphorylation of PER in cytosol Michaelis constant for phosphorylation of PER in cytosol Michaelis constant for dephosphorylation of PER in cytosol Rate constant for entry of cytoplasmic PER into nucleus Rate constant for transport of nuclear PER into cytosol Maximum rate of degradation of phosphorylated PER Michaelis constant for degradation by enzymes Rate constant for association of PER and CLK Rate constant for dissociation of PER–CLK complex Maximum rate of clk transcription Hill coefficient of inhibition of clk transcription Inhibition constant for repression of clk expression by CLK Maximum rate of degradation of clk mRNA Michaelis constant for degradation of clk mRNA Rate constant for synthesis of CLK Maximum rate for phosphorylation of CLK in cytosol Maximum rate for dephosphorylation of CLK in cytosol Michaelis constant for phosphorylation of CLK in cytosol Michaelis constant for dephosphorylation of CLK in cytosol Rate constant for entry of cytoplasmic CLK into nucleus Rate constant for transport of nuclear CLK into cytosol Maximum rate of degradation of phosphorylated CLK Maximum rate of degradation of PER–CLK complex Nonspecific degradation rate constant for the jth quantity

Because parameter values of the model cannot be reliably known from experimental data, we choose their appropriate values based on previous models [29–31] with modifications so as to satisfy that (1) the model must generate circadian oscillations with a period close to 24 h, and (2) can correct the phase relationship between mRNAs and proteins. With these considerations, a set of standard parameter values for simulation are listed in Table 1. The system of ordinary differential equations was solved numerically by using a 4th order Runge–Kutta algorithm with a fixed step size of 10−2 . All simulations were implemented using Matlab (The MathWorks, Natick, MA, USA). 3. Results 3.1. The model can simulate Drosophila circadian oscillations under DD For a model of the circadian clock, we should first test whether its dynamics behavior is consistent with experimental observations. With parameter values in Table 1, our model can reproduce experimental observed sustained oscillations with a circadian period 23.5 h (see Fig. 2). Because the effects of light were not simulated, these oscillations are endogenous and correspond to circadian rhythms that persist in DD. The time course of per and clk mRNAs and proteins are plotted in Fig. 2. It can be seen that the oscillations of per and clk mRNAs are quite symmetric in appearance, with a gradual rise and fall (Fig. 2(a)). The oscillation of per mRNA peaks at the circadian time (CT)13 given that per mRNA level bottoms at CT0. The cytoplasmic PER peaks at CT15. On the other hand, clk mRNA peaks at CT3.5 and subsequently bottoms at CT12.6. These results are consistent with observations that per mRNAs reach peak levels early in the evening (CT12-CT16) [32], while the clk mRNA levels peak at late night to early in the morning (CT23-CT4) and oscillates in anti-phase to per mRNA levels [26,33]. The peak of clk mRNA precedes the peak of CLK protein by about 5.5 h. In contrast, the lag

Fig. 2. Sustained circadian oscillations for the concentrations of mRNAs and proteins with parameter values in Table 1. (a) Time courses of per and clk mRNAs. (b) Time courses of PER and CLK proteins. The system of ordinary differential equations was solved numerically by using a 4th order Runge–Kutta algorithm with a fixed step size of 10−2 . All simulations were implemented using Matlab (The MathWorks, Natick, MA, USA).

between the peaks of per mRNA and PER protein is shorter and about 2 h. These results are qualitatively consistent with experiments that there exist delays of several hours between per (clk) mRNA and PER (CLK) [33]. 3.2. Simulation of per01 and clkJrk mutants Starting from standard parameter values (Table 1), we will simulate the phenotypes of gene mutants per01 and clkJrk by changing

2746

J.W. Wang, T.S. Zhou / Physics Letters A 374 (2010) 2743–2749

Fig. 3. Simulation of phenotypes of the per01 and clkJrk mutants. (a) Oscillations of clock components are abolished in the per01 single-mutant flies by reducing v PER from 0.37 h−1 to 0 h−1 . (b) The clkJrk mutant results in complete loss of circadian oscillations with v CLK = 0 h−1 . The remaining parameter values are the same as in Table 1 except for v PER and v CLK .

particular parameters according to the functionality of the mutants. In Drosophila, per01 and clkJrk refer to null mutations in per and clk genes which produce nonfunctional proteins. Therefore, to simulate per01 and clkJrk , the translation rates of proteins PER and CLK were set to 0. Fig. 3(a) and (b) show simulations of per01 and clkJrk . It can be observed that sustained circadian oscillations have been abolished in both per01 and clkJrk mutant flies, roughly consistent with observation that rhythmicity of per and clk mRNAs is blocked in per01 and clkJrk [33]. For per01 , the clk mRNA level is very low, qualitatively consistent with the observation of a low clk mRNA level in per01 flies [4]. The simulated results show high levels of per mRNA (Fig. 3(a)). However, experiments carried in Drosophila have shown that in mutants lacking PER (per01 ), per mRNA level are constitutive and low [34]. These apparent contradiction results can be explained from the structure of our model: per01 -induced absence of PER leads to free CLK from the PER–CLK complex, which, in turn, causes a very high level of CLK. As the activation effects of CLK, the level of per mRNA is higher than its peak in wild type (WT). Consequently, strong repression of clk transcription from CLK itself makes a low level of clk mRNA. Fig. 3(b) shows the phenotype of clkJrk mutant where circadian oscillations in per and clk mRNAs disappeared. These simulated data agree with experimental data that clkJrk flies exhibit low levels of per mRNA [21] and a high level of clk mRNA which is near the WT peak [4]. We can explain the mechanism underlying these resulting data from the model structure as follows. In clkJrk mutant, lack of functional CLK leads to the absence of activation effects of per gene. Consequently, complete loss of levels of per mRNA follows. A high amount of clk mRNA is present because the repression effect from CLK itself is inhibited and the transcription of clk gene is constitutive with constant rate v clk = 0.9 nM h−1 . 3.3. The role of positive feedback loop in circadian rhythm generation It has been suggested that both NFL and PFL are necessary for sustaining circadian oscillations [22,35]. Hastings [22] suggested that an oscillator based on a single NFL would progressively dampen. To test whether the oscillations disappear without PFL re-

Fig. 4. When the PFL is eliminated, sustained oscillations of some components are preserved. (a) and (b) correspond to time course of mRNAs and proteins, respectively, where kclk = 500 nM. Other parameter values are the same as those in Table 1.

lying on repression of clk expression by CLK, we fix kclk = 500 nM. At such a high value of the inhibition constant, the effect of inhibition of CLK on clk transcription is negligible. Such a model without PFL differs from earlier models of the Drosophila circadian oscillator that eliminate positive feedback [14,15] in the sense that those models assumed that the nuclear CLK protein was fixed to a constant. The simulations in Fig. 4 illustrate that oscillations in per mRNA persist although the levels of clk mRNA were kept at a constant high value. The period of these oscillations is close to 28.8 h and is markedly different from 24 h. The anti-phase relationship between PER and CLK proteins was also maintained as in the model with PFL. In the absence of PFL, the oscillations originate solely from the PER negative autoregulation. The above simulations have established that the PFL is not required to sustain circadian oscillations per se. What is the role of the PFL in circadian rhythm generation? The parameter sensitivity analysis [24] was used to further investigate the role of PFL in Drosophila circadian clock. Because oscillation period is one of the quantities of primary interest in oscillatory system, we employ period sensitivity which yields a quantitative measure of the deviations in period resulting from perturbation of system parameters. The relative (or scaled) periodic sensitivity [24,36] was calculated according to Eq. (10):

S (τ ; p j ) =

p j ∂τ ∂ ln τ = · ∂ ln p j τ ∂pj

(10)

where p j is the parameter with parameter index j, τ is the period of the system. A lower sensitivity of a parameter implies a higher robustness of the period to the parameter of a model. As the number of parameters is very large, it is difficult to thoroughly investigate the behavior of the model in such a high dimensional parameter space. To assess the sensitivity of our model to changes in parameter values, we only change one parameter at a time while holding other parameters at their standard values. It should be noted that this period sensitivity is only locally valid with respect to parameter space, that is, in a neighborhood of a specific parameter set. Changes of ±5% in each parameter were made and the resulting oscillations patterns were examined, as shown in Fig. 5. Oscillations were found to persist in all cases. With PFL, the period tends to be more robust to parameter varia-

J.W. Wang, T.S. Zhou / Physics Letters A 374 (2010) 2743–2749

2747

Fig. 5. Period sensitivity with respect to system parameters in the Drosophila circadian model. (a) and (b) correspond to +5% perturbation and −5% perturbation to each parameter, respectively. The remaining parameter values are the same as those in Table 1.

tions for components in NFL (e.g., v per for transcription rate of per gene, v PER for translation rate of PER) (Fig. 5). Without PFL, period sensitivity of the model with respect to parameter variations for components in NFL was not very different from those with PFL. On the other hand, with the elimination of PFL, the period sensitivity of the model with respect to parameter variations for components in PFL was significantly increased, where differences among parameter classes are much clearer for the interacting negative and positive feedback model than for the simpler negative feedback model [10,11]. This analysis indicates that the regulation of per and tim transcription constitutes a robust part of the Drosophila circadian clockwork. The additional branch of PFL incorporated in the more complex Drosophila circadian model increases the robustness of the system to regulatory processes in PFL itself. In modeling gene regulatory network, regulatory processes (e.g., transcriptional control or phosphorylation of proteins) can be described by parameters in the mathematical models. Hence our analysis of period sensitivity also provides clues on the importance of individual regulatory processes on the function of the Drosophila circadian clock. Both models (i.e., with and without PFL) show high sensitivity toward parameter perturbations affecting clk expression and CLK phosphorylation degradation, such as v clk and v mc (Fig. 5). These regulatory mechanisms are of similar importance for the robustness of such two clock architectures (i.e., with and without PFL). The intuitive finding that the oscillation period is highly sensitive to protein phosphorylation degradation (e.g., v mc ) complies with experimental evidence linking human disorders of the sleep– wake cycle (e.g., familial advanced sleep phase syndrome (FASPS)) to the phosphorylation status of PER protein [37,38].

Fig. 6. Effects of light on the Drosophila circadian clock. (a) Entrainment by a 24 h LD cycle. klight varies periodically from the high value of 1.5 nM h−1 during the 12 h long light phase (white rectangle) to 0 nM h−1 during the 12 h long dark phase (black rectangle). (b) Influence of the magnitude of the change in the lightdependent parameter on the phase of per mRNA for LD cycles in which klight passes from 0 nM h−1 to 1.5, 3, 5 and 7 nM h−1 , respectively. (c) Rhythmicity disappears under constant light conditions when klight is set to the constant high value of 1.5 nM h−1 . Other parameter values are the same as those in Table 1.

In Drosophila, the proteins PER and TIM physically interact, and the timing of their association and nuclear localization drive cycles of per and tim transcription through an autoregulatory feedback loop. It has been experimentally shown that light enhances degradation of phosphorylated TIM [27,39] and increase of TIM degradation rate has been used to simulate entrainment to LD cycles in Drosophila circadian models [11,12,18]. When TIM is released from the PER/TIM complex, PER phosphorylation is strongly enhanced [40]. It is plausible that accelerated degradation of highly phosphorylated PER would result, and experiments have demonstrated light-induced enhancement of PER phosphorylation and degradation [41]. Because our Drosophila model does not have a separate variable for TIM, we will simulate the effect of light by adding a first-order degradation term to the right-hand side of each differential equation for cytosolic PER, i.e.,

dx2 dt

A circadian rhythm that persists under constant conditions is called a free-running rhythm which represents one of the defining characteristic of circadian rhythms. A second defining characteristic of circadian rhythms is that they can be reset (a process termed entrainment) by Zeitgebers (e.g., the daily light and temperature cycles) with light as the predominant Zeitgeber. Entrainment results from perception of light by one or more clock components and shifts the circadian clock to an appropriate phase.

x2 L 1P + x2

+ T 2P

x3 L 2P + x3

(2 )

− k1P x2 + k2P x4 − d2 x2 − klight x2 dx3 dt

3.4. Response of the Drosophila circadian clock to light

= v PER x1 − T 1P

= T 1P

x2 L 1P + x2

− T 2P

− d3 x3 − klight x3

x3 L 2P + x3

− v dPC

x3 L d + x3

(3 )

where klight represents the degradation rate induced by light. In conditions of light–dark (LD) cycles, the parameter klight varies periodically (as a square cycle), going from 0 during the dark phase up to a higher constant value kmax during the light phase. Fig. 6(a) shows the entrainment dynamics of circadian oscillations to a 12 : 12 LD cycle (12 h of light followed by 12 h of darkness). It can be seen that the switch from dark to light occurs during the falling portion of the PER time course, as it does ex-

2748

J.W. Wang, T.S. Zhou / Physics Letters A 374 (2010) 2743–2749

perimentally [42]. The shape of the time course of PER in Fig. 6(a) is slightly different from the time course in DD (Fig. 2(b)). Experiments carried out in Drosophila [43] have shown that the peak of per mRNA appears after ∼ 4 h the onset of the dark phase, regardless of the relative duration of the light and dark phases in such LD cycles of 24-h period. Simulations of our model account for such observations, as shown in Fig. 6(b) where four values of kmax are considered for the light-induced enhancement of PER degradation. This result is, to a large extent, independent of the intensity of the increase in PER degradation triggered by light and indicate that the beginning of the dark phase corresponds to the rise portion of the time course of per mRNA. Our model also shows that the circadian rhythmicity disappears under constant light (LL) conditions when klight is kept at a constant high value (klight = 1.5) (Fig. 6(c)). This agrees with experimental data that although LL does not immediately stop the clock, it does affect the activity or stability of some clock molecules, which leads to a gradual dampening of the rhythm to arrhythmicity [43]. 4. Conclusion and discussion We have presented a molecular model for the Drosophila circadian clock to clarify the different roles of negative and positive feedback loops. Although a detailed model is useful for examining elementary molecular processes of the Drosophila circadian clock [15,18], a reduced model as presented here is more appropriate for the investigation of the role of PFL. In accordance with experimental observations, our model with lower number of variables and parameters successfully predicts: (1) sustained circadian oscillations in DD (Fig. 2); (2) phenotypes of per01 and clkJrk mutants (Fig. 3); (3) response to external light clues (Fig. 6). In addition, our simulations in Fig. 4 also illustrate that circadian oscillations of per mRNA, PER and CLK persist if PFL is eliminated by neglecting the auto-repression of CLK on clk expression. This implies that the NFL in which PER represses per expression by interacting with CLK is essential for Drosophila circadian oscillations. Based on this result, we propose an experimentally testable prediction as follows. Based on clk-null mutant animals, a transgenic Drosophila line might be constructed in which clk expression in lateral neurons is driven by a constitutive promoter [44,45]. We predict that circadian oscillations of PER and CLK will stay rhythmic in these neurons. Current description of the molecular mechanisms of circadian clocks among different organism suggests that they almost universally have a complicated structure with interlocking NFL and PFL [46]. In two previous models [14,15], the function of PFL in the Drosophila clock has been discussed. Their simulations suggest that PFL is not required to sustain circadian oscillations. Through parameter sensitivity analysis, however,we found that PFL increases the robustness of the model to regulatory processes in PFL itself. Specially, although with and without PFL the robustness of oscillations toward regulatory processes in NFL remains stable, the robustness of oscillations to regulatory processes in PFL itself was significantly reduced in the absence of PFL. Another possibility of the function of the PFL is that it may regulate expression of clockcontrolled genes (CCGs) outside the core oscillator. In Drosophila, microarray experiments have identified more than 100 CCGs, majority of which are regulated by protein CLK [47]. Therefore, the clk-mediated PFL is likely to be essential for behavioral aspects of Drosophila circadian oscillations. Theoretically speaking, a real biological system cannot be accurately described by a mathematical model. In the simulation processes, our model fails to predict protein cycling when per is expressed under control of constitutively active promoter [48]. Tests carried in Drosophila have shown that in mutants lacking PER (per01 ), per mRNA level are constitutive and low [34]; however,

our simulated results show high levels of per mRNA. These deficiencies may be because of some clearly oversimplifications of the known mechanism of the Drosophila circadian timekeeping mechanism in developing the theoretical model, e.g. per and tim genes are lumped into single variables, and only a single phosphorylation state for PER and CLK was considered. The successful simulations of our model, however, suggest that such a reduced model with essential set of variables, despite its simplicity, capture important features of the processes underlying the Drosophila circadian clock. Intertwined NFL and PFL also characterize circadian rhythms in other organisms such as the cyanobacterium Synechococcus, the filamentous fungus Neurospora and mammals. Therefore, reduced models similar to ours might be useful for gaining intuitive understanding of the role of PFL in circadian rhythm generation in these organisms. Under natural conditions, the Drosophila clockgene-expressing neurons reside in discrete sites of the brain and consist of multiple autonomous single-cell oscillators, and the experimental evidence strongly suggests that robustness in timekeeping precision emerges at the cellular population level and not at the single-cell level [49]. The mathematical model presented here is a good tool to approach questions of coupling and synchronization of circadian neurons in Drosophila brain [50,51] and this will be our future work. Acknowledgements This work was supported by the Natural Science Key Foundation of China (Grant No. 60736028), the Natural Science Foundation of China (Grant No. 10871074, 60704045) and Research Fund of Youth Scholars for the Doctoral Program of Higher Education of China (Grant No. 20070558053). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

J.C. Dunlap, Cell 96 (1999) 271. R.J. Konopka, S. Benzer, Proc. Natl. Acad. Sci. USA 68 (1971) 2112. J.A. Williams, A. Sehgal, Annu. Rev. Physiol. 63 (2001) 729. M.R. Glossop, L.C. Lyons, P.E. Hardin, Science 286 (1999) 766. P.E. Hardin, Curr. Biol. 15 (2005) R714. C. Lee, K. Bae, I. Edery, Mol. Cell. Biol. 19 (1999) 5316. E.N. Brown, H. Luithardt, J. Biol. Rhythms 14 (1999) 609. D.B. Forger, C.S. Peskin, Proc. Natl. Acad. Sci. USA 100 (2003) 14806. E.B. Klerman, M.S. Hilaire, J. Biol. Rhythms 22 (2007) 91. A. Goldbeter, Proc. R. Soc. London B Biol. Sci. 261 (1995) 319. J.C. Leloup, A. Goldbeter, J. Biol. Rhythms 13 (1998) 70. J.J. Tyson, C.I. Hong, C.D. Thron, B. Novak, Biophys. J. 77 (1999) 2411. P. Smolen, D.A. Baxter, J.H. Byrne, J. Neurosci. 21 (2001) 6644. P. Smolen, D.A. Baxter, J.H. Byrne, Biophys. J. 83 (2002) 2349. P. Smolen, P.E. Hardin, B.S. Lo, D.A. Baxter, J.H. Byrne, Biophys. J. 86 (2004) 2786. H.R. Ueda, M. Higiwara, H. Kitano, J. Theor. Biol. 210 (2001) 401. P. Ruoff, M.K. Christensen, V.K. Sharma, J. Theor. Biol. 237 (2005) 41. Z. Xie, D. Kulasiri, J. Theor. Biol. 245 (2007) 290. T.L. Leise, E.E. Moin, J. Theor. Biol. 248 (2007) 48. N. Bagheri, M.J. Lawson, J. Stelling, F.J. Doyle III, J. Biol. Rhythms 23 (2008) 525. S.A. Cyran, A.M. Buchsbaum, K.L. Reddy, M.C. Lin, N.R. Glossop, P.E. Hardin, M.W. Young, R.V. Storti, J. Blau, Cell 112 (2003) 329. M.H. Hastings, Nat. Rev. Neurosci. 1 (2000) 143. D.A. Rand, B.V. Shulgin, J.D. Salazar, A.J. Millar, J. Theor. Biol. 238 (2006) 616. J. Stelling, E.D. Gilles, F.J. Doyle III, Proc. Natl. Acad. Sci. USA 101 (2004) 13210. K. Bae, C. Lee, P.E. Hardin, I. Edery, J. Neurosci. 20 (2000) 1746. I. Edery, L.J. Zwiebel, M.E. Dembinska, M. Rosbash, Proc. Natl. Acad. Sci. USA 91 (1994) 2260. H. Zeng, Z. Qian, M.P. Myers, M. Rosbash, Nature 380 (1996) 129. J.C. Leloup, A. Goldbeter, Proc. Natl. Acad. Sci. USA 100 (2003) 7051. T.L. Leise, E.E. Moin, J. Theor. Biol. 248 (2007) 48. J.C. Leloup, A. Goldbeter, Proc. Natl. Acad. Sci. USA 100 (2003) 7051. S. Becker-Weimann, J. Wolf, H. Herzel, A. Kramer, Biophys. J. 87 (2004) 3023. P.E. Hardin, J.C. Hall, M. Rosbash, Nature 343 (1990) 536. K. Bae, C. Lee, D. Sidote, K.Y. Chuang, I. Edery, Mol. Cell. Biol. 18 (1998) 6142. W.V. So, M. Rosbash, EMBO J. 16 (1997) 7146. S.K. Crosthwaite, J.C. Dunlap, J.J. Loros, Science 276 (1997) 763.

J.W. Wang, T.S. Zhou / Physics Letters A 374 (2010) 2743–2749

[36] D.E. Zak, J. Stelling, F.J. Doyle III, Comput. Chem. Eng. 29 (2005) 663. [37] K.L. Toh, C.R. Jones, Y. He, E.J. Eide, W.A. Hinz, D.M. Virshup, L.J. Ptacek, Y.H. Fu, Science 291 (2001) 1040. [38] T. Ebisawa, M. Uchiyama, N. Kajimura, K. Mishima, Y. Kamei, M. Katoh, T. Watanabe, M. Sekimoto, K. Shibui, K. Kim, Y. Kudo, Y. Ozeki, M. Sugishita, R. Toyoshima, Y. Inoue, N. Yamada, T. Nagase, N. Ozaki, O. Ohara, N. Ishida, M. Okawa, K. Takahashi, T. Yamauchi, EMBO Rep. 2 (2001) 342. [39] M.P. Myers, K. Wager-Smith, A. Rothenfluh-Hilfiker, M.W. Young, Science 271 (1996) 1736. [40] B. Kloss, A. Rothenfluh, M.W. Young, L. Saez, Neuron 30 (2001) 699. [41] C. Lee, V. Parikh, T. Itsukaichi, K. Bae, I. Edery, Science 271 (1996) 1740. [42] C. Lee, K. Bae, I. Edery, Neuron 21 (1998) 857. [43] J. Qiu, P.E. Hardin, Mol. Cell. Biol. 16 (1996) 4182.

2749

[44] Z. Yang, A. Sehgal, Neuron 29 (2001) 453–467. [45] M. Gallego, E.J. Eide, M.F. Woolf, D.M. Virshup, D.B. Forger, Proc. Natl. Acad. Sci. USA 103 (2006) 10618. [46] D. Bell-Pedersen, V.M. Cassone, D.J. Earnest, S.S. Golden, P.E. Hardin, T.L. Thomas, M.J. Zoran, Nat. Rev. Genet. 6 (2005) 544. [47] M.J. McDonald, M. Rosbash, Cell 107 (2001) 567. [48] Z. Yang, A. Sehgal, Neuron 29 (2001) 453. [49] E.D. Herzog, S.J. Aton, R. Numano, Y. Sakaki, H. Tei, J. Biol. Rhythms 19 (2004) 35. [50] J.W. Wang, J.J. Zhang, Z.J. Yuan, A.M. Chen, T.S. Zhou, J. Biol. Rhythms 23 (2008) 472. [51] J.W. Wang, A.M. Chen, J.J. Zhang, Z.J. Yuan, T.S. Zhou, Chin. Phys. B 18 (2009) 1294.