MODELING TRANSCRIPTIONAL CO-REGULATION OF MAMMALIAN CIRCADIAN CLOCK

. The circadian clock is a self-sustaining oscillator that has a period of about 24 hours at the molecular level. The oscillator is a transcription-translation feedback loop system composed of several genes. In this paper, a scalar nonlinear diﬀerential equation with two delays, modeling the transcriptional co-regulation in mammalian circadian clock, is proposed and analyzed. Suﬃcient conditions are established for the asymptotic stability of the unique nontrivial positive equilibrium point of the model by studying an exponential polynomial characteristic equation with delay-dependent coeﬃcients. The existence of the Hopf bifurcations can be also obtained. Numerical simulations of the model with proper parameter values coincide with the theoretical result.

1. Introduction. Circadian clocks are endogenous 24-h oscillations that regulate the temporal organization of physiology, metabolism and behavior [4,11]. Disruption of circadian rhythms leads to various diseases and may reduce lifespan in mice [2,8,9,14]. The central mechanism of the mammalian circadian oscillation is a negative feedback loop that involves the transcriptional activator genes: clock and bmal1, and five repressor genes: period (per1-per3) and cryptochrome (cry1 and cry2) [11,20]. CLOCK and BMAL1 are transcription factors that activate per and cry gene transcriptions. The resulting PER and CRY proteins accumulate, and then inhibit CLOCK/BMAL1-mediated transcription after a certain time delay. Various auxiliary loops also take part in the regulation of mammalian circadian clocks. Previous studies validated that cry1 expression is positively auto-regulated via the inhibition of its repressor rev-erbα, which forms a positive loop [16]. Some nuclear receptors, exhibiting circadian like patterns of expression [18], may also contribute to the circadian clock via other auxiliary loops. These auxiliary loops are believed to receive the environmental factors for adaption purpose, therefore provide multiple entry points to regulate the circadian oscillator. Therefore, an important question arises: How does the system with more than one feedback loop maintain the robustness of the oscillation?
Goldbeter [5] proposed the first molecular model of circadian clock in Drosophila. This model takes account of the core negative feedback loop of per self-repression. Since then, a lot of mathematical models in different organisms based on the negative feedback loop have been developed and studied by many groups [6,10,15,19]. Recently, the circadian clock models with more than one feedback loops have been proposed, and the computational studies have revealed some interesting results [3,9,12,17]. However, the studies on the stability of the multi-loop regulation in circadian clock remain obscure.
Here we consider a system of one differential equation modeling the transcriptional co-regulation in a mammalian circadian clock. In this system, there is a core feedback loop that CRY1 negatively regulates its own expression by inhibiting CLOCK/BMAL1-mediated transcription. A delay naturally appears in this process, representing the durations of post-translational regulations. There occurs another auxiliary loop in this model in which CRY1 positively mediates its own expression via the inhibition of its repressor rev-erbα (see Fig. 1). Similarly, this positive process also contains a time delay. Therefore, we consider a two-loop system with two time delays in this work. The theoretical approach is based on the work of Adimy et al. [1]. Our aim is to show that the auxiliary loops, which are linked to other pathways for adaption purpose, will not abandon the oscillation with reasonable parameters. This paper is organized as follows: Section 2 presents the model of a mammalian circadian clock, a scalar nonlinear delay differential equation with two delays, and investigates the existence of a nontrivial positive equilibrium point. Section 3 analyses the asymptotic stability of this equilibrium point. We first linearize the model about the equilibrium point and obtain a first degree exponential polynomial characteristic equation. Then we determine the conditions for the stability when only one delay is equal to zero, and, eventually, when both delays are nonzero. Besides, we also establish the existence of the Hopf bifurcations, which destabilizes the system and leads to the existence of periodic solutions. Section 4 numerically illustrates the theoretical results and Section 5 discusses the effect of time delays on the period of the model. 2. The model. In mammalian circadian clocks, there is a core negative feedback loop to drive the oscillation. CRY protein combines with PER to form the dimer CRY/PER, which translocates into nucleus to inhibit CLOCK/BMAL1-mediated transcription of cry1, see Fig. 1(a). Besides the core loop, there are also some auxiliary loops in the mammalian circadian system to co-regulate clock. In this work, we focus on an important auxiliary loop via rev-erbα. The transcription of rev-erbα is also mediated by CLOCK/BMAL1 complex, so that can be inhibited by CRY/PER complex. In return, REV-ERBα feeds back to repress the transcription of cry1. The double-inhibition process forms a positive feedback loop (see Fig. 1(a)).
Our previous work about this system [17] illustrated some interesting numerical results. In order to obtain the theoretical insights of it, we first reduce the variables to form a simple one-dimensional system. We denote by x(t) the mRNA concentration of cry1 at time t. The transcription rate of cry1 (the synthesis rate of x(t)) is regulated by two items: one is the negative feedback through CLOCK/BMAL1 complex after a time delay τ 1 , and the other one is the positive feedback through REV-ERBα after a time delay τ 2 (see Fig. 1(b)). The degradation rate of cry1 mRNA is simply assumed to be proportional to itself. Therefore, the equation modeling the co-regulation of multiple loops in mammalian circadian clocks is: where k 1 and M are positive constants describing the negative regulation, k 2 is a positive constant characterizing the strength of positive regulation, and the constant c > 0 reflects the linear degradation rate. The parameter τ 1 > 0 denotes the average time needed by the negative feedback, which represents the durations of post-translational regulations. The parameter τ 2 > 0 denotes the average time needed by the positive feedback.
Obviously, system (2.1) has a unique continuous solution x(t) which is welldefined for all t ≥ 0 and for a continuous initial condition. Furthermore, it is easy to see that, for nonnegative initial conditions, the solutions of (2.1) remain nonnegative for t ≥ 0. In fact, if we let t 0 > 0 be such that x(t) > 0 for t < t 0 and x(t 0 ) = 0, then, from (2.1), we have and the result follows. Now, we begin to consider the existence of equilibrium points for system (2.1). An equilibrium point of system (2.1) is a stationary solution x * of (2.1), that is, Evidently, 0 is not an equilibrium point of system (2.1), so it is enough to search for its nonzero equilibrium points. For this purpose, we need to assume that c > k 2 and recall the well-known Shengjin's formulas.
Using Lemma 2.1, in Eq.(2.2), we calculate In addition, we can easily see that 2) has a positive solution, which is unique if and only if c > k 2 .
This result is summarized in the following proposition. 3. Local asymptotic stability. In this section, we concentrate on the stability of the nontrivial equilibrium point x * . Hence, we assume throughout this section that c > k 2 .
We first start to linearize system (2.1) around x * and deduce the characteristic equation.
3.1. Linearization and characteristic equation. Take y(t) = x(t) − x * , the linearized system of (2.1) about x * is then given by The characteristic equation associated with (3.1) is given by Through studying the sign of the real parts of roots of (3.3), we can analyse the local asymptotic stability of the equilibrium point x * . Here we recall that x * is locally asymptotically stable if and only if all roots of (3.3) have negative real parts, and its stability can only be lost if roots cross the vertical axis, that is, if purely imaginary roots appear. Because of the presence of two different delays, τ 1 and τ 2 , in Eq.(3.3), it is very complicated to analyse the sign of the real parts of eigenvalues, and there is no direct approaches to be considered. In the following, on the base of the analytic methods as in Ruan and Wei [13], we will discuss the stability of the equilibrium point when one delay is equal to zero, and deduce conditions for the stability of the equilibrium point when both time delays are nonzero.
(3.6) One can notice that if ω is a solution of (3.6), then so is −ω. Thus, in the following, we only look for positive solutions ω of (3.6).
Adding the squares of both hand sides of Eq.(3.6), we see that ω must be a root of the following equation: Set F (X) = X 2 − α 2 + (c − β) 2 and make the two assumptions as follows The function F has no positive zeros when (H1) holds. When (H2) holds, the function F has a unique positive zero Substituting ω 0 into Eq.(3.6), we have (3.8)

YANQIN WANG, XIN NI, JIE YAN AND LING YANG
In addition, by differentiating (3.5) with respect to τ 1 , we obtain From (3.9), we deduce that Then, Combining with (3.6), we get Based on the above analysis, we obtain the following Lemma 3.1.
According to Lemma 3.1 and the Hopf bifurcation theorem for delay differential equations, we conclude, the stability of the equilibrium point x * of (2.1) when τ 1 > 0, τ 2 = 0, in the following theorem. 3.4. The case τ 1 > 0, τ 2 > 0. We now return to investigate Eq.(3.3) with τ 1 > 0, τ 2 > 0. In order to study the local stability of the equilibrium point x * of (2.1), we regard τ 2 as a bifurcation parameter, take τ 1 = τ * 1 ∈ (0, ∞), and discuss the distribution of the roots of the following characteristic equation We first verified a result concerning the sign of the real parts of characteristic roots of (3.10) with τ * 1 ∈ (0, τ 0 1 ) in the following lemma. Lemma 3.2. If all roots of Eq. (3.5) have negative real parts for τ * 1 ∈ (0, τ 0 1 ), then there exists a τ * 2 (τ * 1 ) > 0 such that all roots of Eq.(3.10) have negative real parts when τ 2 < τ * 2 (τ * 1 ). Proof. From Theorem 3.1 (i), we know that Eq.(3.5) has no root with nonnegative real part for τ * 1 ∈ (0, τ 0 1 ). Obviously, the left hand side of Eq.(3.10) is analytic in λ and τ 2 . Following Theorem 2.1 of Ruan and Wei [13], as the parameter τ 2 varies, the sum of the multiplicity of zeros of the left hand side of Eq. (3.3) in the open right half-plane can change only if a zero appears on or crosses the imaginary axis.
4. Numerical simulations. In this section, we illustrate the different stability results obtained in the previous sections, mainly in Theorems 3.1 and 3.2. We also focus on periodic solutions appearing through a Hopf bifurcation. Without loss of generality, we take time unit as an hour, and let initial condition be x 0 = 0.1.
From Fig. 5, we can see that system (4.1) turns its instability into stability, or turns its stability into instability on both sides of the bifurcation line. Especially, when we fix τ 1 = 2.4 and take τ 2 = 0.5, 2, 5, 9, 12, 15.5 which correspond to the six different points in Fig. 5, we simulate the stability of system 4.1 ( see Fig. 6).
Furthermore, in order to verify the theoretical result of the bifurcation Fig. 5, we use numerical simulations and obtain all the values of τ 1 , τ 2 which make system (4.1) generate oscillations (see Fig. 7).
From Fig. 7, we can see that the equilibrium point x * of system (4.1) is unstable when (τ 1 , τ 2 ) locates in the black region, that is to say, system (4.1) generates oscillating solutions. The oscillating range of (τ 1 , τ 2 ) is in accordance with the result of theoretical calculations in Fig. 5.

5.
Discussion. In the previous researched circadian clock models with a time delay, it is found that the period of the model monotonously increases with the increase of the delay. In this paper, there are two different delays in our discussing biological clock model. To analyse the characteristic equation with two delays, we first concentrated on the case when one of the delays, τ 2 equals zero and obtained a critical value for the delay τ 1 : when τ 1 < τ 0 1 all roots of the characteristic equation have negative parts and when τ 1 = τ 0 1 purely imaginary roots appear. Then we assumed that τ 1 = τ * 1 < τ 0 1 and considered the delay τ 2 as a parameter. It showed that there exists a τ * 2 (τ * 1 ) such that all roots of the characteristic equation have negative real parts when τ 2 ∈ (0, τ * 2 (τ * 1 )). Finally, we assumed that τ 1 = τ * 1 ∈ (τ 0 1 , τ 0 1 + ǫ) with ǫ > 0 and considered the delay τ 2 as a parameter. We concluded that there exists a τ * * 2 (τ * 1 ) such that the characteristic equation has at least one root with strictly positive real parts when τ 2 ∈ (0, τ * * 2 (τ * 1 )). Consequently, we obtained stability and instability results for the mammalian circadian clock model with two independent delays respectively. The numerical results in Fig. 5 also indicate the existence of a Hopf bifurcation that leads to the emergence of periodic solutions.
When discussing a mammalian circadian clock model with two delays, we found delays can affect the period of the model. When we first fix τ 2 = 29, assuming that τ 1 ranges from 0 to 12, then fix τ 1 = 10, assuming that τ 2 ranges from 0 to 35, we numerically simulated the influence of delays on the period of model (4.1)(see Fig.  8).  Figure 7. Oscillating range of (τ 1 , τ 2 ) for system (4.1). Black regions represent oscillating solutions with periods for system (4.1) when (τ 1 , τ 2 ) locates in the black region.  Figure 8. The effect of time delays on the period of system (4.1).
In figure (a), we fix τ 2 = 29, the black solid line represents the relation between τ 1 and the period. In figure (b), we fix τ 1 = 10, the black solid line represents τ 2 and the period.
From Fig. 8(a), we can see the period is monotonously increasing with the increase of time delay τ 1 when τ 2 = 29, which is in accord with that of a biological clock model with a single delay. However, from Fig. 8(b), we find that the period emerges non-monotonic regular variety with the increase of time delay τ 2 when τ 1 = 10.
All in all, it is seen from Fig. 8 that τ 1 determines the length of time. When τ 1 is fixed, τ 2 is changed, the period can only change in a certain range, which shows that the positive feedback between CRY1 and REV-ERBα in mammalian circadian clocks can benefit the robustness of the period length.