An Effective Algorithm for the Stability and Bifurcation in a DDE Model of Gene Expression

Abstract

The stability and Hopf bifurcation of gene expression models with a mechanism of delayed state feedback are considered. An effective algorithm for the calculations on the delay stable interval of the equilibrium point, the direction, and stability of the bifurcating periodic solution is also proposed. The $\tau$-decomposition strategy is applied to tackle the issue of local stability, and the explicit formula for the delay stable interval is provided. In addition, the asymptotical behaviors of the bifurcation solutions are investigated by the center manifold theorem and normal form theory. The direction and stability of the Hopf bifurcation are determined naturally. In addition, a subtle bilinear form of the adjoint system is proposed to calculate the bifurcation parameters directly. Finally, the correctness and effectiveness of our results and algorithm are verified by typical numerical examples.

1. Introduction

Gene expression refers to the interaction and restriction between genes, which involves the growth and development rules, morphological, and structural characteristics, and biological functions of plants and animals. If people master the mechanism of gene regulation and understand the concept of gene expression, it means they have mastered a key to unlock the secrets of biology. During the process of gene expression, a series of complex interactions occur between genes and proteins, which inevitably involve time delays between transcription and translation, resulting in the equilibrium and oscillation of the system. During this process, genes produce proteins, which in turn inhibit gene expression. We refer to this feedback as gene expression [1].

The model studied in this paper is similar to the literature [2]. As the first step in gene expression, transcription is catalyzed and regulated by DNA-dependent RNA polymerases whose function is to transfer genetic information from DNA strands to RNA strands. Numerous proteins called transcription factors are involved in the regulation of transcription. The mRNA generated by transcription acts as a blueprint for protein synthesis during translation. Translation is another and final step in gene expression that leads to protein synthesis. The protein then diffuses back into the nucleus, inhibiting transcription. Specifically, cells can regulate the levels of mRNA and protein concentrations by turning on and off specific gene transcription as a process, also known as feedback repression, in which gene expression is regulated by its own protein products. This feedback mechanism occurs when the protein returns to the nucleus to stop the transcription of its own mRNA by binding to the promoter site of the gene. Previous findings have shown that there are time delays in this feedback mechanism. These delays naturally occur as transcriptional delays (the time required for a gene to be copied into mRNA) and translation delays (the time required for the ribosome to translate mRNA into protein). Furthermore, it is sufficient for this article to consider only transcriptional time delays. Among them, the diagram of gene transcriptional translation is given in the first section of [2], making it more visual.

How to determine time delay in gene expression and their mechanisms in gene regulatory networks, and how to better regulate them, has naturally become a hot topic for many molecular biologists, systems biologists, and computational biologists. Goodwin [3] first studied the negative feedback gene expression model and simulated the differential equations governing enzyme synthesis by feedback inhibition, and the results showed the existence of nonlinear oscillations. Mackey [4] linked the occurrence of diseases with bifurcation dynamics and believed that there was a large class of dynamic diseases, which could control the operation of the system normally in the physiologic boundary area producing pathological behavior. It proved that there were periodic and non-periodic dynamics in the mathematical model of the physiologic system. Hori and Kim [5] investigated the conditions for the existence of periodic oscillations in large-scale cyclic gene regulatory networks, proved the local instability of the equilibrium point using the Poincaré–Bendixson theorem for cyclic systems, and proposed that systems can be composed of any number of genes in a large-scale cyclic gene regulatory network. John [6] considered a protein that repressed the transcription of its genes, noting that time delays cause negative feedback control to repeatedly overshoot and drop to steady state. We were inspired by the literature on the center manifold theorem. Wei and Yu [2] investigated the dynamics of a gene expression model with time delay and established the global existence of periodical solutions by using the center manifold theorem. Cao and Jiang [7] selected the transcription rate as the bifurcation parameter and performed nonlinear analysis by the center manifold theorem to verify the stability of the system at the critical value of the parameter. Wang and Yang [8] considered only time delay as a bifurcation parameter to study the gene regulatory network model, and determined explicit formulas for Hopf bifurcation direction and stability using the center manifold theorem. It was observed that noise plays a role in stabilizing an otherwise unstable oscillatory system. Djilali and Bentout [9] analyzed the behavior of the diffusive predator prey model using the regular form and establish the stability of the flush-periodic and non flush periodic solutions generated by the Hopf bifurcation. Soufiane and Touaoula [10] give sufficient conditions for the global asymptotic stability of free equilibrium related to the fundamental reproduction number. In addition, the global asymptotic stability in the presence of local equilibrium points is proved by using the Liapunov generalized function. Rand [11] gave an explicit expression for the radius of the limit cycle arising from Hopf bifurcation in a class of first-order differential delay equations with constant coefficients by Lindstedt’s method. Verdugo and Rand [1] applied the first-order nonlinear Lindstedt’s method he had previously studied in the model of gene transcription and protein synthesis and obtained a closed approximate expression of vibration amplitude and frequency. Alfifi [12] studied the solution of the one-dimensional reaction-diffusion equation by the Galerkin technique and investigated the effect of the free parameters in this model, and the results showed that they can destabilize or stabilize the solution of the equation. In addition, Lindstedt’s method in the perturbation theory is applied in the study of bifurcation. Verdugo [13] analyzed a delay differential equation with negative feedback, obtained amplitude and frequency using the multiple scale method of nonlinear analysis, and studied the importance of a balanced ratio between the synthesis rate and degradation rate in the presence of the periodical solution. Das and Chatterjee [14] studied the small perturbation problem of time–delay differential equations approximating Hopf bifurcation points and demonstrated the effectiveness of multi-scale methods. Li and Liu [15] proposed that time delay can change the level of Hes1 protein concentration, and pointed out that high steady-state Hes1 is susceptible to time delay, which excites oscillatory behavior. Shih and Yang [16] applied a factorization method to the characteristic polynomial of the bifurcation point of the negative feedback model to obtain the frequency of Hopf bifurcation, and then determined the stability of the periodic solution of Hopf bifurcation. Eva and Ileana [17] took the fractional order as the bifurcation parameter and proposed that Hopf bifurcation occurs only when the number of suppressor genes is odd based on the stability analysis. For further understanding, please refer to [18,19,20,21,22,23,24,25,26,27,28,29].

There are basically two types of genetic network models, namely Boolean (or discrete) models and differential equation models (or continuous models). In this paper, the differential equation model is chosen to describe the concentration of gene products in terms of variables more accurately and to provide more detailed information on the behavior of nonlinear dynamics.

Inspired by the above work, this paper studies a gene expression model:

$$\begin{array}{c}\dot{M}\left(t\right)=\frac{k}{1+{(P(t-\tau )/p_0)}^n}-\mu M\left(t\right)\hfill \\ \dot{P}\left(t\right)=M\left(t\right)-dP\left(t\right)\hfill \end{array}$$

(1)

where $M\left(t\right)$ and $P\left(t\right)$ represent the density of mRNA and protein, k is the initiation rate of transcription in the absence of the associated protein, $\mu$ and d denote the degradation rates of mRNA and protein, respectively, $\tau$ is the node delay caused by transcription, and n is called Hill coefficients.

The model is derived from [23] and can be traced back to [30,31]. The time dependence of the variables is shown explicitly in [31], but the delay differential equation in the text is not solved analytically, but only numerically. It is only based on some numerical simulations and does not give a rigorous mathematical proof of the Hopf bifurcation, and the direction of the Hopf bifurcation is not discussed in the paper. In addition, the effect of transcriptional parameters on the system was not considered in [23].

In this paper, we discuss the effect of transcription parameters for time delay, which we use as a static delay state feedback control gain. We adopt the $\tau$-decomposition method, which can quickly and accurately solve the equilibrium point and critical delay of the system. The results show that the critical value of delay is related to the static state delay feedback control. In the nonlinear part of the system, the center manifold theorem is combined with this paper, different from the nonlinear methods in other papers, such as Lindstedt’s method, perturbation method, or multiple scale method; these methods are more complicated and require a great deal of calculation. This paper encapsulates and simplifies the formulation based on the center manifold theorem. It can make the analysis of the nonlinear part clearer and more concise.

The remainder of this article is organized as follows: In Section 2, we make a preliminary analysis of the system and give a brief description of the $\tau$-decomposition method. In Section 3, we have two conclusions, one is the linear part, and we use $\tau$-decomposition to analyze the distribution of the characteristic equations associated with this model and obtain the existence of local Hopf bifurcations. In addition, the other is the nonlinear part, which uses the center manifold theorem to determine the direction of the Hopf bifurcation and the stability of the bifurcation period. In Section 4 section, numerical examples and simulation plots are given to confirm our theoretical results. Finally, in Section 5, we discuss our findings.

2. Problems Statement and Preliminaries

In this section, the main focus is on the introduction of the basic theory of the equations. The main goal of this section is to give a friendly introduction to the mathematical theory behind the rest of the paper. In this section, we will focus on the $\tau$-decomposition method. Before introducing the method, we also focus on giving an explanation of the linearization part of the system.

2.1. Initial Knowledge Preparation

It is not difficult to prove that the system has a unique equilibrium point $(M_e,P_e)$. The equilibrium point is found by making both equations of the system equal to zero. The calculation process uses the zero point theorem, derivatives, etc. The specific procedure is not described in this section. Introducing two new variables $m\left(t\right)=M\left(t\right)-M_e$ and $p\left(t\right)=P\left(t\right)-P_e$, system (1) becomes

$$\begin{array}{c}\dot{m}\left(t\right)=\frac{k}{1+{((p(t-\tau )+P_e)/p_0)}^n}-\mu (m\left(t\right)+M_e)\hfill \\ \dot{p}\left(t\right)=m\left(t\right)+M_e-d(p\left(t\right)+P_e)\hfill \end{array}$$

(2)

It can be seen that Equation (2) is modeled based on Equation (1) by moving the equilibrium point of the system to the origin.

In general, it is very important to consider the linear part of the system, and we study the local stability of the system by the linear equation of the system and consider the Hopf bifurcation on this basis. Therefore, for the study of the linear part of the system, it is necessary to turn system (2) into system (3), and we have the following matrix to represent the linear part of the system:

$$\begin{array}{c}\dot{m}\left(t\right)=-\mu m\left(t\right)-Kp(t-\tau )\hfill \\ \dot{p}\left(t\right)=m\left(t\right)-dp\left(t\right)\hfill \end{array}$$

(3)

within the static delayed state feedback control gain $K=\frac{kn\beta }{P_e{(1+\beta )}^2}$, and $\beta ={(P_e/p_0)}^n$.

The system (3) can be written as the following matrix equation:

$$\dot{u}\left(t\right)=Au\left(t\right)+Bu(t-\tau )$$

(4)

within

$$A=\left[\begin{array}{cc}-\mu & 0\\ 1& -d\end{array}\right],B=\left[\begin{array}{cc}0& -K\\ 0& 0\end{array}\right]\mathrm{and}u\left(t\right)=\left[\begin{array}{c}m\left(t\right)\\ p\left(t\right)\end{array}\right]$$

and the corresponding characteristic polynomial of Equation (4) is

$$P\left(\lambda \right)=det(\lambda I-A-B{e}^{-\tau \lambda })$$

(5)

Substituting the values of the matrix into Equation (5), we can obtain the characteristic polynomial of the system, which is

$$P\left(\lambda \right)={\lambda }^2+(\mu +d)\lambda +\mu d+K{e}^{-\tau \lambda }$$

(6)

A sufficient and necessary condition for the stability of the above linear system is that all eigenvalues lie on the left half-plane of the complex plane. In terms of the calculation of the stability interval of the system with delay, we consider the $\tau$-decomposition strategy to be an efficient computational method.

2.2. On the $\tau$-Decomposition Strategy

For the analysis of local asymptotic stability of time delay systems, the $\tau$-decomposition strategy is a better choice. This strategy takes the time delay as the variable parameter, first divides the delay $\tau$-axis into intervals with essentially the same stability characteristics within each interval, and then determines the stability switch by the crossover behavior at the endpoints of the interval.

The $\tau$-decomposition strategy has two main points. The first is to compute the pure imaginary roots (PIR) corresponding to the characteristic polynomial.

The characteristic equation can be simply written as follows by introducing polynomials $A(·)$ and $C(·)$ to represent the characteristic polynomial

$$P\left(\lambda \right)=A\left(\lambda \right)+C\left(\lambda \right)e^{-\tau \lambda }=0$$

(7)

Let $\lambda =\pm \omega j$ be substituted in Equation (7), and we have the following two equations:

$$\left\{\begin{array}{c}A\left(\omega j\right)+C\left(\omega j\right)e^{-\tau \omega j}=0\hfill \\ A(-\omega j)+C(-\omega j)e^{\tau \omega j}=0\hfill \end{array}\right.$$

deriving the formula about $e^{\tau \omega j}$

$$e^{\tau \omega j}=-\frac{C(-\omega j)}{A(-\omega j),}$$

then

$${\left|A\right(\omega j\left)\right|}^2={\left|C\right(\omega j\left)\right|}^2$$

(8)

where Equation (8) indicates that the two modal values are equal.

Another purpose is to determine whether the intersection of endpoints is stable or not. This is related to the direction of the crossover on each PIR. Assume that the positive roots ${\omega }_1,{\omega }_2,\dots ,{\omega }_p$ of $P\left(\lambda \right)=0$ are distinct and let ${\omega }_1>{\omega }_2>\dots >{\omega }_p>0$. Then, the crossing direction at the largest root ${\omega }_1$ is always to the right; at ${\omega }_2$, it is always to the left, and so on. If there is only one positive root ${\omega }_1$, we see that all crossings are to the right. If there are two positive roots ${\omega }_1$ and ${\omega }_2$, the crossings are to the right at ${\omega }_1$ but to the left at ${\omega }_2$ [18].

3. Main Results

In this section, we study two aspects, which are local stability analysis and Hopf bifurcation. In one of them, in the study of local stability, we obtain the values of specific critical time lags. In addition, in the study of Hopf bifurcation, we focus on the direction and type as well as the stability of the solution.

3.1. Local Stability and Hopf Bifurcation

In this subsection, we focus on the stability of the system at the equilibrium point and the conditions required for a Hopf bifurcation to occur. The discussion of local stability first requires us to obtain the value of the critical time delay. As described above in the first work on the $\tau$-decomposition method, we need to compute the pure imaginary roots. For the pure imaginary roots of the system, we give the following lemma.

Lemma 1.

Assuming that $K>\mu d$, the transcendental equation has only one purely imaginary root (PIR), where

$${\omega }_0={\left(\frac{1}{2}(-{\mu }^2-d^2+\sqrt{{({\mu }^2-d^2)}^2+4K^2})\right)}^{\frac{1}{2}}$$

(9)

Proof of Lemma 1.

The characteristic polynomial of the system (3) is

$$P\left(\lambda \right)={\lambda }^2+(\mu +d)\lambda +\mu d+K{e}^{-\tau \lambda }$$

Introducing A and C concerning $\lambda$ to represent the characteristic polynomial,

$$A\left(\lambda \right)={\lambda }^2+(\mu +d)\lambda +\mu d,C\left(\lambda \right)=K$$

Since conjugated complex roots appear symmetrically, we only need to consider positive pure imaginary roots, set $\lambda =\omega j$, and we can obtain

$$A\left(\omega j\right)=-{\omega }^2+(\mu +d)\omega j+\mu d,C\left(\omega j\right)=K$$

We use the method to obtain the following equation:

$${\omega }_0^4+{\mu }^2{\omega }_0^2+d^2{\omega }_0^2+{\mu }^2{d}^2-K^2=0$$

and the solution of this equation is

$${\omega }_0={\left(\frac{1}{2}(-{\mu }^2-d^2+\sqrt{{({\mu }^2-d^2)}^2+4K^2})\right)}^{\frac{1}{2}}$$

Thus, we obtain the result in (9).  □

On the other hand, the cross direction around ${\omega }_{0}j$ is always to the right, which is explicitly mentioned in the exposition of the $\tau$-decomposition method. Therefore, as for the local stability of system in (3), it is stated below.

Theorem 1.

The equilibrium point $(M_e,P_e)$ of system (3) is exponentially stable for $\tau \in [0,{\tau }_0)$, and ${\tau }_0$ is a Hopf bifurcation point, within

$${\tau }_0=\frac{{\theta }_0}{{\omega }_0}$$

(10)

in which ${\theta }_0=\angle ({\omega }_0^2+(\mu +d){\omega }_{0}j-\mu d)$, ∠ is arg (argument of a complex number). Moreover, ${\tau }_0$ is a Hopf bifurcation point for the system in (3).

Proof of Theorem 1.

Lemma 1 tells us that there is only one pair of $j{\omega }_0$ in the characteristic polynomial, which is the root of the characteristic equation, and then we have

$$-{\omega }_0^2+(\mu +d){\omega }_{0}j+\mu d+K{e}^{-{\tau }_0{\omega }_{0}j}=0$$

We can obtain an expression for $e^{{\tau }_0{\omega }_{0}j}$

$$e^{{\tau }_0{\omega }_{0}j}=\frac{K}{{\omega }_0^2-(\mu +d){\omega }_{0}j-\mu d}$$

Furthermore, the following formula can be obtained:

$${\tau }_0{\omega }_0=-\angle ({\omega }_0^2-(\mu +d){\omega }_{0}j-\mu d)$$

Therefore,

$${\tau }_0=\frac{\angle ({\omega }_0^2+(\mu +d){\omega }_{0}j-\mu d)}{{\omega }_0}$$

Thus, we obtain the result in (10). On the other hand, we will explain the reason why ${\tau }_0$ is a Hopf bifurcation point.

At $\tau ={\tau }_0$, it only corresponds to a pair of pure imaginary roots $\pm {\omega }_{0}j$. In addition, when $\tau \to {\tau }_0^{-}$, all roots lie on the left half plane, but when $\tau ={\tau }_0$, the root lies on the imaginary axis, at this time, all other roots lie on the left half plane, that is, Re$\left({\tau }_0\right)<0$. Moreover, since there is only one pair of the pure imaginary root, the crossing direction must be to the right, so we think ${\alpha }^{\prime }\left(0\right)>0$.  □

3.2. Direction and Stability of Bifurcating Periodic Solutions

In this section, we will consider the bifurcation solution of the system (1). In addition, we use the time delay $\tau$ as a bifurcation parameter. The system Equation (1) can be rewritten by Taylor expansion as follows:

$$\begin{array}{c}\begin{array}{c}\hfill \dot{m}\left(t\right)=-Kp(t-\tau )-\mu m\left(t\right)+h_2{p}^2(t-\tau )+h_3{p}^3(t-\tau )\end{array}\hfill \\ \dot{p}\left(t\right)=m\left(t\right)-dp\left(t\right)\hfill \end{array}$$

(11)

We could write system (11) as the matrix form

$$\dot{u}\left(t\right)=Au\left(t\right)+Bu(t-\tau )+F$$

(12)

within,

$$F=\left[\begin{array}{c}{f}_1\hfill \\ f_2\hfill \end{array}\right]=\left[\begin{array}{c}{h}_2{p}^2(t-\tau )+h_3{p}^3(t-\tau )\hfill \\ 0\hfill \end{array}\right]$$

within $h_2=\frac{\beta n(\beta n-n+\beta +1)}{2{(\beta +1)}^3{p}_e^2}$, $h_3=\frac{\beta n[6\beta n(1+\beta )(n-1)-(n-1)(n-2){(1+\beta )}^2-6{\left(\beta n\right)}^2]}{6{(\beta +1)}^4{p}_e^3}.$

In order to apply Hopf bifurcation theorem, it is necessary to turn Equation (12) to abstract differential equations, and the notation $u_t\left(\theta \right)$ is introduced:

$$u_t\left(\theta \right)=u(t+\theta )=\left[\begin{array}{c}m(t+\theta )\hfill \\ p(t+\theta )\hfill \end{array}\right],-\tau \le \theta \le 0$$

(13)

Equation (13) can induce the differential operator:

$${\displaystyle L{u}_t\left(\theta \right)=\left\{\begin{array}{c}\frac{d{u}_t\left(\theta \right)}{d\theta },-\tau \le \theta <0\\ \\ {\int }_{-\tau }^{0}d\eta \left(\theta \right)u_t\left(\theta \right),\theta =0\end{array}\right.}$$

(14)

where L is a continuous linear function mapping C into $R^2$, and $\eta$ is an 2 × 2 matrix function of bounded variation

$$d\eta \left(\theta \right)=\left[A\delta \right(\theta )+B\delta (\theta +\tau \left)\right]d\theta$$

and, within $\delta (·)$, is the Dirac delta function.

It should be noted that, in the form of the following inner product [19], we can derive the conjugate linear operator of this linear operator

$$<{\mu }^{*},\upsilon >={\overline{{\mu }^{*}}}^T\left(0\right)\upsilon \left(0\right)+{\int }_{-\tau }^0{\overline{{\mu }^{*}}}^T(\xi +\tau )B\upsilon \left(\xi \right)$$

(15)

and the conjugate linear operator is shown below

$${\displaystyle L^{*}{u}_t^{*}\left(\theta \right)=\left\{\begin{array}{c}-\frac{d{u}_t^{*}}{d\theta },-\tau \le \theta <0\\ \\ {\int }_{-\tau }^{0}d{\eta }^{\mathrm{T}}\left(\theta \right)u_t^{*}(-\theta ),\theta =0\end{array}\right.}$$

(16)

with

$$d{\eta }^{\mathrm{T}}\left(\theta \right)=[A^{\mathrm{T}}\delta \left(\theta \right)+B^{\mathrm{T}}\delta (\theta +\tau )]d\theta$$

Therefore, Equation (14) serves as the boundary value condition of the differential operator, and the following abstract differential equation can be derived

$$\dot{u}\left(t\right)=L{u}_t+F{u}_t$$

(17)

within, the nonlinear part in Equation (17) goes

$$F\left(u_t\right)=\left\{\begin{array}{c}\left[\begin{array}{c}0\\ 0\end{array}\right],-\tau \le \theta <0\\ \left[\begin{array}{c}{f}_w\\ 0\end{array}\right],\theta =0\end{array}\right.$$

(18)

with $f_w=h_2{p}^2(t-\tau )+h_3{p}^3(t-\tau )$ in Equation (18).

According to the theory of eigendecomposition, what we need to calculate is the eigenvector $q\left(\theta \right)$ corresponding to $j{\omega }_0$ of the linear operator L. We give the following lemma.

Lemma 2.

$q\left(\theta \right)$ is encoded in Equation (14) as the eigenvector of L, associated with $j{\omega }_0$ and satisfying

$$q\left(\theta \right)=\left[\begin{array}{c}{q}_1\\ 1\end{array}\right]e^{j{\omega }_0\theta }$$

(19)

within $q_1=j{\omega }_0+d$

Proof of Lemma 2.

From Equation (14), we can see that, when the range of $\theta$ is $[-\tau ,0)$, the action of the linear operator L can be colloquially understood as taking the derivative. In addition, in the linear part of the system analyzed earlier, the differential operator corresponds to the eigenvalue of $\pm j{\omega }_0$, we denote the associating eigenvector L with $j{\omega }_0$, and we can obtain the following formula:

$$q\left(\theta \right)=q\left(0\right)e^{j{\omega }_0\theta }$$

and, when $\theta =0$, the action of the linear operator L is equivalent to the differentiation, from which we can derive the following characteristic equation:

$$Lq\left(\theta \right)=j{\omega }_{0}q\left(\theta \right)=j{\omega }_{0}q\left(0\right)e^{j{\omega }_0\theta }$$

The above formula is expanded as

$$Aq\left(0\right)+Bq\left(0\right)e^{-j{\omega }_0{\tau }_0}=j{\omega }_{0}q\left(0\right)$$

thus

$$-\mu q_1-K{q}_2{e}^{-j{\omega }_0{\tau }_0}=j{\omega }_0{q}_1$$

$$q_1-d{q}_2=j{\omega }_0{q}_2$$

Set $q_2=1$, then $q_1=j{\omega }_0+d$,

$$q\left(\theta \right)=q\left(0\right)e^{j{\omega }_0\theta }=\left[\begin{array}{c}{q}_1\\ 1\end{array}\right]e^{j{\omega }_0\theta },q_1=j{\omega }_0+d$$

Thus, we obtain the result in (19).  □

In addition, corresponding to the above discussions, there is the following statement for conjugate eigenvectors.

Lemma 3.

$q^{*}\left(\theta \right)$ is encoded as an adjoint eigenvector of $L^{*}$ in Equation (16) associating with $-j{\omega }_0$ and such that

$$q^{*}\left(\theta \right)=\overline{D}\left[\begin{array}{c}1\\ q_2^{*}\end{array}\right]e^{j{\omega }_0\theta }$$

(20)

within $q_2^{*}=-j{\omega }_0+\mu$, $D=\frac{1}{\overline{<q,q_1^{*}>}}=\frac{1}{\mu +d-2j{\omega }_0-\tau K{e}^{j{\omega }_0{\tau }_0}}$

Proof of Lemma 3.

From Equation (16), we can observe a statement similar to the above, that is, when the range of $\theta$ is $[-\tau ,0)$, the action of the conjugate linear operator $L^{*}$ can be colloquially understood as taking the derivative, and $-j{\omega }_0$ is an eigenvalue of $L^{*}$, we can obtain the following equation:

$$q_1^{*}\left(\theta \right)=q_1^{*}\left(0\right)e^{-j{\omega }_0\theta }$$

and, when $\theta =0$, the action of the conjugate linear operator $L^{*}$ is equivalent to the differentiation, from which we can derive the following characteristic equation:

$$L^{*}{q}_1^{*}\left(\theta \right)=-j{\omega }_0{q}_1^{*}\left(\theta \right)=-j{\omega }_0{q}_1^{*}\left(0\right)e^{j{\omega }_0\theta }$$

Expanding the above formula, we can obtain the adjoint matrix equation

$$A^{\mathrm{T}}{q}_1^{*}\left(0\right)+B^{\mathrm{T}}{q}_1^{*}\left(0\right)e^{j{\omega }_0{\tau }_0}=-j{\omega }_0{q}_1^{*}\left(0\right)$$

thus

$$-\mu q_{01}^{*}+q_2^{*}=-j{\omega }_0{q}_{01}^{*}$$

$$-d{q}_2^{*}-K{q}_{01}^{*}{e}^{j{\omega }_0{\tau }_0}=-j{\omega }_0{q}_2^{*}$$

Setting $q_{01}^{*}=1$, we have $q_2^{*}=-j{\omega }_0+\mu$,

$$q_1^{*}\left(\theta \right)=q_1^{*}\left(0\right)e^{-j{\omega }_0\theta }=\left[\begin{array}{c}1\hfill \\ q_2^{*}\hfill \end{array}\right]e^{-j{\omega }_0\theta }$$

In this case, the $q_1^{*}\left(\theta \right)$ is not the normal orthogonal basis; in order to construct manifold near the origin $C\left(\mu \right)$ on the manifold $C\left(0\right),\mu =0$ coordinate system, we need to normalize it by introducing the inner product form references Equation (15), and the normalized eigenvector $q^{*}\left(\theta \right)=D{q}_1^{*}\left(\theta \right)$ should be such that $<q^{*},q>=1,$ and we have

$$<q,q_1^{*}>={\overline{q_1^{*}}}^T\left(0\right)q\left(0\right)+e^{j{\omega }_0{\tau }_0}{\overline{q_1^{*}}}^T\left(0\right)Bq\left(0\right)=\frac{1}{\overline{D}}$$

thus

$$D=\frac{1}{\overline{<q,q_1^{*}>}}=\frac{1}{\mu +d-2j{\omega }_0-\tau K{e}^{j{\omega }_0{\tau }_0}}$$

Thus, we obtain the result in (20).  □

In the above procedure, we give the adjoint form and the adjoint equation of the solution operator to the similarity transformation, and $q\left(\theta \right)$ and $q^{*}\left(\theta \right)$ can be obtained. Then, we can use the conjugate eigenvector $q^{*}\left(\theta \right)$ to project the central epidemic. After determining the accompanying vectors, the solution space can be divided into central and stable manifolds Like in literature [7], we define

$$z=<q^{*},u_t>$$

(21)

and

$$y(t,\theta )=u_t\left(\theta \right)-z\left(t\right)q\left(\theta \right)-\overline{z}\left(t\right)\overline{q}\left(\theta \right)$$

(22)

The solution of the system is composed of popular linear combinations, z is the projected solution at time t, and y is a high-order term.

Lemma 4.

The evolution of z concerning time t satisfies the following form:

$$\dot{z}\left(t\right)=j{\omega }_{0}z+g_{20}{z}^2+g_{11}z\overline{z}+g_{02}{\overline{z}}^2+g_{21}{z}^2\overline{z}+\cdots$$

(23)

in which

$$g_{20}=\overline{D}{f}_{20},g_{11}=\overline{D}{f}_{11},g_{02}=\overline{D}{f}_{02},and{g}_{21}=\overline{D}{f}_{21}$$

(24)

within

$$\begin{array}{c}{f}_{20}=2h_2{e}^{-2j{\omega }_0{\tau }_0}\hfill \\ f_{11}=2h_2\hfill \\ f_{02}=2h_2{e}^{2j{\omega }_0{\tau }_0}\hfill \\ f_{21}=2[2h_2{e}^{-j{\omega }_0{\tau }_0}{y}_{11}^{\left(2\right)}(-\tau )+h_2{e}^{j{\omega }_0{\tau }_0}{y}_{20}^{\left(2\right)}(-\tau )+3h_3{e}^{-2j{\omega }_0{\tau }_0}]\hfill \end{array}$$

Proof of Lemma 4.

Near the origin O, take q and $\overline{q}$ as local coordinate direction, and $z\left(t\right)$ and $\overline{z}\left(t\right)$ is the manifold manifold in the C. Based on Equations (21) and (23), we can derive the following equation:

$${\displaystyle \begin{array}{ccc}\dot{z}\left(t\right)\hfill & =& <q^{*},{\dot{u}}_t>=<q^{*},L{u}_t+F{u}_t>\hfill \\ \hfill & =& j{\omega }_{0}z\left(t\right)+{\overline{q^{*}}}^T\left(0\right)F{u}_t\left(0\right)\hfill \\ \hfill & =& j{\omega }_{0}z+\overline{D}\left[1\overline{q_2^{*}}\right]\left[\begin{array}{c}{f}_w\\ 0\end{array}\right]\hfill \\ \hfill & =& j{\omega }_{0}z+\overline{D}\left[1\overline{q_2^{*}}\right]\left[F_{20}{z}^2+F_{11}z\overline{z}+F_{02}{\overline{z}}^2+F_{21}{z}^2\overline{z}+\cdots \right]\hfill \end{array}}$$

in which

$$F_{20}=\left[\begin{array}{c}{f}_{20}\hfill \\ 0\hfill \end{array}\right],F_{11}=\left[\begin{array}{c}{f}_{11}\hfill \\ 0\hfill \end{array}\right]$$

$$F_{02}=\left[\begin{array}{c}{f}_{02}\hfill \\ 0\hfill \end{array}\right],F_{21}=\left[\begin{array}{c}{f}_{21}\hfill \\ 0\hfill \end{array}\right]$$

the center manifold theorem, and Equation (22) shows

$$y(t,\theta )=\frac{1}{2}{y}_{20}\left(\theta \right)z^2\left(t\right)+y_{11}\left(\theta \right)z\left(t\right)\overline{z}\left(t\right)+\cdots$$

By definition, we know

$$u_t\left(\theta \right)=z\left(t\right)q\left(\theta \right)+\overline{z}\left(t\right)\overline{q}\left(\theta \right)+\frac{1}{2}{y}_{20}\left(\theta \right)+y_{11}\left(\theta \right)z\overline{z}+\cdots$$

According to Equation (13),

$$u_t\left(\theta \right)=\left[\begin{array}{c}{m}_t\left(\theta \right)\\ p_t\left(\theta \right)\end{array}\right]=\left[\begin{array}{c}m(t+\theta )\\ p(t+\theta )\end{array}\right]$$

then, we can convert the system into the form of u

$$\begin{array}{ccc}m\left(t\right)\hfill & =& m(t+0)=u_t^{\left(1\right)}\left(0\right)\hfill \\ p\left(t\right)\hfill & =& p(t+0)=u_t^{\left(2\right)}\left(0\right)\hfill \end{array}$$

$$\begin{array}{ccc}m(t-\tau )\hfill & =& m(t+(-\tau ))=u_t^{\left(1\right)}(-\tau )\hfill \\ p(t-\tau )\hfill & =& p(t+(-\tau ))=u_t^{\left(2\right)}(-\tau )\hfill \end{array}$$

within

$$u_t^{\left(2\right)}(-\tau )=z{q}^{\left(2\right)}(-\tau )+\overline{z}{\overline{q}}^{\left(2\right)}(-\tau )+\frac{1}{2}{y}_{20}^{\left(2\right)}(-\tau )z^2+y_{11}^{\left(2\right)}(-\tau )z\overline{z}$$

For our system, the nonlinearities have been sorted out earlier, and then we can obtain

$$\begin{array}{cc}\hfill h_2{p}^2(t-\tau )=h_2{\left(u_t^2(-\tau )\right)}^2=& h_2{[z{q}^{\left(2\right)}(-\tau )+\overline{z}{\overline{q}}^{\left(2\right)}(-\tau )+\frac{1}{2}{y}_{20}^{\left(2\right)}(-\tau )z^2+y_{11}^{\left(2\right)}(-\tau )z\overline{z}]}^2\hfill \end{array}$$

$$\begin{array}{cc}\hfill h_3{p}^3(t-\tau )=h_3{\left(u_t^2(-\tau )\right)}^3=& h_3{[z{q}^{\left(2\right)}(-\tau )+\overline{z}{\overline{q}}^{\left(2\right)}(-\tau )+\frac{1}{2}{y}_{20}^{\left(2\right)}(-\tau )z^2+y_{11}^{\left(2\right)}(-\tau )z\overline{z}]}^3\hfill \end{array}$$

so we can figure out the following

$$\begin{array}{ccc}{f}_{20}\hfill & =& 2h_2{e}^{-2j{\omega }_0{\tau }_0}\hfill \\ f_{11}\hfill & =& 2h_2\hfill \\ f_{02}\hfill & =& 2h_2{e}^{2j{\omega }_0{\tau }_0}\hfill \end{array}$$

$$\begin{array}{c}\hfill f_{21}=2[2h_2{e}^{-j{\omega }_0{\tau }_0}{y}_{11}^{\left(2\right)}(-\tau )+h_2{e}^{j{\omega }_0{\tau }_0}{y}_{20}^{\left(2\right)}(-\tau )+3h_3{e}^{-2j{\omega }_0{\tau }_0}]\end{array}$$

Thus, we have

$$\begin{array}{ccc}{g}_{20}\hfill & =& {\overline{q^{*}}}^T\left(0\right)f_{20}=\overline{D}{f}_{20}\hfill \\ g_{11}\hfill & =& {\overline{q^{*}}}^T\left(0\right)f_{11}=\overline{D}{f}_{11}\hfill \\ g_{02}\hfill & =& {\overline{q^{*}}}^T\left(0\right)f_{02}=\overline{D}{f}_{02}\hfill \\ g_{21}\hfill & =& {\overline{q^{*}}}^T\left(0\right)f_{21}=\overline{D}{f}_{21}\hfill \end{array}$$

Thus, we obtain the result in (24).  □

It is easy to see that $g_{2}0,g_{1}1and{g}_{0}2$ can be determined directly from the system parameters. However, the result of $g_{2}1$ is closely related to the calculation of $y_{11},y_{20}$.

Therefore, we summarized the calculations on $y_{20}(-\tau )$ and $y_{11}(-\tau )$ in the following lemma.

Lemma 5

([19]). The coefficients $y_{20}^{\left(2\right)}(-\tau ),y_{11}^{\left(2\right)}(-\tau )$ of the center manifold $y\left(t\right)$ are such that

$$\begin{array}{c}\begin{array}{c}\hfill y_{20}^{\left(2\right)}(-\tau )=\frac{2j}{\omega }\overline{D}{h}_2{e}^{-4j\omega \tau }+\frac{2j}{3\omega }D{h}_2+C_{202}{e}^{-2j\omega \tau }\end{array}\hfill \end{array}$$

(25)

$$\begin{array}{c}\begin{array}{cc}\hfill y_{11}^{\left(2\right)}(-\tau )=& -\frac{2j}{\omega }\overline{D}{h}_2{e}^{-j\omega \tau }+\frac{2j}{\omega }D{h}_2{e}^{j\omega \tau }+C_{112},\hfill \end{array}\hfill \end{array}$$

(26)

in which

$$C_{20}={\left[\begin{array}{cc}2j{\omega }_0+{\tau }_0\mu \hfill & {\tau }_{0}K{e}^{-2j{\omega }_0{\tau }_0}\\ -{\tau }_0\hfill & 2j{\omega }_0+{\tau }_{0}d\end{array}\right]}^{-1}\left[\begin{array}{c}{f}_{20}\hfill \\ 0\hfill \end{array}\right]=\left[\begin{array}{c}{C}_{201}\hfill \\ C_{202}\hfill \end{array}\right]$$

$$C_{11}=-\frac{1}{{\tau }_0}{\left[\begin{array}{cc}-\mu \hfill & -K\\ 1\hfill & -d\end{array}\right]}^{-1}\left[\begin{array}{c}{f}_{11}\hfill \\ 0\hfill \end{array}\right]=\left[\begin{array}{c}{C}_{111}\hfill \\ C_{112}\hfill \end{array}\right].$$

According to Equations (25) and (26), up to this point, we can obtain all expressions of $g_{ij}$. Therefore, we have the following theorem about the nonlinear dynamics of the system.

For applying the Hassard formula,

$$\left\{\begin{array}{cc}\hfill {\mu }_2& =-\frac{\mathrm{Re}\left[C_1\left(0\right)\right]}{{\alpha }^{\prime }\left(0\right)]}\hfill \\ \hfill {\tau }_2& =-\frac{\mathrm{Im}\left[C_1\left(0\right)\right]+{\mu }_2{\omega }^{\prime }\left(0\right)}{{\omega }_0}\hfill \\ \hfill {\beta }_2& =2\mathrm{Re}\left[C_1\left(0\right)\right]\hfill \end{array}\right.$$

(27)

with

$$\begin{array}{cc}\hfill C_1\left(0\right)=& \frac{j}{2{\omega }_0}[4h_2^2{\overline{D}}^2{e}^{-2j{\omega }_0{\tau }_0}-4h_2{\left|\overline{D}\right|}^2-\frac{2}{3}{h}_2{\left|\overline{D}{e}^{2j{\omega }_0{\tau }_0}\right|}^2]\hfill \\ \hfill & +2h_2\overline{D}{e}^{-j{\omega }_0{\tau }_0}{y}_{11}^{\left(2\right)}(-\tau )+h_2{e}^{j{\omega }_0{\tau }_0}{y}_{20}^{\left(2\right)}(-\tau )+3h_3{e}^{-j{\omega }_0{\tau }_0}\hfill \end{array}$$

It is noticed that all the bifurcation parameters can be determined in a systematic way, with the value of

$$-\frac{\partial P/\partial \tau }{\partial P/\partial \lambda }\left|{}_{({\tau }_0,j{\omega }_0)}\right.={\alpha }^{\prime }\left(0\right)+{\omega }^{\prime }\left(0\right)j$$

Theorem 2.

In Formulas (27), ${\mu }_2$ determines the direction of the bifurcation; if $\mu >0$ (<0), the periodic solution is called supercritical (subcritical). ${\tau }_2$ determines the period of the periodic solution of the bifurcation; if ${\tau }_2$ > 0 (<0), then the period increases (decreases). ${\beta }_2$ determines the stability of the periodic solution. If (${\beta }_2$ < 0) (>0), the periodic solution is asymptotically stable (unstable).

4. Numerical Examples

In this section, we will apply a typical example from other literature to verify the correctness and validity of our results. All numerical calculations are performed by MATLAB. Moreover, the analysis shows that our technique is a reliable approach to solve the nonlinear dynamics of gene expression systems.

Example 1.

The typical example comes from [2], and it should be pointed out that this example is a special case in the study of this article

$$\begin{array}{ccc}\dot{m}\left(t\right)\hfill & =& \frac{k}{1+{(p(t-\tau )/40)}^4}-0.16m\hfill \\ \dot{p}\left(t\right)\hfill & =& m-0.1p\hfill \end{array}$$

(28)

$\mu =0.16$ and $d=0.1$ denote the degradation rates of mRNA and protein, respectively, where $p_0=40$ is a reference concentration of protein $n=4$.

The associating characteristic polynomial (6) is

$${\lambda }^2+0.26\lambda +0.016+K{e}^{-\tau \lambda }$$

In the following, we will discuss the special case when $k=1$. Lemma 1 tells us that

$${\omega }_0=0.0990$$

and Theorem (1) leads us to the critical delay

$${\tau }_0=18.2674$$

and

$$Re{\lambda }^{{}^{\prime }}\left({\tau }_0\right)=0.00092$$

Our simulation diagram gives the time delay stability interval, through the previous theoretical calculation, and the interval we calculated is consistent with the interval in the simulation, and both aspects can verify each other.

We take the root trajectory Figure 1, and we find that, at $\tau <18.2674$, the system is stable since the eigenroots all have negative real parts. In addition, a pair of pure imaginary roots crosses the imaginary axis at ${\tau }_0=18.2674$, resulting in the unstable oscillation of the system, that is, Hopf bifurcation. Thus, there exists only one delay interval $[0,18.2674]$, for the system of gene expression.

The corresponding parameters are calculated as follows:

Lemmas 2 and 3 give the adjoint vectors as

$$q\left(\theta \right)=\left[\begin{array}{c}0.1000+0.0990j\\ 1\end{array}\right]e^{j{\omega }_0\theta }$$

$$q^{*}\left(\theta \right)=\left[\begin{array}{c}1\\ 0.2418+0.1192j\end{array}\right]e^{j{\omega }_0\theta }$$

Associating with Lemma 4, we obtain

$$\begin{array}{c}{H}_2=0.000111,H_3=0.000024\hfill \\ g_{20}=-0.000011+0.00029j\hfill \\ g_{11}=0.000142-0.000254j\hfill \\ g_{02}=-0.000242+0.00016j\hfill \end{array}$$

Subsequently, Lemma 5 describes that $y_{20}$ and $y_{11}$ are such that

$$\begin{array}{ccc}{y}_{20}^{\left(2\right)}(-\tau )\hfill & =& 0.001752+0.004121j\hfill \\ y_{11}^{\left(2\right)}(-\tau )\hfill & =& -0.001292\hfill \end{array}$$

Thus, we have

$$g_{21}=-0.000182-0.000049j$$

Therefore, it gives the bifurcation parameters

$${\mu }_2=0.099129,{\tau }_2=0.003557and{\beta }_2=-0.000182$$

Theorem 2 tells us that the Hopf bifurcation occurs at $\tau =18.2674$, and the direction of bifurcation is supercritical.

In addition, we draw stable and unstable phase trajectory in Figure 2 and Figure 3.

Furthermore, we took a representative point inside the stability interval, $\tau =14$; this time delay represents that the system is stable. We drew a phase diagram according to this time lag, and the track of the phase diagram finally converged to the origin, which also shows that the system is stable. Then, it shows that we obtained this time lag, and it is indeed stable.

At the same time, we took another point from the unstable interval. For $\tau =21$ in Figure 4, we chose two different initial states. It is easy to see that both of the two phase trajectories are convergent to the bifurcating periodic solution. It means that both phase trajectories are divergent, and both are further and further away from the equilibrium point. In addition, this means that the equilibrium points are unstable.

Finally, we took the bifurcation point, a special point, and found that the graph converges to the limit ring. The bifurcation means that it starts from the initial state and finally converges more and more to the limit ring, and no longer converges to the equilibrium point. At the beginning, the system has the concept of equilibrium point, and as the time lag changes, it will bifurcate and the original point becomes a circle. As is shown in Theorem 1, the ${\tau }_0=18.2674$ is the Hopf bifurcation point. With $\tau =18.2674$, there is a periodic convergence phase trajectory.

5. Conclusions

In this paper, we propose a unified framework for studying single genetic negative feedback autoregulation systems with time delay. According to the $\tau$-decomposition strategy, the maximum value of the time delay stability bound is given. In addition, we consider the time delay as a bifurcation parameter and prove that Hopf bifurcation occurs when the delay passes this threshold. The direction and stability of the bifurcated periodic orbits are studied using the central manifold theorem, the bilinear form and the projection vector. In fact, the results we discuss can be directly applied to other fields related to two-dimensional delay differential equations. Moreover, a numerical example partially verifies the correctness and reliability of the obtained results.