1. Introduction

Wide-type fruit flies, Drosophila melanogaster, might be the most extensively studied organism in circadian rhythm research. The researches of endogenous activity rhythm on Drosophila generally involve two different kinds of clock genes, called period (per, for short) [14,10] and timeless (tim, for short) [20,27]. Their encoded proteins, PER and TIM, bind to each other [5,0,27,29].

PER protein and per mRNA cycle in a $24$-hour period [7]. When PER protein is at a high level, per mRNA expression is repressed, suggesting that PER is an inhibitor of per mRNA accumulation [7]. The expression of per and tim genes is regulated by dCLOCK and CYC, and PER inhibits the transcription of per and tim by inactivating dCLOCK and CYC [1,3,18]. This negative feedback, introduced by PER inhibiting its own mRNA transcription, is the basis of many classic theoretical models of circadian rhythms [6,17,15,19].

An alternative way to study circadian rhythms is based on a positive feedback, introduced by PER phosphorylation being an activator to PER [26]. Phosphorylation of PER is operated by a double-time gene encoded kinase, DOUBLE-TIME (DBT, for short) [13,16]. As suggested by the dbt mutants phenotypes, PER phosphorylation might be precluded to its degradation. PER and TIM stimulate transcription of per and tim genes by activating dClOCK [2]. Experimental results suggest that per mRNA is stabilized by PER/TIM dimers [24], and PER is stabilized by dimerization with TIM [13,16].

The idea that PER phosphorylation introduces a positive feedback in PER accumulation can be expressed in a model of three-dimensional ordinary differential equations [26] (see (1) below). In [26], by imposing assumptions that the dimerization reactions were fast and dimeric proteins were in rapid equilibrium, they reduced the three-dimensional model to a pair of nonlinear ordinary differential equations of mRNA and total protein concentrations (see (2) below). Then they used the powerful phase plane portraiture to study the simplified two-dimensional model. In this paper, we explore the original three-dimensional model directly. It is shown that the circadian rhythms occur if the model possesses a unique equilibrium which is unstable. Furthermore, we deeply investigate how circadian rhythms are affected by several model parameters, including mRNA translation, mRNA degradation, monomer phosphorylation, protein proteolysis, association of PER/TIM protein and equilibrium constant for dimerization. The results help to explain some former-observed phenomena of circadian rhythms. In particular, our numerical results extremely agree with those given in [26], indicating that their reduction work is greatly reasonable.

---

2. Mechanism and mathematical model

In this section, we restate the model proposed by Tyson et al. [26]. The molecular mechanism for the circadian rhythm in Drosophila is summarized in Figure 1. Here the total PER (monomer $+$ dimer) degradation rate does not increase proportionally with the total PER concentration's increasing.

The mechanism in Figure 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 could not efficiently illustrate the importance of positive feedback in the reaction mechanism. So by noticing that PER and TIM messages and proteins followed roughly similar time courses in vivo, Tyson et al. [26] lumped them into a single pool of clock proteins. In addition, they assumed that the cytoplasmic and nuclear pools of dimeric protein were in rapid equilibrium. Then they established the following differential equations for [mRNA]=$M$, [monomer]=$P_1$, and [dimer]=$P_2$:

$\left\{ \begin{array}{lll} \displaystyle\frac{\mathrm{d}M}{\mathrm{d}t}&=&\displaystyle\frac{v_{m}}{1+(P_{2}/P_{crit})^2}-k_{m}M, \\ \displaystyle\frac{\mathrm{d}P_{1}}{\mathrm{d}t}&=&v_{p}M-\displaystyle \frac{k_{p_{1}}^{'}P_{1}}{J_{P}+P_{1}+rP_{2}}-k_{p_{3}}P_{1}-2k_{a}P_{1}^{2}+2k_{d}P_{2}, \\ \displaystyle\frac{\mathrm{d}P_{2}}{\mathrm{d}t}&=&k_{a}P_{1}^{2}-k_{d}P_{2}-\displaystyle \frac{k_{p_{2}}P_{2}}{J_{P}+P_{1}+rP_{2}}-k_{p_{3}}P_{2}.\\ \end{array} \right.$

(1)

Here monomer was assumed to be phosphorylated more quickly than dimer, i.e., $k_{p_{1}}^{'}\gg k_{p_{2}}$. The parameter $r$ determined the inhibition of dimer to monomer phosphorylation. For convenience, in this paper we follow them by taking $r=2$.

In their work, it was further assumed that the dimerization reactions were fast ($k_{a}$ and $k_{d}$ are large) such that monomers and dimers were always in equilibrium with each other. Then, by equilibrium conditions: $P_{2}=K_{eq}P_{1}^{2}$, $K_{eq}=k_{a}/k_{d}$, they obtained the reduced two-dimensional system:

$\left\{ \begin{array}{lll} \displaystyle\frac{\mathrm{d}M}{\mathrm{d}t}&=&\displaystyle\frac{v_{m}}{1+(P_{t}(1-q)/(2P_{crit}))^2}-k_{m}M, \\ \displaystyle\frac{\mathrm{d}P_{t}}{\mathrm{d}t}&=&v_{p}M-\displaystyle \frac{k_{p_{1}}P_{t}q+k_{p_{2}}P_{t}}{J_{P}+P_{t}}-k_{p_{3}}P_{t}, \\ \end{array} \right.$

(2)

where $P_{t}=P_{1}+2P_{2}=[\, \mbox{total protein}\, ]$, $k_{p_{1}} = k_{p_1}^{'} - k_{p_2} \approx k_{p_{1}}^{'}$, and

$q = q(P_{t}) = \displaystyle \frac{2}{1+\sqrt{1+8K_{eq}P_{t}}}.$

Two widely concerned points of circadian rhythms are whether the endogenous rhythms exist and how long the periods are. Since the mechanism has already been translated into mathematical models, attentions are drawn to examine the existence of periodic orbits and calculate the periods. In their work, system (2) has been thoroughly analyzed. In this paper, we try to study system (1). A typical oscillating solution of system (1) is illustrated in Figure 2, where the corresponding parameter values are chosen from Table 1.

---

3. Method and result

It is well-known that for higher dimensional ordinary differential equations, there is no so-called Poincaré-Bendixson theory: any limit set is a limit cycle if it contains no steady state. So, in order to use the powerful phase plane analysis tools, Tyson et al. [26] reduced (1) into (2) by imposing some assumptions. Fortunately, we observe that (1) is a three-dimensional competitive system in some sense [8,0,22,23]. For $n$-dimensional competitive ordinary differential equations, the dynamics is co-dimensional one. Every limit set lies on a Lipschitz manifold with one dimension lower, and this manifold is homeomorphic to an $(n-1)$-dimensional Euclidean space [28]. As for a three-dimensional competitive system, though there is no phase plane, one has a two-dimensional Lipschitz manifold (qualitatively exists but is unknown), where recurrent motions of system lie in. As a result, any limit set for such a system consists of either limit cycle, or steady state, or steady states connected with homoclinic or heteroclinic orbits. Suppose that all forward orbits for a three-dimensional competitive system are bounded and the system has a unique steady state. Then by the Perron-Frobenius theory, the linearized matrix at the steady state has a negative eigenvalue. The system has a one-dimensional stable manifold which is a strictly monotone curve [21], so it also rules out the third choice for limit set. For more details, please see [23] or the Appendix. We summarize the above discussion into the following theorem which can be found in [30]:

Theorem 3.1. Suppose (1) has a unique steady state $E$. If the linearized matrix of (1) at $E$ has one negative eigenvalue and two positive real part eigenvalues, then (1) has at least one stable limit cycle.

Numerical calculation suggests that (1) has a unique equilibrium in a large region of parameter values. However, limit cycles do not exist all the time. According to Theorem 3.1, when either $K_{eq}$ or $k_a$ is small, (1) has no limit cycles but a unique equilibrium (see Table 2), and the equilibrium appears to be a global stable steady state.

Comparing with the two-dimensional system (2), there are two more parameters $k_{a}$ and $k_{d}$ in system (1). Considering the equilibrium condition, $K_{eq}=k_{a}/k_{d}$, one only needs to detect how periods of (1) are influenced by $k_{a}$. As shown in Table 3 and Figure 3E, if we take $K_{eq}=200$, periodic orbits occur when $k_{a}$ is larger than a critical value $k_a^*=0.9$. In that region, as $k_a$ goes up, the period starts with a rapid decline, and then becomes quite insensitive. At first, we guess that the period is decreasing when $k_a$ is sufficiently large, but numerical calculations tell that it is not the case. In fact, the period even has a tendency to increase when $k_a$ is larger than $2.9\times10^6$ (see Table 3). The similar situations are observed with $K_{eq}=15$ (see Table 4).

Based on Tables 3 and 4, one can choose a suitable value of $k_a$ to calculate the rhythms for wild-type and mutant flies. The numerical results are presented in Table 5, where temperature compensation is found in wild-type flies but not in $per^{L}$ mutant flies.

Table 5 is due to the original three-dimensional system (1). As a comparison, we state Table 6, which is cited from [26] and based on the reduced two-dimensional system (2). Clearly, one can see that Table 5 and Table 6 are almost the same, which indicates that the reduction in [26] is greatly reasonable from this perspective.

In the next section, we will see more about the relation between circadian rhythms and parameters of (1).

---

4. Discussion

In the actual experiment, parameters of the circadian rhythms models are hard to be measured, or even unmeasurable. Parameter values in Table 1 have been chosen to yield a period close to 24-hours and ensure temperature compensation of the wild-type oscillator. The parameter values are arbitrary. Other combinations of parameter values may also yield circadian oscillations with possibly different periods.

It is significant to study how parameters of (1) affect its periodic oscillations. The numerical results are given in Figure 3, where the following parameters are considered: mRNA translation, mRNA degradation, monomer phosphorylation, protein proteolysis, association of PER/TIM protein and equilibrium constant for dimerization.

As shown in Figure 3A, periodic oscillation disappears when the protein synthesis rate $v_p$ is below a critical value. That coincides with the truth that the inhibiting effect of protein synthesis may eventually suppress the circadian rhythmicity [11,4,25]. Moreover, the period decreases when the protein synthesis rate is greater than a certain value. That matches the observations of anisomycin in the mollusk Bulla [12].

The PER/TIM complex formation plays a key role in the model. Circadian rhythm is markedly affected by the dimerization reaction, precisely in the model, by the association rate constant $k_a$ and dissociation rate constant $k_d$. In Figure 3E the period decreases as $k_a$ increases, which coincides with the suggestion that the heterodimeric dimerization is attenuated in the long-period $per^L$ mutant [5]. Here the attenuation is probably due to the competition of PER homodimeric complexes [0]. Note again that it always has $k_d=k_a/K_{eq}$ in this paper. The results of Figure 3F imply that circadian rhythm occurs only when $k_d$ is in a bounded range, and the period can be recognized as an increasing function of the dimer disassociation rate.

In Figure 3B we show how the oscillation is affected by mRNA synthesis. Periodic rhythm requires the mRNA synthesis rate $k_{m}$ to be bounded, implying that the oscillation may be destroyed if mRNA synthesizes either too slow or too fast. Moreover, the period of oscillation becomes shorter as the mRNA synthesis rate goes up. The qualitatively similar results (Figure 3D) are detected when we consider the proteins proteolysis rate $k_{p_3}$, except that periodic oscillation happens even if there is no proteins proteolysis.

According to Theorem 3.1, in Figure 4 we inspect the dependence of oscillations on parameters $K_{eq}$ and $k_{p_1}$. A U-shape region is found, whose boundary is almost the same as the locus of Hopf bifurcation in [26] (see Figure 4 in [26]). Within that region the system exhibits periodic oscillations, and in the outside area the system exhibits no limit cycle but a stable steady state. Figure 3C and 3F help to investigate the variation of period in this U-shape region. In Figure 3C, periodic oscillation requests that $k_{p_1}$ is neither too small nor too large, which means that protein monomers are sufficient but not too unstable. In Figure 3F, periodic oscillation vanishes when $K_{eq}$ is smaller than a critical value, which implies that the proteins tend to dimerize. Meanwhile, when $K_{eq}$ varies within a quite large region beyond a certain value, the period of (1) remains virtually unchanged, suggesting that the wild-type oscillation has temperature compensation (see also Table 5). Furthermore, when $K_{eq}$ decreases in a large region the period increases, which agrees with the consensus that $per^L$ mutant introduces a longer period for the $per^L$-encoded protein to reduce its tendency to form dimers [5,10]. By Table 5, one can also tell that $per^L$ mutant loses temperature compensation.

Appendix. The concentrations of mRNA, monomers and dimers are naturally nonnegative. We therefore focus on the first orthant $\mathrm{R}^3_+ =\{(M, P_1, P_2): M\geq 0, P_1\geq 0, P_2\geq 0\}$. It is easy to see that $\mathrm{R}^3_+$ is a positively invariant set of system (1), i.e., any solution $((M(t), P_1(t), P_2(t))$ of system (1) through a point in the first orthant lies in it when $t \geq 0$.

Let $a > v_{m} / k_{m}, \ \ b > v_{p}\, a / k_{p_{3}}, \ \ c = K_{eq}b^2$ and $B(a, b, c) = \{(M, P_1, P_2): 0 \leq M \leq a, \ 0 \leq P_1 \leq b, \ 0 \leq P_2 \leq c\}$. Denote $(f_1, f_2, f_3)$ the vector field of (1). By estimating the sign of the vector field at vertexes on the boundary of $B(a, b, c)$, one has

$\left\{ \begin{array}{lll} f_1(a, P_1, P_2)& = &\displaystyle\frac{v_{m}}{1+(P_{2}/P_{crit})^2}-k_{m}a < v_{m} -k_{m}a < 0, \\ f_2(M, b, P_2)& = &v_{p}M-\displaystyle\frac{k_{p_{1}}^{'}b}{J_{P}+b+rP_{2}}-k_{p_{3}}b-2k_{a}b^{2}+2k_{d}P_{2} \\ &<&v_{p}a-k_{p_{3}}b < 0, \\ f_3(M, P_1, c)& = &k_{a}P_{1}^{2}-k_{d}c-\displaystyle \frac{k_{p_{2}}c}{J_{P}+P_{1}+rc}-k_{p_{3}}c\\ &<& -\displaystyle \frac{k_{p_{2}}c}{J_{P}+P_{1}+rc}-k_{p_{3}}c < 0. \end{array} \right.$

The vector field for (1) on the boundary of $B(a, b, c)$ is shown in Figure 5, which implies that $B(a, b, c)$ is positively invariant. Note that for any point in $\mathrm{R}^3_+$, one can find such $(a, b, c)$ satisfying that box $B(a, b, c)$ contains the point. It follows immediately that any forward solution of system (1) is bounded. We summarize the above discussion into the following proposition:

Proposition 1. For any $a > v_{m} / k_{m}, \ \ b > v_{p}\, a / k_{p_{3}}$ and $c = K_{eq}b^2$, $B(a, b, c)$ is positively invariant for (1), that is, all forward solutions for (1) are bounded.

From Proposition 1, there are at least one steady state in $B(a, b, c)$. Suppose that the steady state $E$ is unique and there is no zero real part eigenvalue for the linearized matrix at $E$. Then the equilibrium $E$ is either locally asymptotically stable, or has a two-dimensional unstable manifold. The following arguments show that the latter case provides the existence of limit cycles.

By computing the Jacobian matrix of (1), one has

$Df=\left (\begin{array}{lll} -&0&- \\ +&-&+ \\ 0&+&- \end{array} \right),$

where ''$-$'' represents that the entry is strictly negative and ''$+$'' means strict positivity. According to [23], the system is competitive with respect to the cone $K=\{(M, P_1, P_2)\in\mathrm{R}^3: M\geq0, P_1\leq0, P_2\geq0\}$. By applying the theory on competitive systems in [23], we have

Theorem A. Suppose (1) has a unique steady state $E$. If the linearized matrix of (1) at $E$ has one negative eigenvalue and two positive real part eigenvalues, then (1) has at least one stable limit cycle.

Therefore, in order to study the oscillations for (1), one only needs to discuss its steady state and the local stability of the steady state.

---