Mathematical modeling of tumor growth: the MCF-7 breast cancer cell line

1.   Introduction

The most important tumor suppressor p53 has a key role in regulating a number of cellular responses upon various stresses including DNA damage and oncogene activation ^[1,2]. Depending on the type and severity of the damage, p53-mediated responses include cell cycle arrest, cellular senescence and apoptosis, which can be triggered through regulating the expression of different downstream genes ^[3,4]. These responses are closely related to p53 dynamics, where p53 keeps low level in unstressed cells, whereas a high level of p53 induces apoptosis to destroy irreparably damaged cells ^[5,6]. More interestingly, the oscillation dynamics of p53 may initiate cell cycle arrest or apoptosis depending on the number of p53 pulses ^[7,8]. Therefore, it is meaningful to explore p53 dynamics for understanding cellular response in response to different conditions. Furthermore, p53 dynamics can be investigated through simulating the mathematical models of p53 gene regulatory network with p53, it upstream regulators, downstream target genes and feedback loops ^[9,10].

A number of mathematical models depicting p53 gene regulatory network have been established to explore p53 dynamics in response to different conditions from theoretical analyses and numerical simulations ^[7,8,11]. Typically, an integrated model of p53 gene regulatory network with four modules was proposed by Zhang et al., who demonstrated that the number of p53 pulses determined cell fates upon DNA damage ^[7,8]. Based on the integrated model and experiment results, a mathematical model of p53 gene regulatory network regulated by Programmed Cell Death 5(PDCD5) has been developed so as to reveal that PDCD5 modulates cell fate decisions as a co-activator of p53 ^[12,13]. Furthermore, the effects of rate constants in the p53 gene regulatory network on p53 dynamics are explored through bifurcation analyses and potential landscapes ^[14,15]. However, the transcriptional and translational processes in the model are seen as instant. Actually, slow and complex transcriptional and translational processes as the basic steps of gene expression in cells are usually related to long time delay ^[16]. Time delays in various networks, such as neural networks and disc dynamos with viscous friction, have important impact on the dynamics in the network, which are explored through stability and bifurcation analyses ^{[17,18,19,20]}. Also, more researches explore the effect of time delays in some p53 gene regulatory networks on p53 oscillation ^[21]. However, the effect of time delay on p53 oscillation in p53 gene regulatory network regulated by PDCD5 has not been considered.

Inspired by these considerations, we introduce a time delay of the transcriptional and translational processes of Mdm2 activated by p53 into p53 gene regulatory network regulated by PDCD5 in ref. ^[22]. The local stability of positive equilibrium of the model are demonstrated through analyzing the characteristic equations of the corresponding linearized system for two cases without and with time delay. Also, these stability are verified through time series, the distribution of characteristic roots, one-parameter bifurcation diagram and potential landscapes. Then, two-parameter bifurcation diagrams of time delay with the concentration of PDCD5 and four typical rate parameters, $v_{p53}$, $d_{p53}$, $v_{Mdm2}$ and $d_{Mdm2}$ are carried out to explore their effect on the oscillation dynamic of p53.

This paper is organized as follows. The p53 gene regulatory network regulated by PDCD5 is given in Section 2. In Section 3, we analyze the local stability of positive equilibrium of the model by means of stability theory for two cases without and with time delay. Then numerical simulations and bifurcation analyses are given in Section 4. Finally, we conclude the paper in Section 5.

2.   A model of p53 gene regulatory network with time delay

The model of p53 gene regulatory network regulated by PDCD5 is illustrated in Figure 1, which was established in Ref. ^[22]. Double strand breaks (DSBs) induced by various stresses activate ATM monomers, which further phosphorylate p53 at Ser-$15$ in the nucleus(p53) and Mdm2 at Ser-395 in the cytoplasm ($\mathrm{Mdm2_{cyt}^{395P}}$). Phosphorylated p53 in the nucleus inhibits ATM activity through activating it inhibitor Wip1 ^[23,24,25]. Additionally, Mdm2 protein in the cytoplasm($\mathrm{Mdm2_{cyt}}$) transforms with $\mathrm{Mdm2_{cyt}^{395P}}$ and Mdm2 protein in the nucleus ($\mathrm{Mdm2_{nuc}}$). Furthermore, $\mathrm{Mdm2_{cyt}^{395P}}$ stimulates p53 translation through interacting with p53 mRNA while $\mathrm{Mdm2_{nuc}}$ promotes p53 degradation by binding p53 and increasing it ubiquitination ^[26]. Besides, PDCD5 functions as a positive regulator of p53 through accelerating the degradation of $\mathrm{Mdm2_{nuc}}$ and inhibiting the degradation of p53 by $\mathrm{Mdm2_{nuc}}$ ^[22]. Here, phosphorylated p53 accelerates the production of $\mathrm{Mdm2_{cyt}}$ through promoting the expression of Mdm2 gene. As we know, it takes some time from gene transcription to protein production, so we introduce a time delay into the interaction of p53 with $\mathrm{Mdm2_{cyt}}$. The model equations of the p53 gene regulatory network with time delay are given as follows.

where

Here ${x(t)}$, ${y(t)}$, ${z(t)}$, ${w(t)}$ and ${m(t)}$ represent the concentration at time $t$ of active p53, $\mathrm{Mdm2_{cyt}}$, $\mathrm{Mdm2_{cyt}^{395P}}$, $\mathrm{Mdm2_{nuc}}$ and ATM proteins. The biological significance and values of all parameters in the model are given in Table 1.

3.   Local stability analysis

In this section, the local stability of positive equilibrium $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ are analyzed through stability theory for without and with time delay. Let $\overline{x}(t) = x(t)-X^{*}$, $\overline{y}(t) = y(t)-Y^{*}$, $\overline{z}(t) = z(t)-Z^{*}$, $\overline{w}(t) = w(t)-W^{*}$, $\overline{m}(t) = m(t)-M^{*}$, the linearization of system (2.1) at $E^{*}$ is given as follows,

where

The system (3.1) can be written as the form of matrix:

where

Using Laplace transform,

Eq (3.2) can be given as follows:

Without loss of generality, let

Then, the characteristic equation of the linear system (3.1) is give as follows:

That is:

where

3.1.   The case $\tau=0$

For $\tau = 0$, Eq (3.5) becomes

The conditions for local stability of $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ are give in Theorem 3.1 according to the Routh-Hurwitz criterion ^[27].

Theorem 3.1. $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ is locally asymptotically stable when all roots of Eq $(3.6)$ have negative real parts, the following conditions (H) should be satisfied.

(H):

3.2.   The case $\tau\neq0$

For $\tau\neq0$, let $i\omega(\omega > 0)$ be a root of Eq (3.5), so $\omega$ satisfies the following equation:

Separating the real and imaginary parts, we have

Squaring both sides of each equation and then adding them up leads to

Let $z = w^{2}$, $p = a^{2}_{1}-2a_{2}$, $q = a^{2}_{2}+2a_{4}-2a_{1}a_{3}$, $r = a^{2}_{3}-2a_{2}a_{4}+2a_{1}a_{5}-A^{2}_{2}$, $u = a^{2}_{4}-2a_{3}a_{5}-A^{2}_{1}+2A_{0}A_{2}$ and $v = a^{2}_{5}-A^{2}_{0}$. Then Eq (3.8) can be written as

Some lemmas for establishing the distribution of positive real roots of Eq (3.9) are given as follows ^[28].

Lemma 3.2. If $v < 0$, then Eq $(3.9)$ has at least one positive root.

Proof. Denote $h(z) = z^{5}+pz^{4}+qz^{3}+rz^{2}+uz+v$. Since $h(0) = v < 0$ and $\mathop{lim}\limits_{z\rightarrow +\infty}h(z) = +\infty$. Hence, there exists a $z_{0} > 0$ such that $h(z_{0}) = 0$. The proof is complete.

For $v\geq0$, considering the following equation:

Let $z = y-\frac{p}{5}$, then Eq (3.10) is changed into:

where

If $q_{1} = 0$, then the four roots of Eq (3.11) are as follows:

where $\Delta_{0} = p^{2}_{1}-4r_{1}$. Thus $z_{i} = y_{i}-\frac{p}{5}$, $i = 1, 2, 3, 4$ are the roots of Eq (3.10). Then we have the following lemma.

Lemma 3.3. Suppose that $v\geq0$ and $q_{1} = 0$.

(i) If $\Delta_{0} < 0$, then Eq $(3.9)$ has no positive real roots.

(ii) If $\Delta_{0}\geq0$, $p_{1}\geq0$ and $r_{1} > 0$, then Eq $(3.9)$ has no positive real roots.

(iii) If (i) and (ii) are not satisfied, then Eq $(3.9)$ has positive real roots if and only if there exists at least one $z^{*}\in\{z_{1}, z_{2}, z_{3}, z_{4}\}$ such that $z^{*} > 0$ and $h(z^{*})\leq0$.

Proof. (ⅰ) If $\Delta_{0} < 0$, then Eq (3.10) has no real roots. Since $\mathop{lim}\limits_{z \rightarrow +\infty}h^{'}(z) = +\infty$, so $h^{'}(z) > 0$ for $z\in\mathbb{R}$. Hence, $h(0) = v\geq0$ leads to $h(z)\neq0$ when $z\in(0, +\infty)$.

(ⅱ) If $\Delta_{0}\geq0$, $p_{1}\geq0$ and $r_{1} > 0$, then $h^{'}(z)$ has no zero in $(-\infty, +\infty)$. It is similar to (ⅰ) that $h(z)\neq0$ for any $z\in(0, +\infty)$.

(ⅲ) The sufficiency holds obviously. We will prove the necessity. If $\Delta_{0}\geq0$, then Eq (3.11) has only four roots $y_{1}, y_{2}, y_{3}, y_{4}$, so Eq (3.10) has only four roots $z_{1}, z_{2}, z_{3}, z_{4}$ at least $z_{1}$ is a real root. We assume that $z_{1}, z_{2}, z_{3}, z_{4}$ are all real. This means that $h(z)$ has at most four stationary points at $z_{1}, z_{2}, z_{3}, z_{4}$. If it is not true, we will have that either $z_{1}\leq0$ or $z_{1} > 0$ and $min\{h(z_{i}):z_{i}\geq0, i = 1, 2, 3, 4\} > 0$. If $z_{1}\leq0$, then $h^{'}(z)\neq0$ in $(0, +\infty)$. Since $h(0) = v\geq0$ is the strict minimum of $h(z)$ for $z\geq0$ which means $h(z) > 0$ when $z\in(0, +\infty)$. If $z_{1} > 0$ and $min\{h(z_{i}):z_{i}\geq0, i = 1, 2, 3, 4\} > 0$, since $h(z)$ is the derivable function and $\mathop{lim}\limits_{z\rightarrow +\infty}h(z) = +\infty$, then we have $\mathop{min}\limits_{z > 0}h(z) = min\{h(z_{i}):z_{i}\geq0, i = 1, 2, 3, 4\} > 0$. The proof of the necessity is done. We finish the proof.

If $q_{1}\neq0$. Thinking about the solution of Eq (3.11):

that is

Let

So, Eq (3.12) has the following three roots

Let $s^{*} = s_{1}\neq p_{1}$ (since $p_{1}$ is not a root of Eq (3.12)). Equation (3.11) is equivalent to

For Eq (3.13), Eq (3.12) indicates that the formula in square brackets is a perfect square.

If $s^{*} > p_{1}$, then Eq (3.13) is

After factorization, we get

and

Let

Then we obtain the four roots of Eq (3.11)

So $z_{i} = y_{i}-\frac{p}{5}$, $i = 1, 2, 3, 4$ are the roots of Eq (3.10). Thus we can get the following lemma.

Lemma 3.4. Suppose that $v\geq0$, $q_{1}\neq0$ and $s^{*} > p_{1}$.

(i) If $\Delta_{2} < 0$ and $\Delta_{3} < 0$, then Eq $(3.9)$ has no positive real roots.

(ii) If (i) is not satisfied, then Eq $(3.9)$ has positive real roots if and only if there exists at least one $z^{*}\in\{z_{1}, z_{2}, z_{3}, z_{4}\}$ such that $z^{*} > 0$ and $h(z^{*})\leq0$.

Proof. We omit the proof since it is similar to Lemma 3.3.

Finally, if $s^{*} < p_{1}$, then Eq (3.13) is

Let $\overline{z} = \frac{q_{1}}{2(p_{1}-s^{*})}-\frac{p}{5}$, we have the following lemma.

Lemma 3.5. Suppose that $v\geq0$, $q_{1}\neq0$ and $s^{*} < p_{1}$, then Eq $(3.9)$ has positive real roots if and only if $\frac{q^{2}_{1}}{4(p_{1}-s^{*})^{2}}+\frac{s^{*}}{2} = 0$, $\overline{z} > 0$ and $h(\overline{z})\leq0$.

Proof. If Eq (3.13) has a real root $y_{0}$, which satisfies

so

Therefore, Eq (3.13) has a real root $y_{0}$ if and only if $\frac{q^{2}_{1}}{4(p_{1}-s^{*})^{2}}+\frac{s^{*}}{2} = 0$. The rest of the proof is similar to Lemma 3.3, we omit it.

Lemmas 3.2, 3.3, 3.4 and 3.5 give the sufficient and necessary conditions of existence of positive roots of Eq (3.9).

Suppose that Eq (3.9) has five positive roots $z_{k}, i = 1, 2, 3, 4, 5$. Then Eq (3.8) has five positive roots $\omega_{k} = \sqrt{z_{k}}$, $i = 1, 2, 3, 4, 5$. From Eq (3.8), we obtain the corresponding $\tau_{j}^{(k)} > 0$ and the characteristic equation (3.5) has purely imaginary roots.

It is obvious that $\mathop{lim}\limits_{j\rightarrow \infty}\tau_{j}^{(k)} = \infty$ for $k = 1, 2, 3, 4, 5$. Next we can define

On the basis of the above analysis, we obtain the conditions that all roots of Eq (3.5) have negative real parts in the following lemma.

Lemma 3.6. (i) If one of the followings holds:(a) $v < 0$; (b) $v\geq0$, $q_{1} = 0$, $\Delta_{0}\geq0$ and $p_{1} < 0$ or $r_{1}\leq0$ and there exist $z^{*}\in\{z_{1}, z_{2}, z_{3}, z_{4}\}$ such that $z^{*} > 0$ and $h(z^{*})\leq0$; (c) $v\geq0$, $q_{1}\neq0$, $s^{*} > p_{1}$, $\Delta_{2}\geq0$ or $\Delta_{3}\geq0$ and there exist $z^{*}\in\{z_{1}, z_{2}, z_{3}, z_{4}\}$ such that $z^{*} > 0$ and $h(z^{*})\leq0$; (d) $v\geq0$, $q_{1}\neq0$, $s^{*} < p_{1}$, $\frac{q_{1}^{2}}{4(p_{1}-s^{*})^{2}}+\frac{s^{*}}{2} = 0$, $\overline{z} > 0$ and $h(\overline{z})\leq0$, then all roots of the characteristic equation $(3.5)$ have negative real parts when $\tau\in(0, \tau_{0})$.

(ii) If the conditions (a)-(d) of (i) are all not satisfied, then all roots of Eq $(3.5)$ have negative real parts for all $\tau > 0$.

Proof. Lemmas 3.2, 3.3, 3.4 and 3.5 shows that Eq (3.5) has no roots with zero real part for all $\tau > 0$ if the conditions (a)–(d) of (ⅰ) are not satisfied. If one of (a)–(d) holds, when $\tau\neq\tau^{(k)}_{j}$, Eq (3.5) has no roots with zero real part and $\tau_{0}$ is the minimum value of $\tau$, which make equation (3.5) has purely imaginary roots. We get the conclusion of the lemma.

Next, we want to prove the system (2.1) can undergo Hopf bifurcation at $\tau = \tau_{0}$.

Let

be the root of Eq (3.5) satisfying $\alpha(\tau_{0}) = 0$, $\omega(\tau_{0}) = \omega_{0}$. We have the next lemma.

Lemma 3.7. Suppose $h^{'}(z_{0})\neq0$. If $\tau = \tau_{0}$, then $\pm i\omega_{0}$ is a pair of simple purely imaginary roots of Eq (3.5).

Proof. If $i\omega_{0}$ is not simple, then $\omega_{0}$ must satisfy

This implies that

Meanwhile, from Eq (3.7) we konw that $\omega_{0}$ satisfies

Thus, we have

According to Eq (3.9) and $\omega^{2}_{0} = z_{0}$, Eq (3.19) becomes

Notice that $h(z) = z^{5}+pz^{4}+qz^{3}+rz^{2}+uz+v$ and $h^{'}(z) = 5z^{4}+4pz^{3}+3qz^{2}+2rz+u$, Eq (3.20) becomes $((A_{0}-A_{2}z_{0})^{2}+A_{1}^{2}z_{0})h^{'}(z_{0}) = 0$. we obtain a contradiction to the condition $h^{'}(z_{0})\neq0$. This proves the conclusion.

Next, we will get the conditions for $\frac{d Re(\lambda(\tau_{0}))}{d\tau} > 0$ in Lemma 3.8.

Lemma 3.8. If the conditions of Lemma 3.6 (i) are satisfied, then $\frac{d Re(\lambda(\tau_{0}))}{d\tau} > 0$.

Proof. Differentiating both sides of Eq (3.5) with respect to $\tau$ gives

Therefore

If $\frac{dRe\lambda(\tau_{0})}{d\tau} < 0$, then equation (3.5) has a root with positive real part for $\tau < \tau_{0}$ and close to $\tau_{0}$, which contradicts Lemma 3.6 (ⅰ). This completes the proof.

To sum up, we get the following theorems.

Theorem 3.9. Let $\omega_{0}$, $z_{0}$, $\tau_{0}$, and $\lambda(\tau)$ be defined by Eq $(3.15)$ and $(3.16)$, respectively.

(i) If the conditions (a)–(d) of Lemma 3.6 are not satisfied, then all the roots of Eq $(3.5)$ have negative real parts for all $\tau > 0$.

(ii) If one of the conditions (a)–(d) of Lemma 3.6 is satisfied, then all the roots of Eq $(3.5)$ have negative real parts when $\tau\in(0, \tau_{0})$; when $\tau = \tau_{0}$ and $h^{'}(z_{0})\neq0$, then $\pm i\omega_{0}$ is a pair of simple purely imaginary roots of Eq $(3.5)$ and all other roots have negative real parts. In addition, $\frac{dRe\lambda(\tau_{0})}{d\tau} > 0$ and Eq $(3.5)$ has at least one root with positive real part when $\tau\in(\tau_{0}, \tau_{1})$, when $\tau_{1}$ is the first value of $\tau > \tau_{0}$ such that Eq $(3.5)$ has purely imaginary root.

Theorem 3.10. Let $\omega_{0}$, $z_{0}$, $\tau_{0}$, and $\lambda(\tau)$ be defined by Eq $(3.15)$ and $(3.16)$, respectively.

(i) If the conditions (a)–(d) of Lemma 3.6 are not satisfied, then the positive equilibrium point $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ of system (2.1) is asymptotically stable for all $\tau > 0$;

(ii) If one of the conditions (a)–(d) of Lemma 3.6 is satisfied, then the positive equilibrium point $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ of system (2.1) is asymptotically stable when $\tau\in(0, \tau_{0})$;

(iii) If conditions (a)–(d) of Lemma 3.6 is satisfied, and $h^{'}(z_{0})\neq0$, system (2.1) undergoes a Hopf bifurcation at $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ when $\tau = \tau_{0}$.

In summary, Theorems 3.1 and 3.10 give the local stability of positive equilibrium $E^{*} = (X^{*}, Y^{*}, Z^{*}, W^{*}, M^{*})$ of system (2.1) and the conditions undergoing Hopf bifurcation without and with time delay, respectively. In the following, some numerical simulations are given to verify the correctness of these theorems.

4.   Numerical simulation and bifurcation analysis

4.1.   Numerical simulation

In this section, Theorems 3.1 and 3.10 are verified through time series of the concentrations of active p53($x$), total Mdm2($y+z+w$), ATM($m$) and the eigenvalues of Eq (3.5) for $\tau = 0$ and $\tau\neq0$. Besides, one-parameter bifurcation diagram of the concentrations of active p53($x$) with respect to $\tau$ and potential landscapes of active p53($x$) and total Mdm2($y+z+w$) are given to verify the correctness of Theorem 3.10. These figures are realized by MATLAB software.

In case Ⅰ of $\tau = 0$, a positive equilibrium $E^{*} = (3.0492, 0.0770, 0.1184, 0.0858, 2.2877)$ of system (2.1) is obtained for parameters in Table 1. Substituting these parameters into $\Delta_{i}(i = 1, 2, 3, 4, 5)$, we have $\Delta_{1} = 1.9910 > 0$, $\Delta_{2} = 2.2370 > 0$, $\Delta_{3} = 0.6136 > 0$, $\Delta_{4} = 0.0197 > 0$ and $\Delta_{5} = 0.0134 > 0$. The condition (H) in Theorem 3.1 is satisfied, so the positive equilibrium $E^{*}$ is asymptotically stable, which is verified through numerical simulations in Figure 2. Figure 2(a) shows that the concentrations of $x$, $y+z+w$ and $m$ firstly increase and then converge to the positive equilibrium $E^{*}$. Also, all roots of the characteristic equation (3.6) at the positive equilibrium point $E^{*}$ have negative real parts, as shown in Figure 2(b).

For case Ⅱ $\tau\neq0$, substituting the same positive equilibrium $E^{*}$ with $\tau = 0$ and the parameters in Table 1 into Eq (3.8), we get

Equation (4.1) has a positive root $\omega_{0} = 0.2342$, then $z_{0} = \omega_{0}^{2} = 0.0548$ and $\tau_{0} = 2.1217$ from Eq (3.8). Note that $v = -1.2282\times10^{-4} < 0$ and $h^{'}(z_{0}) = 0.0065\neq0$, the conditions (ⅱ) and (ⅲ) in Theorem 3.10 are satisfied, so the positive equilibrium $E^{*}$ is asymptotically stable when $\tau\in(0, \tau_{0})$ and system (2.1) undergoes Hopf bifurcation at $\tau_{0} = 2.1217$, where $E^{*}$ loses stability and a stable limit cycle appears. The correctness of Theorem 3.10 is verified by the following numerical simulations in Figures 3–6. Figure 3(a) describes one-parameter bifurcation diagram of $x$ with respect to $\tau$, where black solid and dashed lines are stable and unstable equilibria, and green solid dots are the maxima and minima of the stable limit cycle. As is shown in Figure 3(a), system (2.1) converges to the positive equilibrium $E^{*}$ for $\tau < \tau_{0}$ and undergoes Hopf bifurcation at HB point with $\tau = \tau_{0}$, then the positive equilibrium $E^{*}$ loses stability and a stable limit cycle appears for $\tau > \tau_{0}$. Also, the stability of the positive equilibrium $E^{*}$ is further verified through the maximum real part of the eigenvalues of the characteristic equation (3.5) in Figure 3(b), where maximum real part is negative for $\tau < \tau_{0}$ and zero at $\tau = \tau_{0}$ then positive for $\tau > \tau_{0}$.

Next, Figures 4 and 5 show time series of the concentrations of active p53($x$), total Mdm2($y+z+w$), ATM($m$) and the eigenvalues of Eq (3.5) for three typical time delays $\tau = 1 < \tau_{0}$, $\tau = 2.1217 = \tau_{0}$ and $\tau = 2.5 > \tau_{0}$ to further verify the stability of system (2.1). As shown in Figure 4, system (2.1) converges to the stable equilibrium $E^{*}$ at $\tau = 1$ (Figure 4(a)) and shows damped oscillation at $\tau = \tau_{0}$ (Figure 4(b)), then is in a continuous oscillating state at $\tau = 2.5$ (Figure 4(c)). Accordingly, all roots of the characteristic equation (3.5) have negative real part at $\tau = 1$ (Figure 5(a)), and a pair of roots with zero real part appear at $\tau = 2.1217$ (Figure 5(b)), and some roots have positive real part at $\tau = 2.5$(Figure 5(c)).

Furthermore, the global stability of p53 dynamics are explored through the potential landscapes of active p53($x$) and total Mdm2($y+z+w$) in Figure 6. The potential landscape $U = -ln(P_{ss})$ proposed by Wang et al. is related with the the steady state probability distribution $P_{ss}$ of the system (2.1) ^[29]. $P_{ss}$ can be obtained through simulating the corresponding stochastic differential equations of system (2.1) for a long time over many independent runs. The stochastic differential equations are obtained through adding Gaussian and white noise $\zeta(t)$ into right hand of each equation of system (2.1), where $< \zeta(t)\zeta(t') > = 2D\delta(t-t')$, and $< \zeta(t) > = 0$, $D$ is noise strength.As shown in the Figure 6(a), the potential landscape has a global minimum that corresponds to the global stable steady state. However, the potential landscapes in Figure 6(b) and 6(c) exhibit closed ring valleys which correspond to the stable limit cycle in system (2.1). The potential landscape at $\tau = 2.1217$ has wider attractive regions than the one at $\tau = 2.5$ since system (2.1) shows damped oscillation at $\tau = 2.1217$ and continuous oscillation at $\tau = 2.5 > \tau_{0}$.

Numerical simulations in this section are in agreement with the theoretical results in Section 3. In the following, we will further explore the effect of time delay $\tau$ and typical parameters on the oscillation dynamics of p53 from bifurcation perspective.

4.2.   Bifurcation analyses

To investigate how time delay $\tau$ affect p53 oscillation, we perform two-parameter bifurcation analyses with respect to time delay $\tau$ and several typical parameters, including the concentration of PDCD5($p$), maximum production and degradation rates of p53 and Mdm2: $v_{p53}$, $d_{p53}$, $v_{Mdm2}$, $d_{Mdm2}$. Figures 7–11(a) show two-parameter bifurcation diagrams, where black solid lines $hb$ represent Hopf bifurcation curves. For each two-parameter bifurcation diagram, several typical one-parameter bifurcation diagrams are given in Figures 7–11(b)–(e), where black solid and dashed lines are stable and unstable steady states, respectively, and green solid dots are maxima and minima of stable limit cycles.

Figure 7(a) shows two-parameter bifurcation diagram with time delay $\tau$ and the concentration of PDCD5($p$). Hopf bifurcation curves $hb_{1}$ and $hb_{2}$ on the figure divide the parameter region into three parts. p53 reaches a lower stable steady state for parameters on the left of the curve $hb1$, where $p$ is less than 0.710 regardless of the value of $\tau$. Then p53 oscillation appears for parameters between these two curves $hb_{1}$ and $hb_{2}$, where p53 always oscillates for any $\tau$ when $0.710 < p < 0.974$ while p53 can oscillate when $\tau$ increases to certain value for $0.974 < p < 11.9$. Also the value of $\tau$ on Hopf bifurcation curve is more larger for larger $p$. However, p53 reaches a higher stable steady state for any $\tau$ when $p > 11.9$. For intuitiveness, one-parameter bifurcation diagrams of p53 concentration with respect to $\tau$ for four typical $p$ are given in Figure 7(b)–(e).

Figures 8–10(a) show similar two-parameter bifurcation diagrams with respect to $\tau$ and $v_{p53}$, $d_{p53}$, $v_{Mdm2}$, respectively, where the Hopf bifurcation curves $hb$ divides the parameter plane into two parts. p53 exhibits oscillation inside the curve $hb$ and reaches a stable steady state outside the curve $hb$. As shown in Figures 8–10(a), time delay can induce p53 oscillation for moderate value of the rate constants while it has no effect on p53 dynamics for smaller and larger value of these parameters. Also, the value of time delay on Hopf bifurcation curve increases with the increase of $v_{p53}$ and $v_{Mdm2}$ while the value of time delay decreases with the increase of $d_{p53}$. Furthermore, from one-parameter bifurcation diagrams in Figures 8–10(b)–(d), p53 reaches a lower stable steady state for smaller $v_{p53}$ and $v_{Mdm2}$ and a higher stable steady state for larger $v_{p53}$ and $v_{Mdm2}$, while the case is opposite for $d_{p53}$.

Besides, two-parameter bifurcation diagram with respect to $\tau$ and $d_{Mdm2}$ is shown in Figure 11(a), where p53 exhibits oscillation for any $\tau$ when $d_{Mdm2}$ is less than 0.0146 and reaches a stable steady state for any $\tau$ when $d_{Mdm2}$ is larger than 0.0355 while p53 changes from a stable steady state to oscillation with the increase of $\tau$ when $d_{Mdm2}$ varies between 0.0146 and 0.0355. These dynamics are shown on one-parameter bifurcation diagrams in Figure 11(b)–(d).

Therefore, bifurcation analyses further explore the effect of time delay $\tau$ on p53 oscillation.

5.   Conclusions

The oscillation dynamics of p53 plays a key role in deciding cell fate after DNA damage. The time delay in transcriptional and translational processes is an important factor to induce oscillation. In this study, we explore the effect of time delay of Mdm2 gene express in p53 gene regulatory network regulated by PDCD5 on p53 oscillation through stability theory and bifurcation analyses. Firstly, we give theorems that provide the conditions on the appearance of p53 oscillation for without and with time delays through stability analysis theory. Then, these theorems are verified through numerical simulations, including time series of p53 concentrations, one-parameter bifurcation diagram of p53 concentration versus time delay, the roots of the characteristic equations, and the potential landscapes of p53 and total Mdm2 concentrations. Furthermore, the effect of time delay and several typical parameters $p$, $v_{p53}$, $d_{p53}$, $v_{Mdm2}$ and $d_{Mdm2}$ in system (2.1) on p53 oscillation are explored through two-parameter bifurcation diagrams, which give boundary curves of p53 oscillation. Our results show that time delay can induce p53 oscillation for moderate value of most parameters while it has no effect on p53 dynamics for smaller and larger value of these parameters. Taken together, analyzing the effect of time delay on p53 oscillation will help us understand the mechanism of p53 oscillation.

In this study, we explore the effect of time delay of Mdm2 gene express in p53 gene regulatory network regulated by PDCD5 on p53 oscillation. Our results have demonstrated that suitable time delay can induce p53 oscillation for moderate value of most parameters. However, the effect of multiple time delays in p53 gene regulatory network on p53 oscillation should be further explored. Besides, noise is common in gene regulatory work, it is worthwhile to explore the combined effect of noise and time delay on p53 dynamics. Also, it is meaningful to investigate the dynamics of the other gene regulatory networks or neural networks with time delay and noise for understanding their biological mechanism.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Nos. 12062017 and 11702149), Natural Science Foundation of Inner Mongolia Autonomous Region of China (Grant 2021ZD01) and Collaborative Innovation Center for Grassland Ecological Security (Jointly Supported by the Ministry of Education of China and Inner Mongolia Autonomous Region). The authors acknowledge the reviewers for their valuable reviews and kind suggestions.

Conflict of interest

The authors declare there is no conflicts of interest.