1. Introduction

An important goal of system biology is to understand the mechanisms that govern the protein regulatory dynamics. The so-called feedback loop is believed to be widely present in various eukaryotic cellular processes and regulates the gene expression of a broad class of intracellular proteins. It is usually understood that in the feedback loop transcription factors are regulatory proteins which activate transcription in eukaryotic cells; they act by binding to a specific DNA sequences in the nucleus, either activating or inhibiting the binding of RNA polymerase to DNA; mRNA is transcribed in the nucleus and is in turn translated in the cytoplasm. Such a feedback loop was described by systems of ordinary differential equations in [6, 7, 8,17,21].

Many scholars have noticed that the genetic regulatory dynamics are dependent on the intracellular transport of mRNA and protein [13,18]. In particular, models of genetic regulatory dynamics are developed to explain a variety of sustained periodic biological phenomena [19]. Hence the possible roles played by the time delays in the signaling pathways become an interesting research topic and we are interested how time delays are associated with periodic oscillations. Regulatory models with constant delays are investigated in the work of [18,24], among many others, where time delay is used as a parameter and is tuned to observe stable periodic oscillations. It was pointed out in [19] that if the simplification of constant time delay is too drastic, another approach is to assume a distributed time delay. In the work of Busenburg and Mahaffy [2], both diffusive macromolecular transport represented as a distributed delay and active macromolecular transport modeled as time delay are considered. It was demonstrated that stable periodic oscillations can occur when transport was slowed by increasing the time delay, and that oscillations were dampened if transport was slowed by reducing the diffusion coefficients.

We follow the idea of the work of [2] to consider both diffusion transport and active transport in genetic regulatory dynamics, while we start from the Goodwin's model [6] in contrast to the reaction-diffusion equations in [2]. Namely, we will provide a new perspective to investigate genetic regulatory dynamics which generalizes the previous extension of Goodwin's model with constant delay and simplifies the approach adopted by [2]. To be more specific, we begin with the prototype model for Hes1 system (see Section 2 for a brief introduction) and discuss how to model the diffusion time if we consider the effect of the fluctuation of concentrations of the substances in the Hes1 system. We show that modeling the diffusion time will lead to a model with threshold type distributed delay, which is also called threshold type state-dependent delay [1,16]. With a general model with threshold type state-dependent delay, we are interested how the model with state-dependent delay differs from those with only constant delays, how the stabilities of the steady state and the periodic oscillations depend on the state-dependent delay.

We organize the remaining of the paper as follows. In Section 2, we model the diffusion time in regulatory dynamics and set up a prototype system of differential equations with state-dependent delay for the Hes1 system. In Section 3, we consider a general model with state-dependent delay and conduct linear stability analysis of the equivalent model and conduct a local Hopf bifurcation analysis of the system using the multiple time scale method. Numerical simulations on the model of Hes1 system with state-dependent delay will be illustrated in Section 4. We conclude the investigation at Section 5.

---

2. State-dependent diffusion time in regulatory dynamics

In this section, we motivate the modeling of diffusion time delay with the Hes1 regulatory system. Hes1 is a protein which belongs to the basic helix-loop-helix (bHLH, a protein structural motif) family of transcription factors. It is a transcriptional repressor of genes that require a bHLH protein for their transcription. As a member of the bHLH family, it is a transcriptional repressor that influences cell proliferation and differentiation in embryogenesis. Hes1 regulates its own expression via a negative feedback loop, and oscillates with approximately 2-hour periodicity [9]. However, the molecular mechanism of such oscillation remains to be determined [9]. Mathematical modeling of Hes1 regulatory system has been an active research area in the last decade. See among many others, [12,20].

The basic reaction kinetics for the regulatory system (Figure 1) can be described by the following ordinary differential equations which is in general called the Goodwin model (See [7,8]). Denote by $x_m(t)$ and $y_p(t)$ the concentrations of the species of Hes1 mRNA and Hes1 protein in a cell, respectively.

$ $$\begin{align}\label{ODE-hes1} \left\{ \begin{aligned} \frac{\mathrm{d}x_m(t)}{\mathrm{d}t}& = -\mu_m x_m(t)+\frac{\alpha_m}{1+\left(\frac{y_p(t)}{\bar{y}}\right)^h}, \\ \frac{\mathrm{d}y_p(t)}{\mathrm{d}t}& = -\mu_p y_p(t)+\alpha_p x_m(t), \end{aligned} \right. \end{align}$$ $

(1)

where $\mu_m$ and $\mu_p$ are the degradation rates of mRNA and protein, respectively, $\alpha_m$ the basal rate of initial transcription without feedback from protein, and $\alpha_p$ the production rate of protein, $\bar{y}$ the concentration of mRNA that reduces the rate of initial transcription to half of its basal value.

Since mRNA is typically transcribed in nucleus and then translated to cytoplasm for protein synthesis, there exists transcriptional and translational delays in the real reaction kinetics. For a better description of the intracellular processes, (1) was extended into the following model of delay differential equation (see, e. g., [12])

$ $$\begin{align}\label{DDE-hes1} \left\{ \begin{aligned} \frac{\mathrm{d}x_m(t)}{\mathrm{d}t}& = -\mu_m x_m(t)+\frac{\alpha_m}{1+\left(\frac{y_p(t-\tau_m)}{\bar{y}}\right)^h}, \\ \frac{\mathrm{d}y_p(t)}{\mathrm{d}t}& = -\mu_p y_p(t)+\alpha_p x_m(t-\tau_p), \end{aligned} \right. \end{align}$$ $

(2)

where $\tau_m$ and $\tau_p$ are transcriptional and translational delays, respectively. This model was able to simulate oscillatory solutions with the period and mRNA expression lag in good agreement with some experimental data. The transcriptional and translational time delay are usually assumed to be constant and it is then a natural problem to model the time delays in order to investigate the effect of the varying time delays on genetic regulatory systems.

For simplicity, we combine the reaction or transcription time with the translation time between the nucleus and the cytoplasm and call the combination is the diffusion time which will be denoted by $\tau$. In such a way $\tau_m$ in (2) is re-defined to be the mRNA transcription time and $\tau_p$ the protein diffusion time. That is, both are assumed to be the part of the time for diffusion processes, and the termination of the diffusion process incurs that of the inhibition of mRNA. That is, we assume that

(A1) $\tau = \tau_m+\tau_p$.

Now we drop the subscripts for the variables $x_m, \, y_p, \, \tau_m$ and $\tau_p$ in the following presentation, and turn to the modeling of the part of diffusion time due to spatial translation. Let $T$ to be the time to homogenize a newly produced solute (mRNA or protein) in the cell solvent, $\epsilon$ be the transcription time which will be assumed to be a constant. From Fick's first law of diffusion, we know that the diffusion time for a newly produced solute to homogenize in a cell is not necessarily a constant but is dependent on the concentration gradient. To model the diffusion time delay, let us suppose that we put a drop (or subtract a drop) of solute in a cup of solvent with concentration $x(t)$ which has taken account of the new drop. The diffusion time $T$ to homogenize depends on the homogenized concentration $x(t-\tau)$ of the solute in the cup without taking account of the new drop. Therefore the diffusion time to homogenize the new drop depends on the concentration difference $x(t)-x(t-\tau)$.

The next problem is how $T$ depends on the concentration difference $x(t)-x(t-\tau)$. In order to derive a reasonable modeling of $\tau$, we virtually assume that the distance of the central locations of concentrations $x(t)$ and $x(t-\tau(t))$ is $L$, where $L>0$ is small. It is known that the diffusion time $\tau$ is inversely proportional to the diffusion coefficient ($D$) and approximately satisfies $T = \frac{L^2}{2D}$. Then by Fick's first law we know that the diffusion flux (${J}$) is proportional to the existing concentration gradient under the assumption of steady state. Namely, $J$ can be approximated by $J = -D \frac{x(t)-x(t-\tau)}{L}$, where the positivity of $J$ is opposite to that of $x(t)-x(t-\tau)$. To obtain an approximation of $\tau$ we reduce $D$ from the aforementioned expressions of $T$ and $J$. Namely, $T $ can be approximated through

$ $$\begin{align} T& = -\frac{L(x(t)-x(t-\tau))}{2J}.\label{Travel-time} \end{align}$$ $

(3)

Then we have $\tau = T+\epsilon = -\frac{L(x(t)-x(t-\tau))}{2J}+\epsilon$. In the following, we assume that

(A2) there exist constants $\epsilon>0, \, c>0$ such that $\tau(t) = c (x(t)-x(t-\tau(t)))+\epsilon$.

Namely, we assume that the diffusion time is linearly dependent on the fluctuation of the concentrations of mDNA. We remark that $x(t)-x(t-\tau)$ may be negative which means that current concentration is less than the historical one and we interpret that in this case the transcription time is reduced by $c(x(t)-x(t-\tau))$ from $\epsilon$ with less concentration. Moreover, if the concentrations are in their equilibrium states, the diffusion time becomes the transcriptional delay $\epsilon$. We also remark that such type of state-dependent delay has appeared in position control [22] modeling the state-dependent traveling time of the echo, and in regenerative turning processes [11] modeling the cutting time of one round of turning.

Now we obtain the following state-dependent delay differential equations,

$ $$\begin{align}\label{SDDE-hes1} \left\{ \begin{aligned} \frac{\mathrm{d}x(t)}{\mathrm{d}t}& = -\mu_m x(t)+\frac{\alpha_m}{1+\left(\frac{y(t-\tau)}{\bar{y}}\right)^h}, \\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& = -\mu_p y(t)+\alpha_p x(t-\tau), \\ \tau(t)& = c(x(t)-x(t-\tau))+\epsilon. \end{aligned} \right. \end{align}$$ $

(4)

which provides another view on the models (1) and (2): If we assume that the transcription and translation of the newly produced solute is instant, namely, $T = 0$, then we have model (1) which is a system of ordinary differential equation. If we assume that the transcription and translation of the newly produced solute needs a constant positive time, then we have model (2) which is a system of delay differential equation with constant delay.

It was shown in [10] that system with the state-dependent delay in the form of $\tau(t) = r(x(t)-x(t-\tau(t)))$ with $r$ being inhomogeneous linear can be transformed into models with both discrete and distributed delay which greatly facilitated the qualitative analysis of such type of models [1, 10]. In the next section we investigate systems (4) based on the discussion of a more broad class of equations.

---

3. Regulatory dynamics with state-dependent delay

Regarding system (4) as a prototype model of a range of genetic regulatory dynamics, we consider the following system,

$ $$\begin{align}\label{SDDE-general-1} \left\{ \begin{aligned} \frac{\mathrm{d}x(t)}{\mathrm{d}t}& = -\mu_m x(t)+f(y(t-\tau)), \\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& = -\mu_p y(t)+ g(x(t-\tau)), \\ \tau(t)& = \epsilon+c(x(t)-x(t-\tau)), \end{aligned} \right. \end{align}$$ $

(5)

where $f, \, g: \mathbb{R}\rightarrow\mathbb{R}$ are three times continuously differentiable functions; $\mu_m$, $\mu_p$, $c$ and $\epsilon$ are positive constants. In the following we denote by $C([a, \, b];\mathbb{R}^N)$ with $a < b$ the space of continuously functions $f: [a, \, b]\rightarrow\mathbb{R}^N$ equipped with the supremum norm and by $C^1([a, \, b];\mathbb{R}^N)$ with $a < b$ the space of continuously differentiable functions. We assume that

B1) System (5) has at least one equilibrium $(r^*, \, \xi^*)$ for $(x, \, y)$ and there exists $\alpha_0>0$ and a neighborhood $D$ of $(r^*, \, \xi^*)$ in the space of continuous differentiable functions $C^1([-\alpha_0, \, 0];\mathbb{R}^2)$ such that for every initial data in $D$ for $(x, \, y)$ and the initial data $\tau_0\in(0, \, \alpha_0)$ satisfying the compatibility conditions

$ $$\begin{align}\label{compatibility} \left\{ \begin{aligned} x'(0)& = -\mu_m x(0)+f(y(-\tau_0)), \\ y'(0)& = -\mu_p y(0)+ g(x(-\tau_0)), \\ \tau_0& = \epsilon+c(x(0)-x(-\tau_0)), \end{aligned} \right. \end{align}$$ $

(6)

system (5) has a unique solution $(x, \, y, \, \tau)$ on $[0, \, +\infty)$.

B2) Let $D$ be as in (B1). For every solution $(x, \, y, \, \tau)$ of system (5) on $[0, \, +\infty)$ with initial data in $D$ for $(x, \, y)$ and initial data $\tau_0\in(0, \, +\infty)$, we have $1/c >\dot{x}(t)$ for all $t>0$.

We rewrite the equation for the delay in (4) as

$ $$\begin{align} \int_{t-\tau(t)}^t\frac{1-c\dot{x}(s)}{\epsilon}\mathrm{d}s = 1.\label{SDDE-eqn-5} \end{align}$$ $

(7)

Let $u(t) = \dot{x}(t)$ and introduce the following transformation inspired by [16]:

$ $$\begin{align}\label{smith-transform} \eta = \int_{0}^t\frac{1-cu(s)}{\epsilon}\mathrm{d}s, r(\eta) = x(t), \xi(\eta) = y(t), k(\eta) = \tau(t). \end{align}$$ $

(8)

Then we have

$ $$\begin{align} \label{const-delay} \left\{ \begin{aligned} \eta-1 & = \int_{0}^{t-\tau(t)}\frac{1-cu(s)}{\epsilon}\mathrm{d}s, \\ r(\eta-1) & = x(t-\tau(t)), \\ \xi(\eta-1) & = y(t-\tau(t)), \\ k(\eta) & = \epsilon+c(r(\eta)-r(\eta-1)). \end{aligned} \right. \end{align}$$ $

(9)

Note that we have

$ $$\begin{align}\label{derive} \left\{ \begin{aligned} \dot{x}(t)& = \frac{\mathrm{d}r}{\mathrm{d} \eta} \cdot \frac{\mathrm{d}\eta}{\mathrm{d}t} = \dot{r}(\eta)\cdot\frac{1-c\dot{x}}{\epsilon}, \\ \dot{y}(t) & = \frac{\mathrm{d}\xi}{\mathrm{d} \eta} \cdot \frac{\mathrm{d}\eta}{\mathrm{d}t} = \dot{\xi}(\eta)\cdot\frac{1-c\dot{x}}{\epsilon}, \end{aligned} \right. \end{align}$$ $

(10)

Then by (9) and (10), system (4) is transformed into

$ $$\begin{align}\label{tranformed-general} \left\{ \begin{aligned} \frac{\mathrm{d}r}{\mathrm{d} \eta}& = \epsilon\frac{-\mu_m r(\eta)+f(\xi(\eta-1))}{1-c(-\mu_m r(\eta)+f(\xi(\eta-1)))}, \\ \frac{\mathrm{d}\xi}{\mathrm{d} \eta}& = \epsilon\frac{-\mu_p \xi(\eta)+g(r(\eta-1))}{1-c(-\mu_m r(\eta)+f(\xi(\eta-1)))}, \\ k(\eta) & = \epsilon+c(r(\eta)-r(\eta-1)), \end{aligned} \right. \end{align}$$ $

(11)

where the assumption that $1/c>\dot{x}(t)$ for all $t\in\mathbb{R}$ has been converted by the first equation in (10) into $\dot{r}(\eta)>-\epsilon/c$. We remark that if we set $c = 0$ in system (11), we obtain

$ $$\begin{align}\label{tranformed-general-c-0} \left\{ \begin{aligned} \frac{1}{ \epsilon} \frac{\mathrm{d}r}{\mathrm{d} \eta}& = {-\mu_m r(\eta)+f(\xi(\eta-1))}\\ \frac{1}{ \epsilon} \frac{\mathrm{d}\xi}{\mathrm{d} \eta}& = {-\mu_p \xi(\eta)+g(r(\eta-1))}\\ k(\eta) & = \epsilon, \end{aligned} \right. \end{align}$$ $

(12)

which is equivalent to (5) with $c = 0$:

$ $$\begin{align}\label{SDDE-general-c-0} \left\{ \begin{aligned} \frac{\mathrm{d}x(t)}{\mathrm{d}t}& = -\mu_m x(t)+f(y(t-\tau)), \\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& = -\mu_p y(t)+ g(x(t-\tau)), \\ \tau(t)& = \epsilon, \end{aligned} \right. \end{align}$$ $

(13)

in the sense of the time-domain transformations leading to system (11).

Comparing system (11) with system (12), we see that both systems have the same set of equilibria. We also observe that, if $1-c\dot{x}(t)>0$ is uniformly bounded above for $t>0$, then (7) describes the situation that the diffusion process is slow when $\epsilon$ is large and is fast when $\epsilon$ is small. Moreover, from system (12) we know that the parameter $\epsilon$ may serve as a time scale which is highly desirable to investigate the intercelluar regulatory dynamics. In the subsequent sections we study how the dynamics of system (11) varies with respect to the parameter $\epsilon\in (0, \, +\infty)$. We remark that limiting profiles of periodic solutions as $\epsilon\rightarrow 0$ was investigated in [15] for the equation $ \epsilon\dot{x}(t) = -x(t) + f(x(t-1)). $

---

3.1. Linear stability analysis

Now we study the stability of the equilibrium of system (11). Note that the variable $k$ has been decoupled from the differential equations. Therefore, we assume that $(r^*, \, \xi^*)\in\mathbb{R}^2$ is the equilibrium of system (11) such that

$ $$\begin{align*} \left\{ \begin{aligned} -\mu_m r^*+f(\xi^*))& = 0, \\ -\mu_p\xi^*+g(r^*)& = 0. \end{aligned} \right. \end{align*}$$ $

We translate the equilibrium of (11) to the origin by letting $(u, \ v)$ be such that $r = u+r^*$, $\xi = v+\xi^*$ and let $F, \, G: \mathbb{R}\rightarrow\mathbb{R}$ satisfy that

$ $$\begin{align*} F(v(\eta-1)) = f(v(\eta-1)+\xi^*)-f(\xi^*), \;\;\;\; G(u(\eta-1)) = g(u(\eta-1)+r^*)-g(r^*). \end{align*}$$ $

Then we have $F'(0) = f'(\xi^*)$, $G'(0) = g'(r^*)$, and

$ $$\begin{align}\label{eqn-5-2} \left\{ \begin{aligned} \frac{1}{\epsilon}\frac{\mathrm{d}u}{\mathrm{d}\eta}& = \frac{-\mu_m u+F(v(\eta-1))}{1-c[-\mu_m u+F(v(\eta-1))]}, \\ \frac{1}{\epsilon}\frac{\mathrm{d}v}{\mathrm{d}\eta}& = \frac{-\mu_p v+G(u(\eta-1))}{1-c[-\mu_m u+F(v(\eta-1))]}. \end{aligned} \right. \end{align}$$ $

(14)

Linearization of (14) at $(0, \, 0)$ with $x = (u, \, v)$ is

$ $$\begin{align}\label{linearized} \dot{x}(\eta) = \epsilon M x(\eta)+\epsilon Nx(\eta-1), \end{align}$$ $

(15)

where $M = $$\begin{bmatrix} -\mu_m&0\\ 0&-\mu_p \end{bmatrix}$$ $ and $N = $$\begin{bmatrix} 0&f'(\xi^*)\\ g'(r^*)&0 \end{bmatrix}$$ $. The characteristic equation of (15) is $\det(\lambda I+\epsilon M-\epsilon N e^{-\lambda}) = 0, $ which is equivalent to

$ $$\begin{align}\label{charac-1} (\lambda+\epsilon\mu_m)(\lambda+\epsilon\mu_p)-\epsilon^2 f'(\xi^*)g'(r^*)e^{-2\lambda} = 0. \end{align}$$ $

(16)

Let $z = 2\lambda$, $\tau = 2\epsilon$ and $f'(\xi^*)g(r^*) = -h$, then the characteristic equation (16) becomes

$ $$\begin{align}\label{charac-2} z^2+\tau(\mu_m+\mu_p)z+\tau^2\mu_m\mu_p+\tau^2h e^{-z} = 0, \end{align}$$ $

(17)

which is the same equation investigated in [3]. Then by the results in [3] we have,

Lemma 3.1. Consider the characteristic equation (16), where $\mu_m$ and $\mu_p$ are positive constants. Then for every $\epsilon\in (0, \, +\infty)$, the equation

$ $$\begin{align*} \beta^2-{\epsilon}^2 (\mu_m\mu_p) = \epsilon(\mu_m+\mu_p)\beta \cot2\beta, \end{align*}$$ $

has a unique solution for $\beta$ in $(0, \, \frac{\pi}{2})$, denoted by $\beta(\epsilon)$ and the following statements hold.

ⅰ) If $ \mu_m\mu_p\geq -f'(\xi^*)g'(r^*)\geq 0, $ then the equilibrium of (11) is asymptotically stable.

ⅱ) If $ \mu_m\mu_p < -f'(\xi^*)g'(r^*), $ then there exists a unique $\epsilon_0\in (0, \, +\infty)$ such that $ (\mu_m+\mu_p) \beta(\epsilon_0)+\epsilon_0 f'(\xi^*)g'(r^*)\sin 2\beta(\epsilon_0) = 0, $ and $\pm i\beta(\epsilon_0)$ is a pair of purely imaginary eigenvalues. Moreover, for every $\epsilon\in (0, \, \epsilon_0)$, the equilibrium of (11) is asymptotically stable; for every $\epsilon\in (\epsilon_0, \, +\infty)$, the equilibrium of (11) is unstable, and there exists a sequence $\{\epsilon_k\}_{k = 0}^\infty$ of critical values of $\epsilon$ for which (16) has a corresponding sequence of purely imaginary eigenvalues $\{i(\beta(\epsilon_0)+k\pi)\}_{k = 1}^\infty$ where $ \epsilon_k = \frac{\epsilon_0(\beta(\epsilon_0)+k\pi)}{\beta(\epsilon_0)}\in S_k, $ with $S_0 = \left\{\lambda\in\mathbb{C}: \mathrm{Im}{\lambda}\in (0, \, \frac{\pi}{2})\right\}, $ $ S_k = \left\{\lambda\in\mathbb{C}: \mathrm{Im}{\lambda}\in ((k-\frac{1}{2})\pi, \, (k+\frac{1}{2})\pi)\right\}, \, k = 1, \, 2\cdots.$

Lemma 3.2. Consider the characteristic equation (16), where $\mu_m$ and $\mu_p$ are positive constants. If

$ \mu_m\mu_p < -f'(\xi^*)g'(r^*), $

then there exists a unique $\epsilon_0\in (0, \, +\infty)$ such that

$ (\mu_m+\mu_p) \beta(\epsilon_0)+\epsilon_0 f'(\xi^*)g'(r^*)\sin 2\beta(\epsilon_0) = 0, $

and $\pm i\beta(\epsilon_0)$ is a pair of purely imaginary eigenvalues. Moreover, $\beta(\epsilon_0)$ satisfies

$ $$\begin{align*} \tan 2\beta(\epsilon_0) = \frac{\sqrt{l}(\mu_m+\mu_p)}{1-l\mu_m\mu_p}, \end{align*}$$ $

where

$ $$\begin{align*} l = \frac{\mu_m^2+\mu_p^2+\sqrt{(\mu_m^2-\mu_p^2)^2+4f'^2(\xi^*)g'^2(r^*)} }{2(f'^2(\xi^*)g'^2(r^*)-\mu_m^2\mu_p^2)}, \epsilon_0 = \sqrt{l}\beta(\epsilon_0). \end{align*}$$ $

Proof. By Lemma 3.1 $\beta(\epsilon_0)$ satisfies

$ $$\begin{align}\label{eqn-tan-epsilon} \left\{ \begin{aligned} \beta^2(\epsilon_0)-{\epsilon_0}^2 (\mu_m\mu_p)-\epsilon_0(\mu_m+\mu_p)\beta(\epsilon_0) \cot2\beta(\epsilon_0)& = 0, \\ (\mu_m+\mu_p) \beta(\epsilon_0)+\epsilon_0 f'(\xi^*)g'(r^*)\sin 2\beta(\epsilon_0)& = 0. \end{aligned} \right. \end{align}$$ $

(18)

Eliminating the trigonometric terms in (18), we obtain that

$ \epsilon_0^4(\mu_m^2\mu_p^2-f'^2(\xi^*)g'^2(r^*)))+\epsilon_0^2(\mu_m^2+\mu_p^2)\beta^2(\epsilon_0)+\beta^4(\epsilon_0) = 0, $

which can be regarded as a quadratic equation of $\epsilon_0^2$ and has a unique nonnegative root since $ \mu_m\mu_p\leq -f'(\xi^*)g'(r^*).$ We have

$ $$\begin{align*} \epsilon_0^2 = l\beta^2(\epsilon_0), l = \frac{\mu_m^2+\mu_p^2+\sqrt{(\mu_m^2-\mu_p^2)^2+4f'^2(\xi^*)g'^2(r^*)} }{2(f'^2(\xi^*)g'^2(r^*)-\mu_m^2\mu_p^2)}. \end{align*}$$ $

Bringing $\epsilon_0 = \sqrt{l}\beta(\epsilon_0)$ into the first equation of (18), we have

$ $$\begin{align*} \tan 2\beta(\epsilon_0) = \frac{\sqrt{l}(\mu_m+\mu_p)}{1-l\mu_m\mu_p}. \end{align*}$$ $

---

3.2. Local Hopf bifurcation

In this section we use multiple time scale method to find the normal form of system (11) at the first critical value $\epsilon^*$ of $\epsilon$. Namely, we let $\epsilon^* = \epsilon_0$ where $\epsilon_0$ is defined at Lemma 4.1. Since we are interested in the dynamics when $\epsilon$ varies in $(0, \, \infty)$, we study the directions of the Hopf bifurcation and the stability of the bifurcating periodic solutions near the equilibrium. We first check the transversality condition. Put $\lambda = \alpha+i\beta$ with $\alpha, \, \beta\in\mathbb{R}$ in the characteristic equation (16) and we obtain

$ $$\begin{align}\label{re-um-3} \begin{cases} \begin{aligned} \alpha^2-\beta^2+\epsilon(\mu_m+\mu_p)\alpha+\epsilon^2\mu_m\mu_p-\epsilon^2f'(\xi^*)g'(r^*)e^{-2\alpha}\cos 2\beta& = 0, \\ 2\alpha\beta+ \epsilon(\mu_m+\mu_p)\beta+\epsilon^2f'(\xi^*)g'(r^*)e^{-2\alpha}\sin2\beta& = 0. \end{aligned} \end{cases} \end{align}$$ $

(19)

By Implicit function theorem, we can show that $\alpha$ and $\beta$ are continuously differentiable with respect to $\epsilon$ near $(\alpha, \, \beta, \, \epsilon) = (0, \, \beta(\epsilon_0), \, \epsilon_0)$. We take derivatives respect to $\epsilon$ on both sides of each of the equations at (19) and then let $\alpha = 0$. We obtain,

$ {\small $$\begin{align}\label{re-um-2} \left\{ \begin{aligned} \frac{\mathrm{d}\alpha}{\mathrm{d} \epsilon}&\left(\epsilon^*(\mu_m+\mu_p)+2{\epsilon^*}^2 f'(\xi^*)g'(r^*)\cos 2\beta\right)+ \frac{\mathrm{d}\beta}{\mathrm{d} \epsilon}\left(-2\beta+2{\epsilon^*}^2 f'(\xi^*)g'(r^*)\sin 2\beta\right)\\ & = -2\epsilon^*\mu_m\mu_p+2f'(\xi^*)g'(r^*) \cos2\beta, \\ \frac{\mathrm{d}\alpha}{\mathrm{d} \epsilon}&\left(2\beta-2{\epsilon^*}^2 f'(\xi^*)g'(r^*)\sin 2\beta\right)+ \frac{\mathrm{d}\beta}{\mathrm{d} \epsilon} \left(\epsilon^*(\mu_m+\mu_p)+2{\epsilon^*}^2 f'(\xi^*)g'(r^*)\cos 2\beta\right)\\ & = -(\mu_m+\mu_p)\beta-2\epsilon^*f'(\xi^*)g'(r^*) \sin2\beta. \end{aligned} \right. \end{align}$$ } $

(20)

Noticing that if $\alpha = 0$, $\epsilon = \epsilon^*$, (19) implies that

$ $$\begin{align}\label{re-um-4} \begin{cases} {\epsilon^*}^2f'(\xi^*)g'(r^*) \cos2\beta& = {\epsilon^*}^2\mu_m\mu_p-\beta^2, \\ {\epsilon^*}^2f'(\xi^*)g'(r^*) \sin2\beta& = -{\epsilon^*}(\mu_m+\mu_p)\beta. \end{cases} \end{align}$$ $

(21)

Combining (20) and (21), we obtain that

$ $$\begin{align}\label{phi-direction} \frac{\mathrm{d}\alpha}{\mathrm{d} \epsilon}\, |_{\epsilon = \epsilon^*, \, \lambda = \pm i\beta} = \frac{\frac{2\beta^2}{\epsilon^*}({\epsilon^*}^2(\mu_m^2+\mu_p^2)+2\beta^2)}{(\epsilon^*(\mu_m+\mu_p)+2{\epsilon^*}^2\mu_m\mu_p-2\beta^2)^2+(2\beta+2\beta\epsilon^*(\mu_m+\mu_p))^2}>0. \end{align}$$ $

(22)

To illustrate the stability of system (5), we have,

Theorem 3.3. Consider system (5) with an equilibrium $(r^*, \, \xi^*)$, where $f, \, g: U\subset C^1([-\alpha_0, \, 0];\mathbb{R}^2)\rightarrow\mathbb{R}$ are three times continuously differentiable and $U$ an open neighborhood of $(r^*, \, \xi^*)$; $\mu_m, \, \mu_p$ and $c$ are positive constants. Assume (B1) and (B2) hold. Then the following statements hold.

ⅰ) If $ \mu_m\mu_p\geq -f'(\xi^*)g'(r^*)\geq 0, $ then the equilibrium of system (5) is asymptotically stable.

ii) If $ \mu_m\mu_p < -f'(\xi^*)g'(r^*), $ then there exists a unique $\epsilon_0\in (0, \, +\infty)$ such that for every $\epsilon\in (0, \, \epsilon_0)$, the equilibrium $(r^*, \, \xi^*)$ is asymptotically stable; for every $\epsilon\in (\epsilon_0, \, +\infty)$, the equilibrium $(r^*, \, \xi^*)$ is unstable, and there exists a sequence $\{\epsilon_k\}_{k = 0}^\infty$ of critical values of $\epsilon$ for which system (5) undergoes Hopf bifurcation.

Proof. By the technique of formal linearization [4,23], we know that the equilibrium stability of system (5) is equivalent to that of system (13) which has the same equilibrium and characteristic equation for system (11). Then by Lemma 3.1, the stability of $(r^*, \, \xi^*)$ of system (11) implies that of system (5). Existence of Hopf bifurcation follows from the Hopf bifurcation theorem established in [1].

Let $x(\eta) = (u(\eta), \, v(\eta))$ and system (11) is written

$ \frac{1}{\epsilon}\dot{x}(\eta) = $$\begin{bmatrix} S((x(\eta), x(\eta-1))\\ H((x(\eta), x(\eta-1)) \end{bmatrix}$$ , $

where $(S, \, H)$ denotes the right hand side of (14). We have the Taylor expansion of (14) up to the cubic terms,

$ {\small $$\begin{align}\label{Taylor} \left\{ \begin{aligned} \frac{1}{\epsilon}\frac{\mathrm{d}u}{\mathrm{d}\eta} = &-\mu_mu(\eta)+f'(\xi^*)v(\eta-1)+ c\mu_{m}^2 u^2(\eta)-2c\mu_mf'(\xi^*)u(\eta)v(\eta-1)\\ &+\left(\frac{1}{2}f''(\xi^*)+cf'^2(\xi^*)\right)v^2(\eta-1)+3c^2\mu_m^2f'(\xi^*) u^2(\eta)v(\eta-1) \\ &-c^2\mu_m^3 u^3(\eta) -(c\mu_m f''(\xi^*)+3c^2\mu_m f'^2(\xi^*)) u(\eta) v^2(\eta-1) \\&+ \left(\frac{1}{6}f'''(\xi^*)+cf''(\xi^*) f'(\xi^*)+c^2f'^3(\xi^*)\right)v^3(\eta-1), \\ \frac{1}{\epsilon}\frac{\mathrm{d}v}{\mathrm{d}\eta} = &-\mu_p v(\eta)+g'(r^*)u(\eta-1)+\frac{1}{2}\left[-2c\mu_p f'(\xi^*)v(\eta) v(\eta-1)+2c\mu_m\mu_p u(\eta) v(\eta)\right.\\ &\left. +g''(r^*)u^2(\eta-1)+2cf'(\xi^*)g'(r^*)u(\eta-1)v(\eta-1)-2c\mu_m g'(r^*)u u(\eta-1) \right]\\ &+\frac{1}{6}\left[g'''(r^*) u^3(\eta-1)-6c^2\mu_m^2\mu_p u^2(\eta) v(\eta)+6c^2\mu_m^2g'(r^*) u^2(\eta) u(\eta-1)\right.\\ &\left.-3c\mu_mg''(r^*)u(\eta) u^2(\eta-1)+3cg''(r^*)f'(\xi^*)u^2(\eta-1)v(\eta-1)\right.\\ &\left.-3(c\mu_p f''(\xi^*)+2c^2\mu_pf'^2(\xi^*)) v(\eta) v^2(\eta-1)+3(cf''(\xi^*)g'(r^*)\right.\\ &\left. +2c^2f'^2(\xi^*)g'(r^*))u(\eta-1)v^2(\eta-1)+12c^2\mu_m\mu_p f'(\xi^*)u(\eta) v(\eta) v(\eta-1)) \right.\\ &\left.-12c^2\mu_mg'(r^*)f'(\xi^*) u(\eta) u(\eta-1)v(\eta-1) \right]. \end{aligned} \right. \end{align}$$ } $

(23)

---

3.3. Normal form computation via multiple time scale

In this section we compute the normal form of system (14) using its Taylor expansion at (23). We seek a uniform second-order approximate solution in the form

$ $$\begin{align}\label{mts-1} x(\eta;\, s) = sx_1(T_0, \, T_2)+s^2x_2(T_0, \, T_2)+s^3x_3(T_0, \, T_2)+\cdots, \end{align}$$ $

(24)

where $T_0 = \eta$, $T_2 = s^2\eta$, and $s$ is a nondimensional book keeping parameter. In the following, we write $x_i = (u_i, \, v_i), \, i = 1, \, 2, \, \cdots$. It is not necessary to relate the approximate solution with the time scale $T_1 = s\eta$ because a second order approximation will require independence of $T_1$ in the solution (See, e.g., [14], page 166).

$ $$\begin{align}\label{mts-2} \frac{\mathrm{d}}{\mathrm{d}\eta} = \frac{\partial}{\partial T_0}+s^2\frac{\partial }{\partial T_2}+\cdots = D_0+s^2 D_2+\cdots. \end{align}$$ $

(25)

By (24) and (25) we have

$ $$\begin{align}\label{mts-3} x(\eta-1)& = x(T_0-1, \, T_2-s^2) = \sum\limits_{n = 1}^\infty s^n x_n(T_0-1, \, T_2)-s^3 D_2 x_1(T_0-1, \, T_2)\cdots. \end{align}$$ $

(26)

Let $\epsilon = \epsilon^*+s^2\delta$ where $\epsilon^*$ is the critical value for Hopf bifurcation and $\delta$ a small detuning parameter. Bringing (24) and (26) into (23), we have the left hand side been transformed into

$ $$\begin{align*} &\frac{1}{(\epsilon^*+s^2\delta)}(D_0+s^2 D_2)\sum^{\infty}_{n = 1}s^n x_n(T_0, \, T_2) \\ = &\frac{1}{(\epsilon^*+s^2\delta)}\left[s D_0 x_1(T_0, \, T_2)+s^2D_0x_2(T_0, \, T_2)\right.\\ & \left.+ s^3D_0 x_3(T_0, \, T_2)+ s^3 D_2 x_1(T_0, \, T_2)+s^4 D_2 x_2(T_0, \, T_2)+\cdots\right], \end{align*}$$ $

and the right hand side becomes

$ $$\begin{align*} M\sum^{\infty}_{n = 1}s^n x_n (T_0, \, T_2) +N\left(\sum\limits_{n = 1}^\infty s^n x_n(T_0-1, \, T_2)-s^3D_2 x_1(T_0-1, \, T_2)+\cdots\right)+\cdots. \end{align*}$$ $

Comparing the like powers of $s$ up to the cubic in the transformed equation of (23) yield the following equations, where for notational simplicity we omitted the time scales of the variables except for the delayed ones and use dashed lines to separate entries of the column matrices.

$ $$\begin{align} &D_0 x_1-\epsilon^*Mx_1-\epsilon^*Nx_1(T_0-1, \, T_2) = 0, \label{s-eqn-1} \end{align}$$ $

(27)

(28)

(29)

We know that when $\epsilon = \epsilon^*$ system (14) undergoes Hopf bifurcation and (29) has periodic solution of the form

$ x_1 = A(T_2)\theta e^{iw^*T_0}+\bar{A}(T_2)\bar{\theta}e^{-iw^*T_0}, $

where $iw^*$ is a purely imaginary eigenvalue with eigenvector $\theta = [\theta_1, \, \theta_2]^T\neq 0$ which satisfies

$ $$\begin{align*} \left(iw^*I+\epsilon^*\begin{bmatrix} \mu_m&0\\ 0&\mu_p \end{bmatrix} -\epsilon^*\begin{bmatrix} 0&f'(\xi^*)e^{-iw^*}\\ g'(r^*)e^{-iw^*}&0 \end{bmatrix}\right) \left(\begin{matrix}\theta_1\\ \theta_2 \end{matrix}\right)& = 0. \end{align*}$$ $

Leting $\theta_1 = 1$ in the above equation we have

$ $$\begin{align*} \theta = \left[ \begin{matrix} 1\\ \frac{e^{i w^*}(i w^*+\epsilon^*\mu_m)}{\epsilon^* f'(\xi^*)} \end{matrix} \right]. \end{align*}$$ $

Substitute $x_1$ into (28), we have

(30)

Equation (28) has a particular solution of the form

$x_2 = aA^2e^{i2w^*T_0}+bA\bar{A}+\bar{a}\bar{A}^2e^{-i2 w^*T_0}, $

where $a = (a_1, \, a_2)\in\mathbb{C}^2, \, b = (b_1, \, b_2)\in\mathbb{C}^2.$ Bringing the undetermined form of $x_2$ into (30), we obtain

$ $$\begin{align}\label{eqn-5-12} &(aA^2e^{i2w^*T_0}\cdot i2w^*+\bar{a}\bar{A}^2e^{-i2w^*T_0}(-i2w^*))\\&+\epsilon^*\left[\begin{matrix} \mu_m&0\notag\\ 0&\mu_p \end{matrix}\right] (aA^2e^{i2w^*T_0}+bA\bar{A}+\bar{a}\bar{A}^2e^{-i2w^*T_0})\notag\\ &-\epsilon^*\left[\begin{matrix} 0&f'(\xi^*)\\ g'(r^*) &0 \end{matrix}\right] (aA^2e^{i2w^*(T_0-1)}+bA \bar{A}+\bar{a} \bar{A}^2 e^{-i2w^*(T_0-1)})\notag \end{align}$$ $

(31)

(32)

Equating the coefficients of the terms with $e^{i2w^*T_0}$ from both sides of (31), we obtain

$ $$\begin{align}\label{eqn-5-13} \left\{ \begin{aligned} &(i2w^*)A^2a_1+\epsilon^*\mu_mA^2a_1-\epsilon^*f'(\xi^*)A^2e^{-i2w^*}a_2 \\ = & \epsilon^*c\mu_m^2A^2-2\epsilon^*c\mu_mf'(\xi^*)A^2\theta_2e^{-iw^*} +\frac{\epsilon^*}{2}(f''(\xi^*)+2cf'^2(\xi^*))A^2\theta_2^2e^{-i2w^*}, \\&(i2w^*)A^2a_2+\epsilon^*\mu_pA^2a_2-\epsilon^*g'(r^*)A^2e^{-i2w^*}a_1 \\ = & -\epsilon^*c\mu_pf'(\xi^*)\theta^2_2A^2e^{-iw^*}+\epsilon^*c\mu_m\mu_pA^2\theta_2 +\frac{\epsilon^*}{2}g''(r^*)A^2e^{-2iw^*}\\ &+\epsilon^*cf'(\xi^*)g'(r^*)A^2\theta_2e^{-2iw^*} -\epsilon^*c\mu_mg'(r^*)A^2e^{-iw^*}. \end{aligned} \right. \end{align}$$ $

(33)

Equation (33) has a solution for $a$ if $i2w^*$ is not an eigenvalue of (15).

(B3) $i2w^*$ is not an eigenvalue of (15).

We have

$ $$\begin{align}\label{eqn-b-1-2} \left\{ \begin{aligned} a_1 = &\frac{\epsilon^*}{ (i2w^*+\epsilon^*\mu_m)(i2w^*+\epsilon^*\mu_p) -{\epsilon^*}^2f'(\xi^*)g'(r^*)e^{-i4w^*} }\\ &\times\bigg[\left( (c\mu_m^2-2c\mu_mf'(\xi^*)\theta_2e^{-iw^*} +\frac{1}{2}(f''(\xi^*)+2cf'^2(\xi^*))\theta_2^2e^{-i2w^*} \right) \\ & \times (i2w^*+\epsilon^*\mu_p) +\bigg(-c\mu_pf'(\xi^*)\theta^2_2 e^{-iw^*}+ c\mu_m\mu_p \theta_2 +\frac{1}{2}g''(r^*) e^{-2iw^*}\\ & + cf'(\xi^*)g'(r^*) \theta_2e^{-2iw^*}-c\mu_mg'(r^*) e^{-iw^*}\bigg)(\epsilon^*f'(\xi^*)e^{-i2w^*}) \bigg], \\ a_2 = &\frac{\epsilon^*}{ (i2w^*+\epsilon^*\mu_m)(i2w^*+\epsilon^*\mu_p) -{\epsilon^*}^2f'(\xi^*)g'(r^*)e^{-i4w^*} }\\ &\times\bigg[\left( (c\mu_m^2-2c\mu_mf'(\xi^*)\theta_2e^{-iw^*} +\frac{1}{2}(f''(\xi^*)+2cf'^2(\xi^*))\theta_2^2e^{-i2w^*} \right) \\ & \times (\epsilon^*g'(r^*)e^{-i2w^*} ) +\bigg(-c\mu_pf'(\xi^*)\theta^2_2 e^{-iw^*}+ c\mu_m\mu_p \theta_2 +\frac{1}{2}g''(r^*) e^{-2iw^*}\\ & + cf'(\xi^*)g'(r^*) \theta_2e^{-2iw^*}-c\mu_mg'(r^*) e^{-iw^*}\bigg)(i2w^*+\epsilon^*\mu_m) \bigg]. \end{aligned} \right. \end{align}$$ $

(34)

Equating the coefficients of the $A\bar{A}$ terms from both sides of (31), we obtain,

$ $$\begin{align*} \left\{ \begin{aligned} \epsilon^*\mu_mb_1-\epsilon^*f'(\xi^*)b_2 = &\, \, 2\epsilon^*c\mu_m^2-2\epsilon^*c\mu_mf'(\xi^*)(\bar{\theta_2}e^{iw^*}+\theta_2e^{-iw^*})\\ &\, \, +\epsilon^*(f''(\xi^*)+2cf'^2(\xi^*))\theta_2\bar{\theta_2}, \\ -\epsilon^*g'(r^*)b_1+\epsilon^*\mu_pb_2 = &\, \, -\epsilon^*c\mu_pf'(\xi^*)\theta_2\bar{\theta_2}(e^{iw^*}+e^{-iw^*})+\epsilon^*c\mu_m\mu_p(\theta_2+\bar{\theta_2})\\ &+\epsilon^*g''(r^*)+\epsilon^*cf'(\xi^*)g'(r^*)(\theta_2+\bar{\theta_2})-\epsilon^*c\mu_mg'(r^*)\\ &\times(e^{iw^*}+e^{-iw^*}). \end{aligned} \right. \end{align*}$$ $

Noticing that $\theta_2 = \frac{e^{i w^*}(i w^*+\epsilon^*\mu_m)}{\epsilon^* f'(\xi^*)}$ and hence $\theta_2e^{-iw^*}+\bar{\theta}_2e^{iw^*} = \frac{2\mu_m}{f'(\xi^*)}$, we have

$ $$\begin{align*} \left\{ \begin{aligned} b_1 = &\frac{1}{{\epsilon^*}^2f'^2(\xi^*)(\mu_m\mu_p-f'(\xi^*)g'(r^*))} \bigg[ \mu_pf''(\xi^*){w^*}^2+\mu_pf''(\xi^*){\epsilon^*}^2\mu_m^2\\ &+2f'^2(\xi^*)c\mu_p{w^*}^2-2f'^2(\xi^*)c\mu_p{w^*}^2\cos w^*-2c(f'^3(\xi^*)g'(r^*)\\ &+\mu_m\mu_pf'^2(\xi^*)){\epsilon^*}{w^*}\sin w^* +f'^3(\xi^*)g''(r^*){\epsilon^*}^2 \bigg], \\ b_2 = &\frac{1}{{\epsilon^*}^2f'^2(\xi^*)(\mu_m\mu_p-f'(\xi^*)g'(r^*))} \bigg[ \mu_mf'^2(\xi^*)g''(r^*){\epsilon^*}^2+ f''(\xi^*)g'(r^*){w^*}^2\\ & +f''(\xi^*)g'(r^*)\mu_m^2{\epsilon^*}^2+2cf'^2(\xi^*)g'(r^*){w^*}^2-2c\mu_m\mu_pf'(\xi^*){w^*}^2\cos w^*\\ &-2c(\mu_mf'^2(\xi^*)g'(r^*)+\mu_m^2\mu_pf'(\xi^*)){\epsilon^*}{w^*}\sin w^* \bigg]. \end{aligned} \right. \end{align*}$$ $

Substituting the expressions of $x_1$ and $x_2$ into (4.9), we obtain,

(35)

Note that $iw^*$ is an eigenvalue of the homogenous system corresponding to (35) and the right hand side of the nonhomogeneous system (35) has terms with $e^{iw^*T_0}$. To remove possible secular terms in the solutions, we seek a particular solution of the form $x_3(T_0, \, T_2) = \phi(T_2)e^{iw^*T_0}$, where the existence of $\phi$ requires solvability of the algebraic equation

$ $$\begin{align}\label{eqn-s-515} \left[iw^*I-\epsilon^*M-\epsilon^*Ne^{-iw^*}\right]\phi = \chi. \end{align}$$ $

(36)

where $\chi$ stands for the coefficients of the $e^{iw^*T_0} $ term in (35). System (36) is solvable for $\phi$ if and only if $\chi$ is orthogonal to the null space of the conjugate transpose of the coefficient matrix. That is, for every vector $d$ which satisfies

$ $$\begin{align*} \left[-iw^*I-\epsilon^*M^T-\epsilon^*N^Te^{iw^*}\right]d = 0, \end{align*}$$ $

we require

$ $$\begin{align}\label{eqn-s-516} \bar{d}^T\chi = 0, \end{align}$$ $

(37)

which is an ordinary differential equation of $A$ and is the normal form. We can obtain $d$ as following with the condition $\bar{d}^T\theta = 1$ imposed for uniqueness:

$ $$\begin{align*} d = \frac{1}{-2iw^*+\epsilon^*(\mu_m+\mu_p)} \begin{bmatrix} {-iw^*+\epsilon^*\mu_p}\\ \epsilon^*e^{iw^*}f'(\xi^*) \end{bmatrix}. \end{align*}$$ $

Bringing the expressions of $x_1$ and $x_2$ into system (35), we have the following expression for $\chi:$

(38)

By (37) and (38) and noticing that $\bar{d}^T\theta = 1$ and

$ iw^*\theta-\epsilon^*M\theta-\epsilon^*N\theta e^{-iw^*} = 0, $

we have the following normal form:

(39)

We arrive at,

Theorem 3.4. Assume that (B1), (B2) and (B3) hold. The normal form of system (5) near its equilibrium is (39).

We remark that if the state-dependent translation time $T$ is not included in the diffusion time $\tau$, namely, $c = 0$, the normal form at (39) becomes the normal form for system (5) with constant delay $\epsilon$:

$ $$\begin{array}{l} A' = \frac{{i{w^*}{e^{i{w^*}}}\delta }}{{{\epsilon ^*}({e^{i{w^*}}} + {\epsilon^*}{{\bar d}^T}N\theta )}}A\\ + \frac{{{\epsilon^*}}}{{({e^{i{w^*}}} + {\epsilon^*}{{\bar d}^T}N\theta )}}{{\bar d}^T}\left[ \begin{array}{l} f''({\xi ^*})({a_2}{{\bar \theta }_2} + {b_2}{\theta _2}) + \frac{1}{2}f'''({\xi ^*})\theta _2^2{{\bar \theta }_2}\\ g''({r^*})({a_1} + {b_1}) + \frac{1}{2}g'''({r^*}) \end{array}$$ \right]{A^2}\bar A. \end{array} $

(40)

---

4. Hes1 system (4) and numerical simulation

In this section, we consider the prototype system (4) and illustrate the results of section 3 by numerical simulation. In the following, a vector is said to be positive if all coordinates are positive. To have assumptions (B1) and (B2) hold for system (4), we show that,

Theorem 4.1. Consider system (4) with positive parameters $\mu_m$, $\mu_p$, $\alpha_m$, $\alpha_p$, $h$, $\bar{y}$, $\epsilon$ and $c$. Suppose that $\alpha_m\leq\frac{1}{c}$. Then system (4) has a positive equilibrium $(x^*, \, y^*)$ and there exists a neighborhood $U$ of $(x^*, \, y^*)\subset C^1([-\alpha_0, \, 0];\mathbb{R}^2)$ with $\alpha_0>0$ such that, for every initial data for $(x, \, y)\in U$, and initial data for $\tau$ in $(0, \, \alpha_0)$ satisfying the compatibility conditions at (6), there exists a unique positive solution $(x, \, y, \, \tau)$ on $[0, \, +\infty)$ with $\dot{x}(t) < \alpha_m$ for every $t\geq 0$.

Proof. Since all parameters of system (4) are positive, direct computation shows that there exists a positive equilibrium. Then there exists a neighborhood $U$ of $(x^*, \, y^*)\subset C^1([-\alpha_0, \, 0];\mathbb{R}^2)$ such that every $(x, \, y)\in U$ is positive and $\dot{x}(s) < \frac{1}{c}$ for all $s\in [-\alpha_0, \, 0]$.

Consider the initial value problem associated with system (4). We obtain by the work of [23] that the initial value problem has a unique solution $(x, \, y, \, \tau)$ on $[0, \, t_e)$ where $[0, \, t_e)$ is the maximal existence interval with $t_e>0$ or $t_e = \infty$ (c.f., [5]).

Notice that if $\tau(s) = 0$ for some $s\in (0, \, t_e)$ then from the third equation of system (4) we have $\tau(s) = \epsilon>0$ which contradicts $\tau(s) = 0$. Therefore, we have $\tau(s)>0$ for every $s\in [0, \, t_e)$.

Next we show that $\dot{x}(t) < \frac{1}{c}$ for every $t\in (0, \, t_e)$. Assume that there exists a minimal $s_0\in (0, \, t_e)$ such that $\dot{x}(s_0) = \frac{1}{c}$. Then by the third equation of system (4) we have for every $t\in [0, \, s_0]$,

$ \dot{\tau}(t) = \frac{\dot{x}(t)-\dot{x}(t-\tau(t))}{\frac{1}{c}-\dot{x}(t-\tau(t))} < 1. $

It follows that for every $t\in [0, \, s_0]$ we have $-\tau_0 < t-\tau(t) < t$. Then we have

$ \dot{x}(s) = -\mu_m x(s)+\frac{\alpha_m}{1+\left(\frac{y(s-\tau(s))}{\bar{y}}\right)^h} < \frac{\alpha_m}{1+\left(\frac{y(s-\tau(s))}{\bar{y}}\right)^h} < \alpha_m, $

which leads to $\frac{1}{c} < \alpha_m$. This is a contradiction. Therefore, $\dot{x}(t) < \frac{1}{c}$ for every $t\in (0, \, t_e)$.

Next we show that $x(t)$ and $y(t)$ are positive on $[0, \, t_e)$. Suppose, for contradiction, that there exists $t_0\in (0, \, t_e)$ which is the minimal time such that $x(t)$ or $y(t)$ vanishes. If $x(t_0) = 0$, then we have

$ \dot{x}(t_0) =-\mu_m x(t_0)+\frac{\alpha_m}{1+\left(\frac{y(t_0-\tau(t_0))}{\bar{y}}\right)^h} = \frac{\alpha_m}{1+\left(\frac{y(t_0-\tau(t_0))}{\bar{y}}\right)^h}>0, $

which means that there exists $t_0^*\in (0, \, t_0)$ such that $x(t_0^*) < 0$. This contradicts the minimality of $t_0$.

If $y(t_0) = 0$, then we have

$ \dot{y}(t_0) =-\mu_p y(t_0)+\alpha_px(t_0-\tau(t_0)) = \alpha_px(t_0-\tau(t_0))>0, $

which means that there exists $t_0^{**}\in (0, \, t_0)$ such that $y(t_0^{**}) < 0$. This contradicts the minimality of $t_0$. That is, $x(t)$, $y(t)$ and $\tau(t)$ are positive on $[0, \, t_e)$ and for every $t\in [0, \, t_e)$ we have

$ $$\begin{align}\label{est-x} \dot{x}(t) =-\mu_m x(t)+\frac{\alpha_m}{1+\left(\frac{y(t-\tau(t))}{\bar{y}}\right)^h} < \frac{\alpha_m}{1+\left(\frac{y(t-\tau(t))}{\bar{y}}\right)^h} < \alpha_m. \end{align}$$ $

(41)

It follows that

$ $$\begin{align}\label{est-y} \dot{y}(t) =-\mu_p y(t)+\alpha_px(t-\tau(t)) < \alpha_px(t-\tau(t)) < \alpha_p(\alpha_m t+x(0)). \end{align}$$ $

(42)

By (41) and (42), and by the third equation of system (4), $(x, \, y, \, \tau)$ is uniformly bounded on $[-\alpha_0, \, t_e)$ and hence extendable beyond $t_e$, which contradicts the maximality of $t_e$ and we have $t_e = \infty$.

Now we turn to assumption (B3). Let $\mu_m = 0.03$, $\mu_p = 0.04$, $\alpha_m = 35$, $\alpha_p = 10$, $\bar{y} = 1200$, $h = 5$. Then we have the equilibria $(r^*, \, \xi^*) = (11.97050076, \, 2992.625189)$, $\omega^* = 0.47038322$. We have $f'(\xi^*) =-0.00059384374$, $g'(r^*) = 10$. Then we have

$ \mu_m\mu_p = 1.2\times 10^{-3} < 5.9384374\times 10^{-3} =-f'(\xi^*)g'(r^*). $

Then by Lemma 3.1 and by Theorem 3.3, there exists a unique $\epsilon_0 = 6.86216245$, such that for every $\epsilon\in (0, \, \epsilon_0)$, the equilibria $(r^*, \, \xi^*)$ is stable; for every $\epsilon\in [\epsilon_0, \, \infty)$ the equilibria $(r^*, \, \xi^*)$ is unstable and $\epsilon_0$ is a critical value of $\epsilon$ for which system (4) undergoes Hopf bifurcation. For $\epsilon = \epsilon_0$ and $\lambda = 2iw^*$, we have

$ $$\begin{align*} & (\lambda+\epsilon\mu_m)(\lambda+\epsilon\mu_p)-\epsilon^2f'(\xi^*)g'(r^*)e^{-2\lambda}\big|_{\, \lambda = 2iw^*, \, \epsilon = \epsilon_0}\\ = &-0.9140361052+0.1856539388\, i\neq 0. \end{align*}$$ $

That is, (B3) holds. By Theorem 3.4, the normal form of system (4) near $(r^*, \, \xi^*)$ is,

$ $$\begin{align*} A' = &(0.01841158248+0.04829902976\, i)\delta A \\ & + \left((2.114544332\, c^2+0.0008578251748\, c-0.001233336633)\right.\\ & \left.-i(1.469928514\, c^2+0.002534237744\, c-0.003599996653)\right)A^2\bar{A}, \end{align*}$$ $

where the real part of the coefficient of the cubic term $ A^2\bar{A}$ is

$ 2.114544332\, c^2+0.0008578251748\, c-0.001233336633, $

which is negative for $0 < c < c_0: = 0.02394886242, $ and is positive for $c>c_0.$

If $0 < c < c_0$, system (4) undergoes supercritical Hopf bifurcation near $(r^*, \, \xi^*)$ as $\epsilon$ crosses the critical value $\epsilon_0$. Namely, if $0 < \epsilon < \epsilon_0$, the equilibrium is stable (see Figure 2(a)); if $\epsilon>\epsilon_0$, there is a stable periodic solution near $(r^*, \, \xi^*)$ (see Figure 2(b)).

If $c>c_0$, system (4) undergoes subcritical Hopf bifurcation near $(r^*, \, \xi^*)$ as $\epsilon$ crosses the critical value $\epsilon_0$. Namely, if $0 < \epsilon < \epsilon_0$, the equilibrium is stable (see Figure 3(a)) and there is a unstable periodic solution near $(r^*, \, \xi^*)$. See Figure 3(a), where the equilibrium $(r^*, \, \xi^*) = (11.97050076, \, 2992.625189)$ is asymptotically stable with $\epsilon = \epsilon_0-\delta$, $c = c_0+0.001$ (the solid curve); when initial value is far enough from the equilibrium $(r^*, \, \xi^*) = (11.97050076, \, 2992.625189)$, solution with nonpositive initial values (the dashed curve) may converge to a negative equilibrium. If $\epsilon>\epsilon_0$, the equilibrium $(r^*, \, \xi^*)$ is unstable (see Figure 3(b)).

---

5. Concluding remarks

Starting from model (5) of differential equations with threshold type state-dependent delay, we obtained model (11) with constant delay equivalent to (5) under assumptions (B1), (B2) and (B3), by using a time domain transformation. System (11) provides a possibility to investigate the dynamics of the original system near the equilibrium points. We remark that the normal form (39) of the general type of system (11) obtained by using the method of multiple time scales is also applicable to the following type of system:

$ $$\begin{align}\label{tranformed-general-c-1} \left\{ \begin{aligned} \frac{1}{ \epsilon} \frac{\mathrm{d}r}{\mathrm{d} \eta}& = {-\mu_m r(\eta)+f(\xi(\eta-1))}\\ \frac{1}{ \epsilon} \frac{\mathrm{d}\xi}{\mathrm{d} \eta}& = {-\mu_p \xi(\eta)+g(r(\eta-1))}, \end{aligned} \right. \end{align}$$ $

(43)

which is equivalent to system (13). Even though many specific types of models of genetic regulatory dynamics have been developed in recent years (see [24], among many others), the general normal form we developed applies to a broad class of models.

Modeling the state-dependent delay of the genetic regulatory dynamics gives rise to the parameters $c>0$ and $\epsilon>0$ in model (11) with constant delay, where $c$ takes account of the state-dependent fluctuation of the diffusion time for homogenization of the substances produced in the regulatory network and $\epsilon$ measures the basal diffusion time. Since $c$ is the coefficient of the concentration disparity $x(t)-x(t-\tau)$, it is an analogue of the diffusion coefficients in [2]. The basal diffusion time $\epsilon$ is an analogue of the active macromolecular transport modeled as a time delay in [2].

From the analysis of model (11), we find that the characteristic equation is independent of $c$ while the basal diffusion time $\epsilon$ determines (with other parameters fixed) the stability of the equilibrium state and the approximate period of the bifurcating periodic solutions. However, the state-dependent effect represented by $c$ may give rise to both supercritical and subcritical Hopf bifurcations, while the former has been observed in [24], the latter is rarely seen in the literature. Note that the work of [24] is a theoretic analysis with normal form computation for the model with constant delay considered in [18], where the nonlinear feedbacks are Holling type Ⅲ functions and only numerical analysis was conducted.

A subcritical Hopf bifurcation usually means a sudden change of the dynamics when the system state is far enough from the stable equilibrium state. In terms of the model for the state-dependent delay $ \tau(t) = \epsilon+c(x(t)-x(t-\tau(t))), $ the subcritical Hopf bifurcation occurs when the contribution of homogenization to the diffusion time is strong enough and the basal diffusion time $\epsilon$ is less than its critical value. This observation may shed lights on the mutation phenomenon in certain genetic regulatory dynamics.

We remark that the numerical results in Section 4 is consistent with the observations from [2] described in the Section 1: when the analogous diffusion coefficient $c$ is large enough, the Hes1 system undergoes subcritical Hopf bifurcation with respect to $\epsilon$ and cannot sustain stable periodic oscillations; when the analogous diffusion coefficient $c$ is small, the Hes1 system undergoes supercritical Hopf bifurcation with respect to $\epsilon$ and can sustain stable periodic oscillations. The state-dependent delay provides an approach to take into account of both the diffusive and active macromolecular transports without the complication of involving spatial information which typically will lead to reaction-diffusion equations.

We also remark that the attempt to model the state-dependent delay with an inhomogeneous linear function of the state variables reveals that the concentration gradients in the genetic regulatory dynamics may change the Hopf bifurcation direction modeled with constant delays. Since diffusion of substances in live cell may be confined or otherwise free and is a complex process [25], improved models for the state-dependent delay in the diffusion process may provide better means for understanding the underlying mechanisms of genetic regulatory dynamics. For instance, what if the stationary state of $\tau$ depends on $c$ so that $c$ also changes the stability of the equilibria? It also remains to be investigated that to which extent the model of regulatory dynamics with linear type of state-dependent delay reflects that with a nonlinear one.

---

Acknowledgment

The author would like to thank an anonymous referee for the detailed and constructive comments which greatly improved the paper.

---