Vector control for the Chikungunya disease

1.   Introduction

Numerous studies have shown that P53 plays a key role in determining cell survival and death. More than 50 percent of cancer patients are detected with mutations in the P53 gene ^[1]. High levels of the mutant p53 protein were found in mouse tumors as early as 1983 ^[2]. In subsequent studies, abnormal expression of P53 protein was also found in a variety of human tumors ^[3,4,5]. In normal human cells, P53 is generally kept at a low level due to the downregulation of Mdm2 (murine double minute2) ^[6], which also helps cells avoid premature aging and apoptosis ^[7]. In general, with cells stimulated by hypoxia, DNA damage and impair of telomere function, P53 is rapidly produced and activated ^[8]. As a transcription factor, P53 prevents (or at least alleviates) damage caused by mutations, thereby regulating the expression of genes involved in various cellular functions, including cell cycle arrest, DNA damage repair, and apoptosis ^[9]. Previous investigation has proved that the specific function of P53 closely related to its level. When a weak stimulus is felt, P53 rapidly rises to a moderate level, thereby making cell cycle arrest and avoiding the inheritance of abnormal DNA replication to the next generation. Until the damage was eliminated, the expression level of P53 returned to the normal level ^[10,11,12]. Moreover, a large damage can cause the level of P53 rises to a higher state in order to induce apoptosis ^[13,14]. As we all know, P53 promotes transcription of its target gene Mdm2 (a P53-specific E3 ubiquitin ligase), which contributes to the degradation of P53 in return ^[15,16]. In fact, many cancer treatment programs are related to the network that is based on the P53-Mdm2 negative feedback loop ^[17,18]. Clinically, when the P53 protein is mutated, for example, mutations in the P53 protein make cells resistant to chemotherapy drugs ^[19,20]. However, when the P53 mutant protein is knocked out or interfered with, chemotherapy-induced apoptosis can be restored ^[21,22].

In recent years, many scholars dedicated to investigate the P53-Mdm2 feedback network. In 2000, the researchers theoretically modeled the P53 damped oscillations and also found that the negative feedback of P53-Mdm2 with time delay can led to the occurrence of P53 oscillation ^[23]. In 2005, Ma and his collaborators explored the mechanism of P53 pulse generation in theory at the IBM Watson Research Center in the United States ^[24]. In 2007, the Tysons team constructed a theoretical model of the P53 network response to DNA damage and explored the dynamics of P53 pulse-determining cell fate whose view is that the number of P53 pulses determining cell fate ^[25]. The more perfect result is that for mild lesions, a small amount of P53 pulses cause cell cycle arrest, and cells can return to normal after repair is completed. Whereas for severe damage, persistent P53 pulses kill the irreparable cell by activating downstream apoptosis program ^[25]. In 2009, Zhang et al. constructed a more comprehensive P53 network model and concluded that the fate of cells is closely related to the P53 pulse ^[26]. In 2011, the study of Zhang et al. shows that a sequential predominance of distinct feedback loops may elicit multiple-phase dynamical behaviors ^[27]. Only the research of Tingzhe Sun et al. shows that the kinetics of P53 under non-stressed conditions is triggered by an excitatory mechanism, and cells can only respond fully in the face of severe damage ^[28]. Although a large number of studies have been done, most of them were based on stress (DNAdamage, radiation, hypoxia, etc.). However, researches on non-stress situations are rarely involved. Moreover, as for these systems, the transcription, translation and integration into the polymer are usually seen as instant. In fact, the impact of the time delay of transcription and translation on the system in real biological processes is not negligible ^[29], especially multiple simultaneous time delays occur at the same time. Admittedly, there are several mathematical models of p53 that include time delays ^[7,28,30,31], which just illustrated that the objective existence of time delay is recognized. These studies have perfectly studied the dynamic mechanism of the p53 pulse from different aspects. But without exception, they did not specially focus on the role of time delay on the system. For example, Purvis et al. proposed a delayed computational model and identified a sequence of precisely timed drug additions that alter p53 pulses to instead produce a sustained p53 response that suggested that the drug is efficient and can directly influencing cellular fate decisions ^[7]. Even though there are three delays including Nutlin-3 activity delay, p53-induced expression of Mdm2 and Wip1 (Wild type p53-induced phosphatase 1) delays, the effect of time delay on the dynamic process is not received attention. Zhang et al. established a complex model without time delay to study the signal threshold to elicit cell death and the time delay between signal and response ^[30]. It should be noticed that the time delay there refers to the time required from the onset of external signal to the apoptosis triggered by the system. A mathematical model with two time delays that respectively represent the time consumed on the transcription of Mdm2 and wip1 activated by p53 is proposed by Sun et al. The model described the negative feedback loop Wip1 –$\shortmid$ ATM (mutated in the disease ataxia telangiectasia) $\rightarrow$ P53 $\rightarrow$ Wip1 except to the core negative feedback loop formed by P53 and Mdm2 to study the excitability of p53 pulses and the effect of Wip on the p53 pulses ^[28]. Chong et al. developed a delayed mathematical model that incorporated the molecular interactions in the core regulation of p53 and the apoptosis initiation module involving Puma, Bcl2 and Bax and demonstrated how molecular interactions and stress signaling of the p53 network dictate cell fate decisions ^[31]. Although there are 10 time delays, their dynamic effects on the p53 pulse are also lacked. Here, we consider the following facts. P53 is regulated by ARF (ADP-ribosylation factor), which exerts an anti-cancer role by combining with Mdm2 (or $Mdm2_ {p}$), and is considered to be the main regulator of P53, so it is well known to promote the activation and stability of P53 ^[32]. Previous studies have shown that Growth factors in serum not only induce the transcriptional activation of Mdm2 ^[33], but the activation of Akt (protein kinase B) which phosphorylates Mdm2 promoting $Mdm2_{p}$ nuclear accumulation ^[34]. P53 trans-activates PTEN (gene of phosphate and tension homology deleted on chromsome ten), a tumor suppressor that dampens Akt activation by dephosphorylating PIP3 ^[34,35]. Especially, this paper places emphasis on the effect of positive feedback loop P53 $\rightarrow$ PTEN –$\shortmid$ Akt $\rightarrow$ Mdm2 –$\shortmid$ P53 and two time delays that respectively represent the p53-induced expression of Mdm2 and PTEN delays on the dynamic of the system. On the above interaction facts among biological molecules, we focus on the dynamics of P53 under non-stressed conditions based on the model which was given by Xinyu Tian, Bo Huang, Xiao-Peng Zhang et al. in 2017 ^[36].

The novelty of this paper is that most studies have considered the dynamic behavior of P53 under stress, while our study considered the non-stress state and studied the delay in detail. First of all, the effects of time delay on the P53 oscillation is analyzed in two different cases including there are one time delay and two time delays co-exist by using Hopf bifurcation theory. Previous studies have well demonstrated that negative feedback loop with time delay can produce oscillations and the network which includes coupled positive and negative feedback loops can also generate oscillations when negative feedback plays the main role ^[37,38]. It is quite important that we find the relative length of the double time delays is likely to control the dynamic properties of the system in the coupled positive and negative feedback loops. In addition, we consider the influences of several important chemical reaction rates that relating to P53 directly or indirectly.

This paper is organized as follows. In Section 2, we introduced the materials and method. In Section 3, we give numerical simulation results of the influences of time delays and key parameters on the P53 oscillation behaviors. Discussion is made in Section 4. Conclusion is drawn in Section 5. The theoretical derivations about the properties of Hopf bifurcations are presented in Supplementary.

2.   Materials and method

As can be seen from Figure 1. Growth factors promote Mdm2 accumulation and Akt activation, while active Akt can promote the transfer of Mdm2 into the nucleus and phosphorylate Mdm2. The phosphorylated Mdm2 can be denoted by $Mdm2_{p}$, which has a more pronounced ability to degrade P53 ^[34]. $Mdm2_{p}$ can be further phosphorylate to $Mdm2_{pp}$. The latter ability of $Mdm2_{pp}$ to degrade P53 is very weak ^[39]. The increase of the total amount of Mdm2 in the nucleus indirectly leads to a decrease in the level of P53. In verse, P53 promotes the transcription and translation of Mdm2. Thus, a negative feedback loop of P53-Mdm2 is formed. On the other hand, PTEN can inhibit the activity of Akt and convert it to inactive $Akt_0$ to reduce the conversion of Mdm2 to $Mdm2_p$ and thus indirectly promote the accumulation of P53. P53 can also promotes PTEN transcription. So the interations construct a positive feedback loop P53-PTEN-ATK-Mdm2. ARF binds to two forms of Mdm2 and inhibits the degradation of P53. Among these regulation process, it should be noticed that from DNA to protein, all transcription, transport and translation take time ^[40,41,42], which is ignored by the previous theoretical study ^[36]. We considered that both Mdm2 and PTEN could be upregulated by a means of a transcription activation from P53 ^[15,43,44]. Based on the related literatures mentioned above, we developed a delayed ODEs system to model the P53-PTEN-Akt-Mdm2 network, which is shown as (2.1).

Next, we mainly imply the stability and bifurcation theory to analyze the dynamic behavior of system (2.1). First, by selecting time delays as the bifurcation parameters and analyzing the distribution of the roots of the corresponding characteristic equation in detail, we mainly investigate the stability of positive equilibrium and the existence of Hopf bifurcation. Further, their direction and properties of the oscillatory dynamic are also researched by using normative theory and central manifold method. Moreover, the dynamic influences of four important parameters, namely S (concentration of serum), $k_{3}$ (Rate constant for $Mdm2_{p}$ dephosphorylation), $k_{10}$ (Basal expression rate of PTEN) and $k_{14}$ (Rate constant for PTEN-induced Akt dephosphorylation) on the P53-PTEN-Akt-Mdm2 network are discussed. In addition, some numerical simulations are executed by using the software Mathematica 11 to verify the theoretical predictions. Finally, the discussion and conclusion are drawn.

where $\tau_{1}$ and $\tau_{2}$ respectively represent transcriptional delays of Mdm2 and PTEN triggered by P53, and the meaning of all the variables are given in Table 1 and related parameters are given in Table 2.

3.   Results

In this section, we will focus on the effects of two time delays including $\tau_1$ and $\tau_2$ and four important parameters related to P53 including S, $k_3$, $k_{10}$, $k_{14}$ on the system oscillations. All the numerical simulations are performed via Mathematical 11.

3.1.   Effects of time delay $\tau_1$ on the system

For the legibility of the text, we only give short explanation for the theoretical results and the expatiatory and tedious theoretical derivations and results are put in the Supplementary Section. Then numerical simulations of time delay $\tau_1$ on the system (2.1) are given when $\tau_{2} = 0$.

3.1.1.   Theoretical analysis about Hopf bifurcation and its properties

There is a positive equilibrium solution $E^{*} = (x^{*}, y^{*}, z^{*}, u^{*}, v^{*}, p^{*}, q^{*})$ by setting the right-hand side of the equation equal to zero in (2.1). Then we can get the linearized equation near the equilibrium as follows

where $c_{i} (i = 1, ..., 20)$ are decided by the parameters in system (2.1) and the positive equilibrium solution $E^{*} = (x^{*}, y^{*}, z^{*}, u^{*}, v^{*}, p^{*}, q^{*})$ and the actual formulas are given in the Supplementary Section.

Furthermore, the characteristic equation of the system (3.1) is given as

where $m_{i} (i = 1, ..., 14)$ are depended on the parameters $c_{i} (i = 1, ..., 20)$ and the actual formulas can be found in the Supplementary Section.

Suppose $\pm i \omega (\omega > 0)$ is a pair of pure virtual roots of the characteristic equation (3.2), then $\omega$ satisfies the following condition

Through a large amount of calculation which is displayed in the Supplementary Section, it is found that there exists unique positive root $w_0$ of equation (3.3). Moreover,

Define $\tau _1^0 = \min \left\{ {\tau _1^{(j)} > 0} \right\}_{j = 0}^{ + \infty }$. Thus, when $\tau = \tau_1^{0}$ equation (3.2) has a pair of pure imaginary roots $\pm i{w_0}$.

Furthermore, it is verified that the condition

is holds.

It indicates that the root of characteristic equation(3.2) changes though negative to positive that corresponds to the positive equilibrium of system (2.1) is changed from stability to instability as the time delay ${\tau _1}$ increases and crosses the threshold $\tau _1^0$. What is more, Hopf bifurcation of the positive equilibrium solution occurs at the threshold ${\tau _1} = \tau _1^0$. Concretely, the system is asymptotically stable when ${\tau _1} \in [0, \tau _1^0)$, while it is unstable when ${\tau _1} > \tau _1^0$. Thus, $\tau _1^0$ is a vital critical value and represents the length of time delay need for the system arising oscillation behavior.

In addition, according to the normal form of functional differential equation, the following four index to depict the properties of the above Hopf bifurcation are obtained.

Among them, ${\mu _2}$ determines the direction of Hopf bifurcation. If ${\mu _2} > 0$, then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exist for $\tau > \tau _1^0$, and while it is subcritical when ${\mu _2} < 0$ and the bifurcating periodic solutions exist for $\tau < \tau _1^0$. The stability of the bifurcating periodic solutions is depended on ${\beta _2}$. Especially, the bifurcating periodic solutions in the center manifold are stable when ${\beta _2} < 0$, while they are unstable if ${\beta _2} > 0$. The period of the bifurcating periodic solutions is determined by ${T_2}$. If ${T_2} > 0$, the period increases. When ${T_2} < 0$, the period decreases.

3.1.2.   Numerical simulation on the oscillation derived from Hopf bifurcation

In this section, we will perform numerical simulation to study the affect of time delay $\tau_{1}$ on the system (2.1). On the one hand, it will further verify the correctness of the above theoretical derivation. On the other hand, it will give a more intuitively understand on the importance of time delay $\tau_{1}$. Except for the special parameter changes, all others are taken from Table 2. It is easy to calculated that there is a approximate positive equilibrium $E^*(0.56548, 0.37000, 0.11052, 0.00451, 0.39694, 0.18240, 0.33258)$. The critical value of $\tau _1$ is calculated as $\tau_1^0 = 0.70454767059$ by the formula (3.4). The four important indexes are also calculated as follows according to the formulas (3.6).

In order to vividly show the affect of $\tau_1$, we used different values to obtain a series of images of Mdm2+$Mdm2_p$, P53, PTEN and Akt concentrations over time as shown in Figure 2. We can see that the system is stable when $\tau_1 = 0$. But as its increases, the system begins to change from stable to progressively stable. When $\tau_1$ reaches the value $\tau_1^0 \approx {\rm{0}}{\rm{.70455}}$, the system experiences Hopf bifurcation and changes from asymptotic stability to sustained oscillation, which is consistence with the above theoretical conclusion. As $\tau_1$ continues to increase, the system still keeps oscillating continuously. At the same time, the magnitude and period of the oscillation of the system become larger. As identified by the above calculated values $\mu _2 > 0$, $T _2 > 0$ and $\beta _2 > 0$, which suggested that the Hopf bifurcation is subcritical and the period and amplitude of the oscillations increases as $\tau_1$ increases. Therefore, the behavior of the system can be indirectly controlled by regulating the time delay $\tau_1$. Moreover, we plot their images by setting discrete time delays $\tau_1 = 0.75$, $\tau_1 = 0.8$, $\tau_1 = 0.85$, $\tau_1 = 0.95$, $\tau_1 = 1.05$ separately as shown in Figure 3. From Figure 3(a) and (b) we can see that both of the magnitude and period of P53 and $Mdm2_{tot}$ increases as time delay ${\tau_1}$ gradually varies from small to big. Figure 3(c) and (d) respectively shows that magnitude and period of PTEN and Akt are also added as time delay ${\tau_1}$ increases. Interestingly, the oscillation valleys of PTEN nearly remain unchanged. But for the the Akt, its oscillation peaks almost stayed the same. For more visual observation, the change in the period and amplitude are drawn in Figure 4, which illustrated that both of the period and amplitude of four molecule including P53, $Mdm2_{tot}$, PTEN and Akt are increased as the time delay $\tau_1$ increases. The increase degrees of the period for the four molecules are almost the same as shown in Figure 4(a), which is coordinate to the experiments. However, the increase degrees of the amplitude for the four molecules are different as shown in Figure 4(b), which is $Mdm2_{tot}$, P53, PTEN and Akt in the order from small to large and the biological meanings are mysterious.

3.2.   Effects of $\tau_1$ and $\tau_2$ on the system

In the Supplementary, we have only proved the existence of Hopf bifurcation when ${\tau _1}$ varies. According to the above subsection, we identified that system (2.1) undergoes a supercritical Hopf bifurcation at the positive equilibrium $E^*$ when ${\tau_1} = \tau_1^0$. However, we do not consider the time delay ${\tau_2}$ of P53-promoting transcription of PTEN for the convenience of the study. In fact, the time delay does exist. Therefore, it is necessary to study the co-regulation of ${\tau_1}$ and ${\tau_2}$. In this subsection, we take $\tau_1$ to satisfy the oscillation condition and then make the $\tau_2$ varies. By comparing Figure 5 (b) and (c), (a) and (e), (d) and (f), we get that the system oscillation is enhanced by $\tau_1$, which is agree with the above analysis. Intestinally, if we comparing Figure 5 (a) and (b), (c) and (d), (e) and (f), the inverse phenomenon occurs, which is he system oscillation is suppressed by $\tau_2$. Therefore, the relative sizes between $\tau_1$ and $\tau_2$ determines the degree of oscillation. To our knowledge, transcription and splicing of intron sequences increase the time required for the extension and splicing of genes to mature RNA ^[40,48,49]. Thus, an effective way to vary the time delay is changing the number of introns within the gene. Takashima et al. previously showed that deletion of all three introns within the Hes7 gene reduces the time delay by 19 min and completely abolishes oscillatory expression ^[50]. Based on conclusion that each consecutive intron splicing in mammalian cells required about 0.4 to 7.5 minutes ^[51,52], the number of introns within the genes of Mdm2 and PTEN can be calculated in order to achieve the desired length of time delay $\tau_1$ and $\tau_2$, which requires specific experimental verification. Time delays $\tau_1$ and $\tau_2$ can flexibly regulate the oscillation dynamics of P53. Experiments have shown that P53 oscillation dynamics is closely related to many diseases including cancer ^[9,25,26,27]. In particular, P53 oscillation can promote the cells to repair damage. And if the damage is too severe to repair, it can also promote the cell trigger the programmed cell death. In this way, the cells can avoid the wrong genetic information from being inherited to the next generation of cells. Therefore, the above research results provide a method for regulating P53 oscillation dynamics through time delays, which may provide a new insight for the treatment of cancer and other diseases.

3.3.   Dependence of dynamics on the important chemical reaction rates related to P53

In this subsection, in order to more comprehensive understand the P53 oscillation dynamic, a Sobol sensitivity analysis is performed to examine the influence of each parameter $k_{i} (i = 1, 2, ..., 14)$ on the P53 gene regulatory network using the ODE sensitivity package in R software. The evolution processes of the parameter sensitivity indices including the first order and total Sobol' indices for the P53 gene regulatory system (2.1) are shown in Figure 6. The first order Sobol' indices measures the pure interaction influence of the variables, while the total sensitivity index measures the influence of the variables including all interactions of any order and it reflects the individual influence of all factors plus every interaction effect ^[53]. Using the total Sobol' indices for the concentration of MDM2 denoted by $x$, a clear order of series parameters have been established: $k_3$ has the largest total effect, followed by S, $k_7$, $k_{11}$, $k_5$, $k_8$, $k_2$, $k_{12}$, $k_{14}$, $k_{10}$ and so on. In addition, considering the influence of some parameters that have important biological significance and can cause the essential changes in the system, this subsection will focus on the four parameters S, $k_3$, $k_{10}$, and $k_{14}$ to elaborate the cooperative effects among them and two time delays on the system dynamics. In concretely, growth factor S in serum not only induce the transcriptional activation of Mdm2, but the activation of Akt which, phosphorylates Mdm2, promoting its nuclear accumulation and further enhance inhibition to P53 ^[36]. Both of Mdm2 and $Mdm2_p$ have the ability to induce the P53 degradation. $k_3$ represents the rate constant for $Mdm2_p$ dephosphorylation. $Mdm2_p$ has the stronger degradation ability on P53 than Mdm2 ^[24,54]. In other words, the greater $k_3$, the greater dephosphorylation rate of $Mdm2_p$, this will convert $Mdm2_p$ into Mdm2 and increase P53 indirectly. Thus, S and $k_3$ can indirectly affect the kinetics of P53. $k_{10}$ is the basal expression rate of PTEN and $k_{14}$ refers to the rate constant for PTEN induced Akt dephosphorylation. Both of them reflect the strength of the positive feedback loop P53 $\rightarrow$ PTEN –$\shortmid$ Akt $\rightarrow$ Mdm2 –$\shortmid$ P53 including P53 and PTEN. The larger the $k_{10}$, the more PTEN is produced, which in turn promotes the positive feedback loop of P53 leading to an increase in P53. The larger the $k_{14}$, the less Akt is phosphorylated, the more P53 produced due to the weak repression of P53 from Mdm2. Obviously, both of the values $k_{10}$ and $k_{14}$ are important to the system dynamic behavior.

Now, we will discuss the the cooperative effects of S, $k_3$, $k_{10}$ and $k_{14}$ with the time delays respectively, and the simulation results are obtained in Figure 7, 8, 9, and 10. Here, we set the time delay ${\tau _1} = 0.8$ and only with ${\tau _2}$ varies. Especially, we consider three values ${\tau _2}$, which is 0, 0.1, and 0.15 in turn. On the one hand, if we horizontally observe the Figure 7, it is obviously that the parameter S can induce and block the oscillation of the four molecules in all three sets of the time delays. The oscillation range occurs at a neighborhood of $S\approx0.035$. On the other hand, if we vertically observe the Figure 7, it can be find that the time delay $\tau_2$ has the inhibitory role on the system oscillation in all three sets of the time delays, which is in accordance with the above results. S represents the Serum concentration that promotes the concentration of Mdm2 so that S indirectly against the accumulation of P53 concentration. Too much or too little serum concentration will not cause the produce of system oscillation, only moderate serum concentration will induce it. For the parameter $k_3$, $k_{10}$ and $k_{14}$, we can get the same result with S. First, observing the Figures 8, 9, and 10 horizontally, we will find that regardless which value set of the time delays, the three parameters $k_3$, $k_{10}$ and $k_{14}$ can induce the oscillation of the four molecules or make it disappear. Then, looking at the Figures 8, 9, and 10 vertically, it can also be found that no matter which value of the parameters $k_3$, $k_{10}$ and $k_{14}$ takes, the time delay $\tau_2$ has a suppressive effect on system oscillation. In particular, the parameter range of the system oscillation is $k_3\approx11.2$, $k_{10}\approx0.04$ and $k_{14}\approx24$ in sequence when other parameters keep fixed. In summary, under the influence of time delays, the system has strong sensitivity to four parameters and is also strictly controlled by the them. The $k_3$, $k_{10}$ and $k_{14}$ are all indirectly enhance the P53 concentration. Similarly with S, neither too large nor too small parameters can cause system oscillations, and only appropriately large parameters can do it.

In order to obtain a more visual perspective on the cooperative effect of two time delays and four important parameters, bifurcation diagrams of P53 with respect to four important parameters under different cases of time delay are given in Figure 11. The existence of Hopf bifurcation indicates that the system oscillation occurs. It can be seen that there are three common characters among the Figures 11(a), (b), (c) and (d). One of them is that the Hopf bifurcation occurs only when $\tau_{1} = 0.8, \tau_{2} = 0$, while under other two cases including $\tau_{1} = 0, \tau_{2} = 0$ and $\tau_{1} = 0.8, \tau_{2} = 0.1$ the bifurcation phenomenon can not be observed. This led us to conclude that time delay $\tau_{1}$ has the opposite effect on the P53 oscillation production than $\tau_{2}$. Concretely, $\tau_{1}$ can prompt the P53 oscillation occur, whereas $\tau_{1}$ tends to hold back the oscillation. From a biological point of view, P53 oscillation means that the cell choices to repair the damage in response to stresses. The second common character is that bistability occurs only when $\tau_{1} = 0.8, \tau_{2} = 0$, while under other two cases including $\tau_{1} = 0, \tau_{2} = 0$ and $\tau_{1} = 0.8, \tau_{2} = 0.1$ the bistability phenomenon can not be found. The bistability of P53 is an important feature, which corresponding to the switch of cell fate decision. More specifically, the low steady-state means that the cell stays on a normal cell cycle state, while the high steady-state means that the cell triggers appoptosis when upon extra stimulations. Both of the DNA damage repair and appoptosis are the preventive measures of cells to avoid disease or cancers ^{[8,9,10,11,12,13,14]}. Therefore, $\tau_{1}$ provides two ways to prevent the occurrence of diseases. The third common character is that time delays $\tau_{1}$ and $\tau_{2}$ is helpful to accumulate or silent the activity of P53 as can be seen by comparing the pink and blue lines, which suggested that the time delays $\tau_{1}$ and $\tau_{2}$ have the ability to flexibly control the state of P53. On the other hand, there are also two main differences among the Figures 11(a), (b), (c) and (d). The first one is the parameter $s$ has the opposite role than other three parameters. As $s$ increases, the P53 concentration is declined. But as $k_{3}, k_{10}$ and $k_{14}$ increases, the P53 concentration is raised. The second one is that the switch in the Figure 11(c) is irreversible, while others are reversible because the parameter $k_{10}$ can not be negative. This indicates that the apoptosis triggered by the cooperation of $k_{10}$ and $\tau_{1}$ is irredeemable death. The above analysis again shows that p53 oscillation is regulated by two time delays and four important parameters. The specific conclusion is a prediction, which needs further verification by subsequent biological experiments.

4.   Discussion

In this paper, we studied a mathematical model of the core regulatory network of cell fate decision that responds to growth factor instead of DNA damage. It contains seven major components and two key time delays. In our model, the dynamics of the P53-PTEN-Akt-Mdm2 network is studied using Hopf bifurcation theory and numerical simulation. The results help us understand how normal cells respond to growth factor stimuli, thereby preventing the occurrence and development of disease. Here, we just studied the effect of time delay on the P53-PTEN-Akt-Mdm2 network, whereas the downstream networks including survival and apoptosis module have not been considered. Beyond that, It has been reported that aberrant expression of miRNAs and a global decrease of their levels are often observed in multiple types of human cancer cells. In addition, the P53 pathway are often crosstalk with E2F and NF-kB signalling pathway. Therefore, more extensively network should be further investigated for comprehensive understand the potential mechanism of P53 related network in the future.

5.   Conclusion

In this paper, we proposed a delayed mathematical model to describe the P53-PTEN-Akt-Mdm2 network and mainly focused on Hopf bifurcation and oscillation behavior of the system. Firstly, we analyse the role of time delays in the system and find that the time delay could cause Hopf bifurcation. In order to facilitate the research, we let $\tau_2 = 0$ and the threshold value $\tau _1^0 \approx {\rm{0}}{\rm{.70455}}$ for generating the Hopf bifurcation is obtained by calculation. Then, we obtain that the period and amplitude increase theoretically and numerically with the rising of $\tau_1$. Besides, for the first step of numerical simulation, we verify the existence of Hopf bifurcation consistent with theoretical derivation. Next, we take $\tau_2$ into account that $\tau_{1}$ and $\tau_{2}$ and infer that $\tau_1$ and $\tau_2$ can cause oscillation and suppressed oscillation respectively. More important is that we found that the relative size of $\tau_2$ and $\tau_2$ controls the pulse shape of P53. The dynamic properties of the whole system could be subverted for their minor change. Furthermore, we study the influences of four important parameters, namely S(concentration of serum), $k_{3}$(Rate constant for $Mdm2_{p}$ dephosphorylation), $k_{10}$(Basal expression rate of PTEN) and $k_{14}$(Rate constant for PTEN-induced Akt dephosphorylation) in the system. The results indicate that we obtained is the stability and the amplitude of the periodic solution could also be destroyed by their minor changes. We mainly focus on the dynamic properties of P53 and its related components under non-stress conditions, therefore this research might be helpful to further understand the mechanism of P53 and its related components under the normal circumstances and provide a new perspective for preventing cancer.

Acknowledgments

The authors are grateful to the anonymous referees and editors for their valuable suggestions which helped us to improve the quality of this work. This research was funded by the National Natural Science Foundation of China (11762022) and the youth academic and technical leaders of Yunnan Province (2019HB015).

Conflict of interest

The authors declare that they have no conflict of interest.

Supplementary

1.   The existence of Hopf bifurcations

For the convenience of presentation, we set x, y, z, u, v, p and q to represent Mdm2, Mdm2$_{p}$, MA, $M_{p}$A, P53, PTEN and Akt respectively. And parameters of the model are replaced by the form in Table 1, our system will be converted to

Obviously, it is easy to get a positive equilibrium solution by setting the right-hand side of the equation equal to zero in (6.1), which is quite consistent with the actual biological background and significance. For the convenience of research, we set $\tau_{2} = 0$ and $E^{*} = (x^{*}, y^{*}, z^{*}, u^{*}, v^{*}, p^{*}, q^{*})$ as a positive balance point of the system(6.1), and then we can get the linearization equation near the equilibrium as follows

where

Furthermore, we can get the characteristic equations of equations (6.2)

where

In order to theoretically analyze the sufficient conditions for generating oscillations, we assume that $iw(w > 0)$ is a root of equation (6.3). Then, bring $iw$ into equation (6.3) to get:

Separating the real part and the imaginary part we can get

Then there will be

The above equations are squared and then the equations of the equal sign are added separately. Then, there will be

Moreover, it is easy to verify that

Therefore, equation (6.7) exists unique positive root $w_0$. From (6.6) we can get

where j = 0, 1, 2, 3$\cdots$. And we define $\tau _1^0 = \min \left\{ {\tau _1^{(j)} > 0} \right\}_{j = 0}^{ + \infty }$. Thus, when $\tau = \tau_1^{0}$ equation (6.3) has a pair of pure imaginary roots $\pm i{w_0}$. Next we will prove that the root of equation (6.3) has a strict negative real part at ${\tau _1} \in \left[ {0{\rm{ }}, {\rm{ }}\tau _1^0} \right)$, and there are at least two roots containing positive real parts when ${\tau _1} \in (\tau _1^0{\rm{ }}, {\rm{ }}\tau _1^m)$, where $m$ is the minimum positive integer satisfying for $\tau _1^0{\rm{ < }}\tau _1^m$.

Lemma 6.1 ^[55] For the exponential polynomial $P(\lambda {\rm{ }}, {\rm{ }}{e^{ - \lambda {\tau _1}}}, {\rm{ }}{e^{ - \lambda {\tau _2}}} \cdots {e^{ - \lambda {\tau _m}}})$ in the process of $({\tau _1}, {\tau _2}, \cdots, {\tau _m})$ changes, only when $P(\lambda {\rm{ }}, {\rm{ }}{e^{ - \lambda {\tau _1}}}, {\rm{ }}{e^{ - \lambda {\tau _2}}} \cdots {e^{ - \lambda {\tau _m}}})$ has a zero point on the imaginary axis, or when it has a zero point through the imaginary axis, the sum of the zeros in the half-open plane is likely to change. Where

It can be seen from the lemma that it is only necessary to prove that the root of the equation (6.3) has a strict negative real part at ${\tau _1} = 0$. And then as ${\tau _1} = \tau _1^0$, a pair of pure imaginary roots $\pm i{w_0}$ are derived, therefore, for the second half of our results we only need to prove the following formula:${\left. {\frac{{d{\mathop{\rm Re}\nolimits} (\lambda ({\tau _1}))}}{{d{\tau _1}}}} \right|_{{\tau _1} = \tau _1^0}} > 0$. The left and right ends of the second equation of equation (6.5) are derived for A respectively, we can get

Let $A = {m_3} - 3{m_7}w_0^2 + 5{m_{11}}w_0^4 + 7w_0^6$, $B = 2{m_5}{w_0} - 4{m_9}w_0^3 + 6{m_{13}}w_0^5$, ${\rm{C = cos}}{w_0}\tau _1^0 - {\rm{i}}\sin {w_0}\tau _1^0$, $D = {m_2} - {m_6}w_0^2 + {m_{10}}w_0^4 - {m_{14}}w_0^6$, $E = {m_4} - {m_8}w_0^3 + {m_{12}}w_0^5$, $F = {m_4} - 3{m_8}w_0^2 + 5{m_{12}}w_0^4$, $G = 2{m_6}{w_0} - 4{m_{10}}w_0^3 + 6{m_{14}}w_0^5$. Furthermore, ${\left. {{{\left[ {\frac{{d\lambda ({\tau _1})}}{{d{\tau _1}}}} \right]}^{ - 1}}} \right|_{{\tau _1} = \tau _1^0}} = \frac{{\rm{i}}}{{{w_0}}}\left[ {\frac{{ - (A + B) - C(F + G)}}{{C(D + E)}} + \tau _1^0} \right]$. Then by calculating we can get ${\mathop{\rm sgn}} \left\{ {{{\left. {\frac{{d{\mathop{\rm Re}\nolimits} (\lambda ({\tau _1}))}}{{d{\tau _1}}}} \right|}_{{\tau _1} = \tau _1^0}}} \right\} = {\mathop{\rm sgn}} \left\{ {{{\left. {{\mathop{\rm Re}\nolimits} {{\left[ {\frac{{d\lambda ({\tau _1})}}{{d{\tau _1}}}} \right]}^{ - 1}}} \right|}_{{\tau _1} = \tau _1^0}}} \right\} > 0$.

Therefore, we prove that the root of characteristic equation(6.5) crosses the virtual axis from left to right as $\lambda (\tau _1^0) = \pm {\rm{i}}{w_0}$. Obviously, Hopf bifurcation of the positive equilibrium solution occurs at ${\tau _1} = \tau _1^0$. Besides, the system is progressively stable when ${\tau _1} \in [0, \tau _1^0)$ while, if ${\tau _1} > \tau _1^0$ it is unstable. So for our system, $\tau _1^0$ is a vital quantity, it could predict the oscillation of P53 or not.

2.   Properties of Hopf bifurcations

In the previous section, we have obtained that the system (6.2) undergoes a Hopf bifurcation or Turing-Hopf bifurcation at the positive constant $({x^*}, {y^*}, {z^*}, {u^*}, {v^*}, {p^*}, {q^*})$ when ${\tau _1} = \tau _1^0$. In this section, we will study the properties of Hopf bifurcation and the stability of bifurcated periodic solutions by using the normal form theory and central manifold theory due to Hassard, Kazarinoff and Wan ^[56]. Throughout this section, we always assume that system (2.1) undergoes Hopf bifurcation when ${\tau _1} = \tau _1^0$ at the positive equilibrium $({x^*}, {y^*}, {z^*}, {u^*}, {v^*}, {p^*}, {q^*})$, and the corresponding purely imaginary roots of the characteristic equation are $\pm {\rm{i}}{w_0}$. Subsequently, we let $\bar x = x(t) - x^*, \bar y = y(t)-x^*, \bar z = z(t)-z^*, \bar u = u(t)-u^*, \bar v = v(t)-v^*, \bar p = p(t)-p^*, \bar q = q(t)-q^*$. And we still denote $\bar x, \bar y, \bar z, \bar u, \bar v, \bar p, \bar q$ as $x, y, z, u, v, p and q$. Moreover, let ${\tau _1} = \tau _1^0 + \gamma$ normalizing the time delay ${\tau _1}$ by the time-scaling $t \to {t / {{\tau _1}}}$, and the system (2) transform to

where ${c_1}, {c_2}, \cdots, {c_{20}}$ is consistent with the previously defined, and

Let

And define $C = C\left({\left[ { - 1, 0} \right], {R^7}} \right)$, then (6.2) becomes to

where ${L_\gamma }:C \to {R^7}, f:R \times C \to {R^7}$ whose specific form is as follows

where $\phi = {\left({{\phi _1}\left(t \right), {\phi _2}\left(t \right), {\phi _3}\left(t \right), {\phi _4}\left(t \right), {\phi _5}\left(t \right){\phi _6}\left(t \right), {\phi _7}\left(t \right)} \right)^T} \in C$. According to the Riesz representation theorem there exists a $7 \times 7$ matrix function $\eta (\theta, \gamma)$, which is bounded variogram on $\theta \in \left[ { - 1, 0} \right]$ such that ${L_\gamma }\phi = \int_{ - 1}^0 d \eta \left({\theta, \gamma } \right)\phi \left(\theta \right)$ for $\phi \in C\left({\left[ { - 1, 0} \right], {R^7}} \right)$.

In fact, we can choose

where $\delta \left(\theta \right)$ is Dirac delta function. And we define

Obviously, we can transform (6.10) to the following form

where ${U_t}\left(\theta \right) = U\left({t + \theta } \right)$. For $\psi \in {C^1}\left({\left[ { - 1, 0} \right], {{\left({{R^7}} \right)}^*}} \right)$, we define

and

which is a bilinear inner product. Obviously, ${\Lambda _0}$ and $\Lambda _0^*$ are adjoint operators for each other. In addition, $\pm {\rm{i}}{w_0}\tau _1^0$ are the eigenvalues of ${\Lambda _0}$. Therefore, they are also eigenvalues of $\Lambda _0^*$. We let $q\left(\theta \right)$ be the eigenvector of ${\Lambda _0}$ corresponding to ${\rm{i}}{w_0}\tau _1^0$ and ${q^*}\left(s \right)$ be the eigenvector of $\Lambda _0^*$ corresponding to $- {\rm{i}}{w_0}\tau _1^0$, which meet the following conditions

From ${\Lambda _0}q\left(0 \right) = {\rm{i}}{w_0}\tau _1^0q\left(0 \right)$ and $\Lambda _0^*{q^*}\left(0 \right) = - {\rm{i}}{w_0}\tau _1^0{q^*}\left(0 \right)$, it is easy to deduce

From

to make $\left\langle {{q^*}, q} \right\rangle = 1$, we let $\bar G = \frac{1}{{(1 + {v_1}\bar v_1^* + {v_2}\bar v_2^* + {v_3}\bar v_3^* + {v_4}\bar v_4^* + {v_5}\bar v_5^* + {v_6}\bar v_6^*) + \tau _1^0{e^{ - {\rm{i}}{w_0}\tau _1^0}}{c_4}{v_4}}}$. And then $G = \frac{1}{{(1 + {{\bar v}_1}v_1^* + {{\bar v}_2}v_2^* + {{\bar v}_3}v_3^* + {{\bar v}_4}v_4^* + {{\bar v}_5}v_5^* + {{\bar v}_6}v_6^*) + \tau _1^0{e^{{\rm{i}}{w_0}\tau _1^0}}{c_4}{v_4}}}$.

The above ${\bar v_1}$ and $\bar v_1^*$ represent the conjugate plural of ${v_1}$ and $v_1^*$ respectively and the other analogy. Next we will compute the coordinate to describe the center manifold ${C_0}$ at $\gamma = 0$ using the way of Hassard et al ^[56].

We let ${U_t}$ be the solution of (6.15) at at $\gamma = 0$ and define

On the center manifold ${C_0}$, we can regard $W\left({t, \theta } \right)$ as

In fact, $z\left(t \right)$ and $\bar z\left(t \right)$ are local coordinates for center manifold $C_0$ in the direction of ${q^*}$ and ${\bar q^*}$. It is easy to know that $W\left({t, \theta } \right)$ is real only when ${U_t}\left(\theta \right)$ is real. We just consider it is a real solution of (6.10), and then

where

The following formula can be obtained from equations (6.19) and (6.20).

Combined with (6.12), (6.22), (6.25), we have

Comparing the coefficients with equation (6.24), we obtain

Since ${W_{20}}\left(\theta \right)$ and ${W_{11}}\left(\theta \right)$ are unknown in ${g_{21}}$, we will continue to solve for ${W_{20}}\left(\theta \right)$ and ${W_{11}}\left(\theta \right)$.

From (6.15) and (6.21) we have

where

Differentiating formula (6.20) for $t$ we can have

Bring (6.27) into (6.29), we obtain

Compare the coefficients of ${z^2}$ and $z\bar z$ in (6.29) and (6.30) to obtain the following equation

and

where $I$ is unit matrix or unit transformation. From (6.27) we can get

for $\theta \in \left[ { - 1, 0} \right)$. Comparing the corresponding coefficients of (6.28), (6.33) we can get the following equation

From (6.31) we can obtain

According the definition of ${\Lambda _0}$ and $q\left(\theta \right)$ and combining (6.34), then

where ${N_2} = {\left({N_2^{(1)}, N_2^{(2)}, N_2^{(3)}, N_2^{(4)}, N_2^{(5)}, N_2^{(6)}, N_2^{(7)}} \right)^T}$ is also a constant vector. Next we just need to calculate ${N_1}$ and ${N_2}$. From (6.31), (6.32) and the define of ${\Lambda _0}$ the following equations can be exported.

where $\eta \left(\theta \right) = \eta \left({0, \theta } \right)$. From (6.27) we have

Substituting (6.38) and (6.41) into (6.39), we obtain

And then since ${\rm{i}}{w_0}\tau _1^0$ is the eigenvalues of ${\Lambda _0}$, which corresponding eigenvector $q\left(0 \right)$. Furthermore, From (6.43) to

Similarly, we can get

Finally, according to the definition of $\eta \left(\theta \right)$, we can solve ${N_1}$, ${N_2}$ and ${g_{21}}$. Besides we can also get the following values

From the above discussion we can get the following results:

(1) The direction of Hopf bifurcation is determined by ${\mu _2}$ : if ${\mu _2} > 0$ (resp.${\mu _2} < 0$), then the Hopf bifurcation is supercritical (resp. subcritical) and the bifurcating periodic solutions exist for $\tau > \tau _1^0$ (resp.$\tau < \tau _1^0$);

(2) The stability of the bifurcating periodic solutions is depended on ${\beta _2}$, the bifurcating periodic solutions in the center manifold are stable (resp. unstable) for ${\beta _2} < 0$ (resp. ${\beta _2} > 0$);

(3) The period of the bifurcating periodic solutions is determined by ${T_2}$ :the period increases if ${T_2} > 0$ (resp. decreases ${T_2} < 0$).