Bifurcation analysis of a two–dimensional p53 gene regulatory network without and with time delay

1.   Introduction

Dynamical analysis of gene regulatory networks (GRNs) characterized by mathematical models plays an important role in understanding the underlying mechanism of the corresponding biological processes and predicting what will happen^[1,2,3]. A GRN is composed of genes, RNAs, and proteins represented by nodes, which are connected through the edges that indicate the transformation, promotion, and inhibition between them^[4]. The expression level of proteins in the GRN is closely related to the mechanisms of a number of biochemical processes, including cell differentiation^[5,6,7], neural plasticity^[8,9,10] and the development of cancer^[11,12,13]. Typically, cell fates in response to stresses are closely related to the dynamics of the tumor suppressor protein p53 in p53 GRN with core p53–Mdm2 feedback loops^[14,15]. p53 maintains a low level under a homeostatic condition while moderate and high stresses make p53 exhibit oscillation and high level, which lead to cell cycle arrest and apoptosis^[16,17], respectively. Therefore, more and more research has been focusing on exploring p53 dynamics numerically and theoretically through the construction of mathematical models with the help of bifurcation diagrams^[18,19].

Bifurcation diagrams of mathematical models is a useful tool for analyzing the dynamics of GRN ^[20]. Bifurcation diagrams of p53 GRN in ^[16,17,21] display various types of bifurcation, such as codimension–1 saddle–node, Hopf bifurcation and codimension–2 Bogdanov–Takens bifurcation. Saddle–node bifurcation can generate the bistability with a low and high expression levels of p53, which results in cell survival and cell apoptosis, respectively. Hopf bifurcation can induce the appearance of a stable limit cycle ^[22] corresponding to the oscillation expression of p53, which result in cell cycle arrest^[23]. Codimension–2 Bogdanov–Takens bifurcation may give rise to the coexistence of a stable steady state and a stable limit cycle^[24], which correspond to cell survival and cell cycle arrest, respectively. Therefore, the analysis of the conditions under which these bifurcations occur can allow for a deeper understanding of cell fate decision in response to different parameters. More research has focused on bifurcation analyses of high–dimensional p53 GRNs through numerical simulations^[25,26,27]. However, theoretical analysis of low-dimensional p53 GRNs contribute to the understanding of cell fate decisions under different conditions^[28,29].

Theoretical analyses of the bifurcation of dynamical systems described by ordinary differential equations play an important role in unveiling their complex dynamic properties. Previously, extensive effort has been devoted to bifurcation analyses of predator–prey models and SIR models of infectious diseases with various factors^[30], while several recent studies have theoretically revealed the bifurcation of GRNs. These studies focus on investigating the existence and stability of possible equilibria of the system and deriving the rigorous mathematical proofs for the existence of bifurcations, such as saddle-node bifurcation, Hopf bifurcation of codimension–1 and Bogdanov–Takes bifurcation of codimension–2 or codimension–3^{[31,32,33,34]}. Besides, Hopf bifurcation may be caused by time delay in GRNs, and it is analyzed by studying the associated characteristic equation of the corresponding linearized system^[18,35]. Although there has been much research on bifurcation analyses for p53 GRNs through numerical simulations, there is scant theoretical analyses of bifurcations of low dimensional p53 GRNs.

In this paper, we investigated the bifurcation of a two–dimensional p53 GRN without and with time delay, as described in ^[14], by performing rigorous mathematical analysis. Firstly, the existence of all possible positive equilibria are investigated by applying Descartes' rule of signs and the local stability of the positive equilibria are analyzed. Then, in the absence of time delay, the conditions for the appearance of codimension–1 saddle–node, Hopf bifurcation and codimension–2 Bogdanov–Takens bifurcation are derived by using Sotomayor's theorem^[36], Hopf bifurcation theorem^[37] and the normal form method, respectively, and the first Lyapunov number is calculated to obtain the stability of the limit cycle. Furthermore, in the presence of time delay, the Hopf bifurcation induced by time delay is analyzed on the basis of the Hopf bifurcation theorem. These theoretical results are numerically supported by bifurcation diagrams and phase portraits, and they can be considered as a complement to existing literature on the dynamics of p53 GRN.

The organization of this paper is as follows. The mathematical model of the p53 GRN is given in Section 2. The existence and local stability of positive equilibria of the p53 GRN without time delay are analyzed in Section 3. Section 4 presents the conditions for the occurrence of codimension–1 saddle–node bifurcation, Hopf bifurcation and codimension–2 Bogdanov–Takens bifurcation of the p53 GRN without time delay through theoretical analyses, which are supported by bifurcation diagrams and phase portraits. Section 5 is devoted to the analysis of Hopf bifurcation induced by time delay. We end the paper with the conclusion in Section 6.

2.   The model of the p53 GRN

In the present work, we consider a core p53 GRN with p53 and its key regulator Mdm2 in ^[14,15], as shown in Figure 1. Figure 1 includes a p53 self–induction positive feedback loop and a negative feedback loop, where p53 elevates the expression level of the Mdm2 protein and Mdm2 promotes the degradation of p53. Here, the degradation of p53 is regulated by the concentration of Mdm2 at some previous time. The rate equations for the concentration of p53 (denoted by $x$) and that of Mdm2 (denoted by $y$) are given by the following delay different equations:

Here, $r_1$ and $r_2$ denote the basal production rates of p53 and Mdm2, respectively. Also the production of both p53 and Mdm2 activated by p53 are modeled by using Hill functions with the production rates $v_1$ and $v_3$, and Michaelis constants $k_1$ and $k_3$, respectively. $d_1$ and $d_2$ are the basal degradation rates of p53 and Mdm2, respectively. Besides, p53 is degraded by Mdm2 at a rate $v_2$ in a Michaelis–Menten function with the Michaelis constant $k_2$. Time delays $\tau_1$ and $\tau_2$ characterize the periods of time for gene expression to protein production of p53 and Mdm2.

3.   Analysis of the positive equilibria

In this section, the existence and stability of the positive equilibria of the system $(2.1)$ are investigated and the qualitative behavior of the system $(2.1)$ are given in the following subsections.

3.1.   Existence of the positive equilibria

In the section, we focus on analyzing the conditions for the existence of a positive equilibrium in the system $(2.1)$ for biological reasons. Assuming a positive equilibrium of the system $(2.1)$ is $E(x_{*}, y_{*})$, which satisfies the following equations:

Obviously, the second equation of Eq $(3.1)$ is equivalent to

Rearranging the first equation of Eq $(3.1)$, we get

where

Obviously, $x_{*}$ is the root of the following equation

If the root $x_{*}$ of Eq $(3.4)$ is positive, $y_{*}$ is positive according to Eq $(3.2)$ with positive rate constants. Therefore, the conditions for the existence of positive roots of Eq $(3.4)$ are suitable for the one of the positive equilibria of the system $(2.1)$.

Applying Descartes' rule of signs to Eq $(3.4)$, we obtained the number of possible positive equilibria in the system $(2.1)$ is concluded in Table 1. Obviously, the system $(2.1)$ must have at least one positive equilibrium for $C_0 C_6 < 0$. Besides, according to Eq $(3.6)$, we conclude that if $C_2 < 0$, then $C_4 < 0$ and if $C_4 > 0$, then $C_2 > 0$. Hence, based on the cases 1–4 in Table 1, a unique positive equilibrium $E(x_{*}, y_{*})$ in the system $(2.1)$ exists under the conditions in the following theorem.

Theorem 3.1. The system $(2.1)$ has a unique positive equilibrium $E(x_{*}, y_{*})$ if one of the following conditions holds:

$(i)\; C_2 > 0, C_3 > 0, C_5 > 0$;

$(ii)\; C_1 > 0, C_3 > 0, C_4 > 0$;

$(iii)\; C_1 > 0, C_2 > 0, C_4 > 0, C_5 > 0$.

Theorem $3.1$ is verified by the nullclines of $x$ and $y$ in Figure 2, where the system $(2.1)$ has a unique positive equilibrium for the condition (i) $C_2 = -0.1197 < 0$, $C_3 = -0.1384 < 0$, $C_5 = -0.000723 < 0$ with $v_1 = 0.18$, $v_2 = 0.003$, $v_3 = 0.6$, $d_1 = 0.034$, $d_2 = 0.02$, $r_1 = 0.01$, $r_2 = 0.001$, $k_1 = 6$, $k_2 = 4$, $k_3 = 4$ in Figure 2 (a) and the condition (ii) $C_1 = 0.0057 > 0$, $C_3 = 0.1080 > 0$, $C_4 = 0.0024 > 0$ with $v_1 = 0.63$, $v_2 = 0.01$, $v_3 = 0.985$, $d_1 = 0.004$, $d_2 = 0.0025$, $r_1 = 0.021$, $r_2 = 0.0001$, $k_1 = 1.5$, $k_2 = 2$, $k_3 = 9$ in Figure 2 (b). Besides, Figure 2 (c) shows that the system $(2.1)$ has three positive equilibria for $C_1 = -0.1585 > 0$, $C_2 = 0.0932 > 0$, $C_3 = 0.0173 > 0$, $C_4 = -0.0107 > 0$, $C_5 = -0.0018 > 0$ in case 9 with $v_1 = 0.18$, $v_2 = 0.01$, $v_3 = 0.55$, $d_1 = 0.034$, $d_2 = 0.03$, $r_1 = 0.011$, $r_2 = 0.001$, $k_1 = 2.4$, $k_2 = 2$, $k_3 = 4$.

Here, the conditions for the existence of the positive equilibria in system $(2.1)$ are given in Table 1. Furthermore, the stability of the positive equilibria is analyzed in the next section.

3.2.   Stability of the positive equilibria

To investigate the stability of any positive equilibrium $E(x_{*}, y_{*})$ of the system $(2.1)$, the corresponding Jacobian matrix $J$ is given by

The characteristic equation of the system $(2.1)$ is

The trace $\mathrm{tr}(J)$ and the determinant $\mathrm{det}(J)$ of the Jacobian matrix $J$ are given by

The local stability of the positive equilibrium $E(x_{*}, y_{*})$ is decided by the signs of $\mathrm{tr}(J)$ and $\mathrm{det}(J)$. Next, we will first study the sign of $\mathrm{det}(J)$. Since

According to $F(x_{*}) = 0$ in Eq $(3.4)$, we conclude that

and the signs of $\mathrm{det}(J)$ and $F^{\prime}(x_{*})$ are opposite. Therefore, the signs of $F^{\prime}(x_{*})$ and $\mathrm{tr}(J)$ decide the stability of the positive equilibrium $E(x_{*}, y_{*})$, which are given in the following theorem.

Theorem 3.2. The stability of the positive equilibrium $E(x_{*}, y_{*})$ under different conditions is described in Table $2$.

However, the signs of $F^{\prime}(x_{*})$ and $\mathrm{tr}(J)$ are not decided explicitly due to their complex expression, so we will give some numerical examples to illustrate the stability of $E(x_{*}, y_{*})$ in the following section. The stability of $E(x_{*}, y_{*})$ can be changed by the bifurcation, which will be explored in the following bifurcation analysis.

4.   Bifurcation analysis of p53 GRE without time delay

Bifurcation changes the stability and the number of equilibria in the system as the parameter varies through a critical value. In this section, we investigate the conditions for the occurrence of codimension–1 saddle–node and Hopf bifurcation with respect to $v_3$ and codimension–2 Bogdanov–Takens bifurcation with respect to $v_3$ and $d_2$ in the system $(2.1)$ without time delay.

4.1.   Saddle-node bifurcation

According to Sotomayor theorem^[36], we shall establish the conditions under which the system $(2.1)$ experiences saddle–node bifurcation at the equilibrium $E(x_{*}, y_{*})$ when the control parameter $v_3$ crosses the critical value $v_3^{S N}$.

The first condition is that $\mathrm{tr}(J(E))|_{v_3 = v_3^{S N}} \neq 0$ and $\mathrm{det}(J(E))|_{v_3 = v_3^{S N}} = 0$, which correspond to

based on Eqs $(3.8)$ and $(3.10)$.

Next, to obtain the transversality conditions, the eigenvectors of the matrices $J(E, v_3^{S N})$ and $J^T(E, v_3^{S N})$ with the eigenvalue $\lambda = 0$ are given as follows, respectively,

Denoting $f(x, y) = (f_1(x, y), f_2(x, y))^T$,

Clearly, the transversality conditions are

Therefore, based on the above analyses, we get the following theorem.

Theorem 4.1. The system $(2.1)$ experiences saddle–node bifurcation at the positive equilibrium $E(x_{*}, y_{*})$ as the parameter $v_3$ crosses the critical value $v_3^{S N}$ if the conditions (SN.1) and (SN.2) hold. Here $v_3^{S N} = \frac{2 k_1^2 v_1d_2 x_{*}(k_2+x_{*})(k_3^2 + x_{*}^2)^2}{2 k_3^2 v_2 x_{*}^2(k_1^2 +x_{*}^2)^2}-\frac{k_2 v_2 d_2 y_{*}(k_3^2 + x_{*}^2)^2}{2 k_3^2 v_2 x_{*}^2(k_2 + x_{*})}-\frac{d_1d_2(k_2+x_{*})(k_3^2 + x_{*}^2)^2}{2 k_3^2 v_2 x_{*}^2}$.

Theorem 4.1 is verified by the bifurcation diagram of $x$ with respect to $v_3$ in Figure 3(a) and phase portraits of $x$ and $y$ for five typical values of $v_3$ in Figure 3(b)–(f) with the same parameters $v_1 = 0.18$, $v_2 = 0.01$, $d_1 = 0.034$, $d_2 = 0.03$, $r_1 = 0.011$, $r_2 = 0.001$, $k_1 = 2.4$, $k_2 = 2$, $k_3 = 4$. In Figure 3(a), black solid and dashed lines represent stable and unstable equilibria, respectively, which meet at two saddle–node bifurcation points $SN_1$ and $SN_2$ with $v_3^{SN1} = 0.4606397$ and $v_3^{SN2} = 0.5991365$, respectively. In Figure 3(b)–(f), the solid lines represent the trajectory running along the arrows and red solid and hollow dots denote the stable and unstable equilibria, respectively.

As shown in Figure 3(a), the system $(2.1)$ undergoes saddle–node bifurcation at $E_1(x_{*}, y_{*}) = (0.7593, 0.5674)$ and $E_2(x_{*}, y_{*}) = (1.3543, 2.0872)$ as $v_3$ passes through $v_3^{SN1} = 0.4606397$ and $v_3^{SN2} = 0.5991365$, respectively, with $\mathrm{tr}(J(E_1, v_3^{SN1})) = -0.0263\neq 0$, $W^T\left[D^2 f\left(E_1; v_3^{ S N1}\right)(V, V)\right] = 0.0118 \neq 0$ and $\mathrm{tr}(J(E_2, v_3^{SN2})) = -0.0315 \neq 0$, $W^T\left[D^2 f\left(E_2; v_3^{ S N2}\right)(V, V)\right] = 0.0104\neq 0$, which meet all conditions in Theorem {4.1}. $v_3^{SN1}$ and $v_3^{SN2}$ divide the region in Figure 3(a) into five parts, in which phase portraits of $x$ and $y$ are illustrated in Figure 3(b)–(f). Only a stable equilibrium appears for $v_3 = 0.45 < v_3^{SN1}$ in Figure 3(b) and $v_3 = 0.61 > v_3^{SN2}$ in Figure 3(f). Two equilibria coexist at $v_3 = v_3^{SN1}$ in Figure 3(c) and $v_3 = v_3^{SN2}$ in Figure 3(e). There are three equilibria for $v_3$ that vary between $v_3^{SN1}$ and $v_3^{SN2}$ in Figure 3(d)).

Furthermore, the stability and property of these equilibria $E(x_{*}, y_{*})$ in Figure 3(b)–(f) are listed in Table 3 to verify three cases in Table 2. The properties of unstable non-hyperbolic equilibria and unstable saddle are consistent with the conditions of cases 4 and 7 in Table 2, respectively. Other equilibria are asymptotically stable nodes, which correspond to case 1 in Table 2.

4.2.   Hopf bifurcation

In this section, we try to explore the conditions under which a positive equilibrium $E(x_{*}, y_{*})$ loses the stability through Hopf bifurcation under some parametric restriction. Here, considering $v_3$ as the bifurcation parameter, we shall establish the conditions under which the system $(2.1)$ experiences Hopf bifurcation at the positive equilibrium $E(x_{*}, y_{*})$ when $v_3$ crosses the critical value $v_3^{HB}$.

The first condition is that the Jacobian matrix $J(E, v_3^{HB})$ has a pair of purely imaginary eigenvalues, that is $\mathrm{tr}(J(E, v_3^{HB})) = 0$ and $\mathrm{det}(J(E, v_3^{HB})) > 0$, which correspond to

Besides, the transversality condition that ensures the changes of stability of the positive equilibrium through non–degenerate Hopf bifurcation is ${\frac{d Re(\lambda_{i})}{d v_3}\mid_{v_3 = v_3^{HB}} \neq 0}$, i.e.,

Lastly, the first Lyapunov number $\Gamma$ at the equilibrium $E(x_{*}, y_{*})$ is computed to analyze the stability of the limit cycle. By the transformation $X = x - x_{*}$, $Y = y - y_{*}$, the system $(2.1)$ becomes

where

According to the formula in ^[21], the first Lyapunov number $\Gamma$ is given as follows

where $\Phi = {a}_{10}{b}_{01}-{a}_{01}{b}_{10}$.

Therefore, we get the following theorem.

Theorem 4.2. The system $(2.1)$ experiences Hopf bifurcation at the positive equilibrium $E(x_{*}, y_{*})$ when $v_3$ crosses the critical value $v_3^{HB}$ with the conditions (HB.1) and (HB.2). A supercritical (subcritical) Hopf bifurcation occurs for the first Lyapunov number $\Gamma < 0 (>0)$ in (HB.3). Here $v_3 ^{HB} = -\frac{d_2(d_1+d_2)(k_2 + x_{*})^2(k_3^2+x_{*}^2)}{k_2 v_2 x_{*}^2}+\frac{2d_2 k_1^2 v_1 x_{*}(k_2 + x_{*})^2(k_3^2+x_{*}^2)}{ k_2 v_2 x_{*}^2(k_1^2 + x_{*}^2)^2} -\frac{r_2 (k_3^2+x_{*}^2)}{ x_{*}^2}.$

Note that the sign of the first Lyapunov number $\Gamma$ cannot be determined directly due to its complex expression; we illustrate Hopf bifurcation of the system $(2.1)$ through the following numerical examples to verify the correctness of Theorem 4.2.

Two supercritical Hopf bifurcation points $HB_{sup1}$ and $HB_{sup2}$, and a subcritical Hopf bifurcation point $HB_{sub}$ are shown in bifurcation diagrams of $x$ with respect to $v_3$ in Figure 4 (a)–(b), where black solid and dashed lines respectively represent stable and unstable equilibria while green and purple lines denote stable and unstable limit cycles. For $v_1 = 0.8153$, $v_2 = 0.2$, $d_1 = 0.015$, $d_2 = 0.04$, $r_1 = 0.03$, $r_2 = 0.002$, $k_1 = 7$, $k_2 = 3$ and $k_3 = 7$, the first critical value $v_3^{HB} = 0.2096$ labeled by $HB_{sup1}$, where the conditions, $E(x_{*}, y_{*}) = (5.1628, 1.8957)$, $F^{\prime}(x_{*}) = -34.7231 < 0$, ${\frac{\text{d}tr(J(E))}{\text{d}v_3} = 0.1612 > 0}$ and $\Gamma = -0.1455 < 0$, imply the appearance of a stable limit cycle. Then the stable limit cycle disappears at another critical value $v_3^{HB} = 0.2589$ labeled by $HB_{sup2}$ with $E(x_{*}, y_{*}) = (3.3111, 1.2332)$, $F^{\prime}(x_{*}) = -19.3163 < 0$, ${\frac{\text{d}tr(J(E))}{\text{d}v_3} = -0.2941 < 0}$ and $\Gamma = -0.9058 < 0$. For $v_3 = 0.25$ between $HB_{sup1}$ and $HB_{sup2}$, a stable limit cycle surrounding an unstable equilibrium is illustrated in the phase portrait of $x$ and $y$ in Figure 4 (c). Beside, Figure 4 (b) shows a subcritical Hopf bifurcation point $HB_{sub}$ at $v_3^{HB} = 0.1852$ with the conditions $E(x_{*}, y_{*}) = (4.4394, 4.1289)$, $F^{\prime}(x_{*}) = -2.1014 < 0$, ${\frac{\text{d}tr(J(E))}{\text{d}v_3} = 2.7028 > 0}$ and $\Gamma = 0.4856 > 0$. An unstable limit cycle surrounding a stable focus coexists with a saddle and a node for $v_3 = 0.1838$ in the phase portrait in Figure 4 (d).

4.3.   Bogdanov – Takens bifurcation

Apart from codimension–1 saddle–node and Hopf bifurcation, the system $(2.1)$ may undergo codimension–2 Bogdanov–Takens bifurcation at the positive equilibria $E(x_*, y_*)$, which will be analyzed by considering $v_3$ and $d_2$ as bifurcation parameters in the next section.

Let $v_3^{BT}$ and $d_2^{BT}$ be two critical values of $v_3$ and $d_2$, at which $\mathrm{det}(J(E))\left|_{(v_3, d_2) = (v_3^{BT}, d_2^{BT})} = 0 \right.$ and $\mathrm{tr}(J(E))\left|_{(v_3, d_2) = (v_3^{BT}, d_2^{BT})} = 0.\right.$ Perturbing $v_3$ and $d_2$ by $v_3 = v_3^{BT} + \mu_1$ and $d_2 = d_2^{BT} + \mu_2$ with $\mu_1$ and $\mu_2$ in a small neighborhood of $(0, 0)$, the system $(2.1)$ becomes

Transforming the equilibrium $E(x_{*}, y_{*})$ to the origin $(0, 0)$ by $z_1 = x - x_{*}$ and $z_2 = y - y_{*}$, the system $(4.1)$ becomes

where

and $z = (z_1, z_2)^T$, $\mu = (\mu_1, \mu_2)^T$. $R_{i}(z, \mu) = O({\|z\|}^3)$ $(i = 1, 2)$ denotes the power series as to $z_1$, $z_2$ with the order 3 and more. The coefficients of $R_{i}(z, \mu)(i = 1, 2)$ and $a_{ij}$, $b_{ij}$ smoothly depend on $\mu_1$ and $\mu_2$. According to $b_{00}(0) = 0$, we rewrite the system $(4.2)$ at $\mu_1 = 0$ and $\mu_2 = 0$ in the following form

where

Let $a_{10}(0) = a_{10}$, $a_{01}(0) = a_{01}$, $b_{10}(0) = b_{10}$ and $b_{01}(0) = b_{01}$. Since $J_0$ has two zero eigenvalues, we get $a_{10} + b_{01} = 0$ and $a_{10}b_{01} = a_{01}b_{10}$. Then, we choose $\alpha_0 = \left(1, -\frac{a_{10}}{a_{01}}\right)^{T}, \; \alpha_1 = \left(0, \frac{1}{a_{01}}\right)^{T}$ and $\beta_0 = \left(1, 0\right)^{T}, \; \beta_1 = \left(a_{10}, a_{01}\right)^{T}$ as the eigenvector and generalized eigenvector with zero eigenvalues for $J_0$ and $J_0^T$, respectively, which satisfy $\langle \alpha_0, \beta_0\rangle =$ $\langle \alpha_1, \beta_1\rangle = 1$ and $\langle \alpha_1, \beta_0\rangle = \langle \alpha_0, \beta_1\rangle = 0$. The linearly independent vectors $\alpha_0$ and $\alpha_1$ form a basis of $\mathbb{R}^2$. Thus, we make the following transformation

i.e., $\gamma_1 = z_1$, $\gamma_2 = a_{10}z_1 + a_{01}z_2$. Then the system $(4.2)$ becomes

where

and $R_{3, 4}(\gamma, \mu) = O({\|\gamma\|}^3)$ ($\gamma = (\gamma_1, \gamma_2)^T$). Due to $a_{10} + b_{01} = 0$ and $a_{10}b_{01} = a_{01}b_{10}$, we obtain $d_{00}(0) = d_{10}(0) = d_{01}(0) = 0$.

Next, making the transformations $\eta_1 = \gamma_1$ and $\eta_2 = \gamma_2 + c_{20}(\mu){\gamma_1}^2 + c_{11}(\mu)\gamma_1\gamma_2 + R_3(\gamma, \mu)$, system $(4.3)$ becomes

where

and $R_5(\eta, \mu) = O({\| \eta\|}^3)$, $\eta = (\eta_1, \eta_2)^T$. Moreover, we have

And we assume that

Then making a coordinate shift $\eta_1 = \omega_1 + \xi(\mu), \; \eta_2 = \omega_2,$ where $\xi(\mu) \approx-\frac{e_{01}(\mu)}{e_{11}(0)},$ the system $(4.4)$ reduces to

where

and $R_{6}(\omega, \mu) = O(\|\omega\|^3)$ ($\omega = (\omega_1, \omega_2)^T$).

Next, introducing a new time variable $\tau_{1}$ by $\mathrm{d} t = \left(1 + \theta(\mu) \omega_1\right) \mathrm{d} \tau_1$, $\theta(\mu) = -f_{02}(\mu)$, system $(4.5)$ becomes

in which

and $R_{7}(\omega, \mu) = O(\|\omega\|^3)$.

If the following condition holds,

then we introduce the new variables ${\tau_2 = \left|\frac{h_{20}(\mu)}{h_{11}(\mu)}\right| \tau_1}$, ${y_1 = \frac{h_{11}^2(\mu)}{h_{20}(\mu)} \omega_1}$, ${y_2 = \operatorname{sign}\left(\frac{h_{20}(\mu)}{h_{11}(\mu)}\right) \frac{h_{11}^3(\mu)}{h_{20}^2(\mu)} \omega_2}$, the system $(4.7)$ is changed to

where

and $R_{8}(y, \mu) = O(\|y\|^3)$, $y = (y_1, y_2)^T$.

If the following transversality condition holds,

the system $(4.1)$ experiences Bogdanov–Takens bifurcation when $\mu = \left(\mu_1, \mu_2\right)$ is in a small neighborhood of $(0, 0)$ based on the results in ^[22]. According to the above discussion, the following theorem is obtained.

Theorem 4.3. The system $(2.1)$ experiences codimension–2 Bogdanov–Takens bifurcation at the positive equilibrium $E(x_{*}, y_{*})$ as $\left(v_3, d_2\right)$ varies near $\left(v_3^{BT}, d_2^{BT}\right)$ and the conditions (BT.1)–(BT.4) are satisfied. The local representations of the bifurcation curves are given as follows:

(i) the saddle–node bifurcation curve $S N = \left\{\left(\mu_1, \mu_2\right) \mid 4 m_{00} - m_{10}^2 = 0\right\}$;

(ii) the Hopf bifurcation curve $H = \left\{\left(\mu_1, \mu_2\right) \mid m_{00} = 0, \; m_{10} < 0\right\}$;

(iii) the homoclinic bifurcation curve $H L = \left\{\left(\mu_1, \mu_2\right) \mid m_{00} = {-\frac{6}{25}} m_{10}^2 + O\left(m_{10}^2\right), \; m_{10} < 0\right\}$.

Theorem 4.3 is verified by codimension–2 bifurcation diagram of $v_3$ and $d_2$ in Figure 5 and phase portraits of $x$ and $y$ in Figure 6 with $r_1 = 0.011$, $r_2 = 0.001$, $v_1 = 0.63$, $v_2 = 0.06$, $k_1 = 4.5$, $k_2 = 2$, $k_3 = 4$ and $d_1 = 0.034$. Through numerical simulation, we find that system $(2.1)$ undergoes Bogdanov–Takens bifurcation at $(v_3, d_2) = (0.3419136, 0.04331208)$, where the conditions (BT.1)–(BT.4) in Theorem 4.3 are as follows,

$e_{11}(0) = -0.002882 \neq 0$, $h_{20}(0) = d_{20}(0) = -0.000125 \neq 0$, and $s = sign(0.042991) = +1$. Moreover, the local representations of the bifurcation curves are given as follows,

(i) the saddle–node bifurcation curve $SN$

(ii) the Hopf bifurcation curve $H$

(iii) the homoclinic bifurcation curve $HL$

Figure 5 illustrates the curves $SN$, $H$ and $HL$ in the ($\mu_1$, $\mu_2$) parameter plane. These curves divide the small neighborhood of the origin (0, 0) in the $(\mu_1, \mu_2)$ plane into four parts, in which the phase portraits of $x$ and $y$ are given in Figure 6 with the trajectory in blue solid lines and the stable and unstable equilibria in red solid and hollow dots, respectively. In order to see the variation around Bogdanov–Takens bifurcation point, phase portraits of $x$ and $y$ are given in insets in Figure 6(a)–(f).

(a) The system ($2.1$) has a cusp of codimension–2 Bogdanov–Takens bifurcation point $E_2^{BT}$ with another stable node for $\left(\mu_1, \mu_2\right) = (0, 0)$ (see Figure 6$(a)$).

(b) The system ($2.1$) has a unique stable node when $(\mu_1, \mu_2)$ lies in region I (see Figure 6$(b)$).

(c) A saddle and an unstable focus coexist with another stable node when $(\mu_1, \mu_2)$ enters into region II from region I through the branch $SN^{-}$ of the curve $SN$ (see Figure 6$(c)$).

(d) The unstable focus becomes stable and is surrounded by an unstable limit cycle when $(\mu_1, \mu_2)$ crosses the subcritical Hopf bifurcation curve $H$ into region III (see Figure 6$(d)$).

(e) A homoclinic loop occurs for $(\mu_1, \mu_2)$ on the curve $HL$ (see Figure 6$(e)$).

(f) The unstable limit cycle disappears and three stable equilibria are left when $(\mu_1, \mu_2)$ crosses the $HL$ curve into region IV (see Figure 6$(f)$).

5.   Hopf bifurcation with time delay

In this section, we will investigate the effect of time delay on the stability of the positive equilibrium $E(x_{*}, y_{*})$ in system ($2.1$).

The system ($2.1$) is linearized at $E(x_{*}, y_{*})$ as follows,

where

The characteristic equation of the linearized system $(5.1)$ is

Assuming that $\lambda = i \omega(\omega > 0)$ is the root of Eq $(5.3)$ and $\tau = \tau_1+\tau_2$, we get the following equation

Separating the real and imaginary parts of Eq $(5.4)$ results in

Thus, $\cos(\omega\tau)$ and $\sin(\omega\tau)$ are given by

which implies that

The discriminant of Eq $(5.7)$ is

Therefore, Eq $(5.7)$ has two different roots $\omega_1^{2}$ and $\omega_2^{2}$, and $\omega_1^{2} + \omega_2^{2} = - (p_{11}^2+p_{22}^2) > 0$, $\omega_1^{2}\omega_2^{2} = (p_{22} p_{11})^2 - {(p_{12}p_{21})}^2$. Therefore, if $(p_{22} p_{11})^2 - {(p_{12}p_{21})}^2 > 0$, Eq $(5.7)$ has a purely imaginary root $i\omega_0$ and

Then, according to Eq $(5.6)$, the critical value of $\tau$ is given as follows,

Let

Next, we will verify the transversality condition ${sign\left\{\left.\left[\frac{d \text{Re}(\lambda(\tau))}{d \tau}\right]\right|_{ \tau = \tau_0}\right\}\neq 0.}$ By differentiating both sides of Eq $(5.3)$ with respect to $\tau$, we get

At $\tau = \tau_0$, we have

Hence,

Finally, we have the following theorem based on the Hopf bifurcation theorem ^[37].

Theorem 5.1. Let $\tau_0$ and $p_{11}, p_{12}, p_{21}, p_{22}$ be defined by Eqs (5.11) and (5.2), If $(p_{22} p_{11})^2 - {(p_{12}p_{21})}^2 < 0$, the positive equilibrium $E(x_{*}, y_{*})$ of system $(2.1)$ is asymptotically stable for $\tau \in(0, \tau_0]$ and the system $(2.1)$ undergoes Hopf bifurcations at $\tau = \tau_0$.

Theorem $5.1$ is verified by the bifurcation diagram of $x$ with respect to $\tau_2$ in Figure 7 (a) and phase portraits of $x$ and $y$ in Figure 7 (b)–(d) with $v_1 = 0.86, v_2 = 0.6, v_3 = 0.89, d_1 = 0.017, d_2 = 0.49, r_1 = 0.17, r_2 = 0.004, k_1 = 1.07, k_2 = 0.03$ and $k_3 = 0.42$. For these parameters, the system ($2.1$) has a unique equilibrium $E(x_{*}, y_{*}) = (0.2126, 0.3785)$ and $(p_{22} p_{11})^2 - {(p_{12}p_{21})}^2 = -0.4940 < 0$, then $\omega_0 = 0.7677$, $\tau_0 = 0.4677$ in Theorem $5.1$. Based on Theorem $5.1$, the positive equilibrium $E(x_{*}, y_{*}) = (0.2126, 0.3785)$ is locally asymptotically stable for $\tau \in(0, \tau_0)$ and the system $(2.1)$ undergoes supcritical Hopf bifurcation at $\tau = \tau_0$, which is accord with bifurcation diagram in Figure 7 (a). In Figure 7 (a), black solid and dashed lines respectively denote stable and unstable equilibria while green dots represent the maxima and minima of stable limit cycle, supcritical Hopf bifurcation $HB_{sup}$ occurs at $\tau_2 = 0.2677 = \tau_0 - \tau_1$ with $\tau_1 = 0.2$. Besides, in the phase diagrams of $x$ and $y$ in Figure 7 (b)–(d), red solid and hollow dots denote stable and unstable equilibria, respectively, and blue and green lines respectively denote the trajectory and a stable limit cycle, the positive equilibrium $E(x_{*}, y_{*})$ is stable for $\tau_2 = 0.1 < \tau_0 - \tau_1$ in Figure 7 (b); and it lose stability at $\tau_2 = 0.2677 = \tau_0-\tau_1$ in Figure 7 (c); then it become unstable one with the appearance of a stable limit cycle at $\tau_2 = 0.35 > \tau_0-\tau_1$ in Figure 7 (d).

6.   Conclusions

In this paper, we are mainly concerned with the stability and bifurcation of system $(2.1)$ without and with time delay $\tau$. For $\tau = 0$, the existence of positive equilibria of system $(2.1)$ are analyzed through the Descartes' rule of signs to obtain the conditions under which the system $(2.1)$ has a unique positive equilibrium in Theorem 3.1, which are verified by the nucllines of $x$ and $y$ in system $(2.1)$ in Figure 2. The stability of positive equilibria of system $(2.1)$ without time delay is presented in Table 2 due to the complex expression of the determinant and trace of the Jacobian matrix. For positive equilibria, selecting $v_3$ as a bifurcation parameter, the conditions of saddle-node and Hopf bifurcation in system $(2.1)$ are given in Theorems 4.1 and 4.2, which gives the first Lyapunov number that determines the stability of limit cycle. These two theorems are verified by codimension-1 bifurcation diagram of $x$ with respect to $v_3$ in Figures 3 and 4, which include saddle-node bifurcation and supercritical and subcritical Hopf bifurcation, respectively. Furthermore, by choosing two parameters $v_3$ and $d_2$ in system $(2.1)$ as bifurcation parameters, we prove that the system exhibits codimension-2 Bogdanov-Takens bifurcation under the conditions in Theorem 4.3, which is obtained by calculating universal unfolding near the cusp. Theorem 4.3 is verified by codimension- 2 bifurcation diagram in Figure 5, which includes Bogdanov-Takens bifurcation point generating the curves of saddle-node, Hopf and homoclinic bifurcation. For $\tau \neq 0$, we find time delay induce the superscribe Hopf bifurcation under the conditions in Theorem 5.1, which is verified by the bifurcation diagram of $x$ with respect to $\tau_2$ and phase portraits of $x$ and $y$.

Our results give rigorous mathematical proofs of the stability and bifurcation for a two-dimensional p53 GRN to expand the understanding of p53 GRN. Besides, the theoretical analyses of GRNs with noise and space will be further to explored^[38,39,40]. In recent years, fractional-order differential equations (FODEs) are used to describe GRNs because they possess memory, after-affects and hereditary properties, which are more compatible with reality than the integer-order differential equations^{[41,42,43,44]}. Therefore, it is worth to explore the stability and bifurcation of GRNs described by three or four dimensional FODEs.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence tools in the creation of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Nos. 12062017 and 12262025), Natural Science Foundation of Inner Mongolia Autonomous Region of China (Grants 2021ZD01), and Program for Innovative Research Team in Universities of Inner Mongolia Autonomous Region of China (No. NMGIRT2008). The authors acknowledge the reviewers for their valuable reviews and kind suggestions.

Conflict of interest

The authors declare that there is no conflict of interest.