Evolutionary game dynamics of cooperation in prisoner's dilemma with time delay

1.   Introduction

In mammalian cells, the tumor suppressor P53 is activated under various cell pressures and ensures an appropriate response from arrest and repair to the induction of senescence and apoptosis ^[1]. One of the main defense mechanisms against cellular stress is the activation of the tumor suppressor P53 ^[2]. However, alterations in the structure and function of the P53 network can lead to serious human diseases, such as cancer ^[3]. P53 itself is regulated by a feedback loop: in the absence of DNA damage, it is ubiquitinated by the E3 ligase Mdm2 and degraded by the proteasome ^[1,4,5]. Upon the DNA damage, post-translational modifications of P53 and Mdm2 prevent their interaction and allow P53 to aggregate in the nucleus, where it binds to the target gene promoters and induces their expression ^[1]. At the moment, a particularly dangerous type of damage is DNA double-strand breaks (DSBs), which can destroy the integrity of the genome and can have harmful consequences if not repaired ^[6]. In order to combat these injuries, cells will evolve sensitive sensing mechanisms that activate the DNA damage response (DDR) and induce transient cell cycle arrest in the G1 or G2 phase or terminal cell fates, such as apoptosis and senescence ^[6,7,8].

The tumor suppressor P53 is highly regulated by a complex network of interactions whose dominant characteristics vary by cell stress, cell type and species ^[9]. Recent studies have shown that the specific function of P53 is also related to its expression level. In normal mammalian cells, P53 is usually maintained at a low level due to the down-regulation of Mdm2, which helps cells avoid premature aging and apoptosis ^[10,11]. The concentration of P53 changes to a medium-level state when receiving a weak stimulus, thereby inducing cell cycle arrest and avoiding the inheritance of abnormal DNA replication to offspring ^[12,13]. When the stimulation was eliminated, the expression level of P53 returned to the normal level ^[14,15]. P53 will have a sustained oscillations state when stimulated by mild or moderate stimulation, which will contribute to cell cycle arrest. However, the concentration of P53 reaches a higher level when subjected to greater stimulation, which in turn induces cell apoptosis ^[16,17].

Time delays generally exist in various systems. From the perspective of dynamic systems, a time delay can cause larger oscillations, which plays an important role in biological systems or gene regulatory networks. It can produce a wealth of dynamic behaviors, such as periodic dynamics and chaotic dynamics ^{[18,19,20,21]}. In general, there is an average transcriptional time delay of 10$\sim$20 minutes between the action of a transcription factor on a gene promoter and the presence of the corresponding mature mRNA in the cytoplasm ^[22]. Similarly, the synthesis of a typical protein from mRNA requires a translation time delay of 1$\sim$3 minutes. In addition, due to the low copy number of the molecule and the fluctuation of the environment ^[23,24], the interference of noise cannot be avoided in the system. Noise can induce bistability, oscillation, and bifurcation, and it plays a key role in a series of processes including skeleton dynamics, cellular polarization, signal transduction and neural activity. From a dynamic point of view, bistability refers to the phenomenon that there exist two stable equilibrium points under a certain set of parameters of a dynamic system. Different initial values will cause the solution curves to tend to different stable equilibrium points. From a biological point of view, bistability is a common phenomenon of cell behavior, mainly showing the coexistence of two phenotypes that can occur under certain conditions. In order to adapt to the harsh living environment, some organisms often choose appropriate forms of expression according to the changes in the environment. For example, Ozbudak et al. discussed the condition of bistability behavior of the lac operon, which is a gene involved in lactose metabolism in Escherichia coli and its expression level can determine the state of the cells ^[25]. They found that in the absence of glucose, the lac operon gene was not induced (low expression levels) at low thio-methylgalactoside (TMG) concentrations (< 3 ${\rm{ \mathsf{ μ}}}$M), but could be induced (high expression level) at high TMG concentrations (> 3 ${\rm{ \mathsf{ μ}}}$M). If starting from the uninduced cells (low TMG concentration), the TMG concentration is steadily raised until the lac operon gene is induced in the cells after more than 30 ${\rm{ \mathsf{ μ}}}$M. Conversely, starting from the induced cells (high TMG concentration), the concentration of TMG was gradually reduced until it was lower than 3 ${\rm{ \mathsf{ μ}}}$M to turn off the expression of the lac operon gene. Hence, when the concentration of TMG is between two thresholds, the bacteria exhibit bistable dynamics.

Based on the above biological background, this paper proposes an improved model by including noise and time delay and studies the dynamic characteristics, including stability and bifurcation. The bifurcation analysis shows that under the appropriate parameter range P53 oscillation can be induced. Since time delay is common in various systems ^[20,21], we study the performance of the system by analyzing the influence of time delay. Taking time delay as the bifurcation parameter, the stability of the system and the existence of Hopf bifurcation are studied by using Hopf bifurcation theory. It is found that time delay is not only an important condition for the generation of Hopf bifurcation, but that it also affects the period and amplitude of the Hopf bifurcation. At the same time, the research also shows that the combination of time delay can not only promote the oscillation of the system, but it also has better robustness, which means that the parameter domain in which oscillation occurs has been broadened due to time delay. In order to study the more abundant dynamic characteristics of the system, the impact of time delay and parameter interaction on its stability is also investigated. In addition, because of the low copy number of molecules and the fluctuation of the environment ^[24], we also consider the impact of noise on the system. It is found that the noise can change the state of the system without delay when the stochastic differential equation is used for numerical simulation. Meanwhile, under the action of a single time delay or multiple time delays, the noise can not only induce the switching of the system state, but it can also promote system oscillation. The above research may help us to further understand the regulation mechanism of P53 in the cell cycle.

2.   Materials and methods

In the absence of pressure, ataxia telangiectasia mutated protein (ATM) will form inactive homodimers. After being induced by DSBs, ATM is rapidly activated by the complex formed by Mre11, Rad50 and NBS1 and decomposed into catalytically active monomers through autophosphorylation ^[26]. The activated ATM (ATM$^*$) phosphorylates the histone variant H2AX, which acts as a scaffold to recruit more proteins related to DDR and damage repair ^[7,27]. The negative feedback loop of Mdm2 and Wip1 is an important regulatory element in the P53 network in the model proposed by ^[7]; the corresponding interaction diagram is shown in Figure 1. P53-Mdm2 is the core negative feedback loop in response to the DNA DSB ^[7,28]. As a transcription factor, P53 can induce the transcription of Mdm2 ^[29], whereas normally functioning Mdm2 can be used as an E3-ubiquitin ligase to target P53 for proteasomal-mediated degradation ^[4,5]. Another feedback loop is P53-Wip1-ATM. The transcription of Wip1 is induced by P53. The phosphatase Wip1 not only directly dephosphorylates ATM$^*$, but it also inhibits ATM$^*$ self-activation ^[7,28,30].

Figure 1.  Interaction diagram of the coupling feedback loops of the P53-Wip1-ATM network and P53-Mdm2 network. The figure includes phosphorylation (P), dephosphorylation (dP) and ubiquitination (Ub). The mRNA (mdm2 and wip1) species are shown in purple, and all other species represent proteins. Black arrows indicate stimulus or inhibition; the red arrows represent DNA damage; the cyan arrows indicate noise effects. Due to the existence of time delay in various systems, we chose to study the performance of the system by analyzing the influence of time delay. In addition, due to the randomness of DNA DSBs, we also considered the influence of noise in the system.

DownLoad: Full-Size Img PowerPoint

However, many works on P53 oscillation have supposed that transcriptional and translational events are instantaneous, whereas the time needed in these gene expression processes has been poorly noticed. It is well known that gene expression refers to the process of synthesizing genetic information from a gene to a functional gene product that is usually a protein or functional RNA such as miRNA and sRNA. Transcription, transcript splicing and processing and translation are essential steps in this synthesis process. Meanwhile, the existing theoretical results suggest that the time delay can usually lead to oscillation and complex dynamic behaviors of gene regulatory networks. This shows that the delay of protein synthesis is a common phenomenon, and it is reasonable and necessary to include the delay factor in the model of the gene regulation network.

In this paper, we introduced four time delays to the model proposed by M$\ddot{o}$nke et al. ^[7] to study the performance of the P53 gene regulation model. The delayed differential equations are as follows:

The variables involved in the model are ATM$^*$, P53, mdm2, Mdm2, wip1 and Wip1. $ATM^*$ represents the activated ATM. P53 is the tumor suppressor P53 protein; mdm2 and Mdm2 are mdm2 mRNA and mdm2 protein, respectively; wip1 and Wip1 are wip1 mRNA and wip1 protein, respectively. Here, $\tau _1$ is the time delay needed in the dephosphorylation of ATM$^*$ depending on the transport of Wip1. $\tau _2$ is the time delay needed in the degradation of P53 depending on the transport of Mdm2. $\tau _3$ is the time delay of Mdm2 protein synthesis. $\tau _4$ is the time delay of Wip1 protein synthesis. Other parameters of this model were mostly selected according to the literature^[7], and all parameter values are shown in Table 1.

Table 1.  Parameter values for the mathematical model.

Parameter	Description	Value	Reference
$A$	Maximum self-activation rate of ATM	$30.5$	[7]
$P$	Dephosphorylation rate of ATM by Wip1	$22$	[28]
$C$	Production rate of P53	$1.4$	[30]
$g$	maximum degradation of P53 by Mdm2	$2.5$	[30]
$d_{AM}$	Degradation rate of Mdm2 by ATM$^{*}$	$20$	[28]
$T_{m}$	Mmaximum production rate of mdm2	$1$	[7]
$T_{M}$	Maximum production rate of Mdm2	$4$	[7]
$T_{\omega}$	Maximum production rate of wip1	$1$	[7]
$T_{W}$	Maximum production rate of Wip1	$1$	[7]
$d_{A}$	Basic dephosphorylation rate of ATM$^{*}$	$0.16$	[7]
$d_{P}$	Basic degradation rate of P53	$0.1$	[31]
$d_{m}$	Basic degradation rate of mdm2	$1$	[7]
$d_{M}$	Basic degradation rate of Mdm2	$2$	[30]
$d_{\omega}$	Basic degradation rate of wip1	$1.3$	[7]
$d_{W}$	Basic degradation rate of Wip1	$2.3$	[28]
$k_{A}$	Michaelis constant for the ATM$^{*}$ self activation	$0.5$	[7]
$k_{WA}$	Michaelis constant for the inhibition of the ATM$^{*}$ self-activation by Wip1	$0.14$	[7]
$k_{MP}$	Michaelis constant of the degradation for P53 by Mdm2	$0.15$	[30]
$k_{Pm}$	Michaelis constant for the production of mdm2 by P53	$1$	[7]
$k_{P{\omega}}$	Michaelis constant for the production of wip1 by P53	$1$	[7]
$R$	Strength of the inhibition of the Mdm2 mediated degradation of P53 by ATM$^{*}$	$2$	[7]
$\gamma$	Michaelis contant for the signal strength of the DSB process	$9$	[32]
$S_{max}$	Maximum signal strength of the DSB process	$0.2$	[33]
${\rm DSB}$	DNA damage in the form of DSBs	$0.1$	Estimate
$\tau_{1}$	Time delay needed in the dephosphorylation of ATM$^*$ depend on the transport of Wip1	0$\sim$5	Estimate
$\tau_{2}$	Time delay needed in the degradation of P53 depend on the transport of Mdm2	0$\sim$5	Estimate
$\tau_{3}$	Time delay needed in the Mdm2 protein synthesis	0$\sim$5	Estimate
$\tau_{4}$	Time delay needed in the Wip1 protein synthesis	0$\sim$5	Estimate

---

 | Show Table

DownLoad: CSV

3.   Results

In this research, particular attention has been paid to the bistability and oscillations of the P53-Mdm2-Wip1 network. Bistability refers to the ability to achieve two different internal states under different stimulus conditions. It is ubiquitous in natural biomolecular networks and has important biological significance ^[34,35]. Cells can switch their state between two internal states to accommodate environmental and intercellular conditions owing to regulatory interactions among cellular components. Oscillatory behaviors are produced by the corresponding interactions among genes, proteins, and metabolites. They are used to control various aspects of cell physiology, including cell growth, maturity, division and death ^[36,37]. Hopf bifurcation, including supercritical and subcritical Hopf bifurcation, is the key principle in the design of biochemical oscillators. Noise is ubiquitous in cells, and it is generated from the fluctuation of the chemical reaction during the transcription and translation process or caused by the fluctuation of the concentration or state of other components in the cell ^[38]. Moreover, research has shown that noise can often cause bistability, oscillations and bifurcations, even if these behaviors cannot be observed in deterministic models ^[13,14]. In addition, noise-induced stochastic resonance is widely used in chemical analysis, instrument analysis and signal detection^[16,17]. Next, we will investigate these dynamics via numerical simulation and bifurcation analysis by using XPPAUT and Matlab software. To be specific, first of all, the influence of important parameters in the model on the system when there is no time delay is investigated. This part essentially considers the model proposed by M$\ddot{o}$nke et al. ^[7]. Second, the influence of a single time delay on system dynamics is studied. Third, the cooperative effect of multiple time delays is discussed. Furthermore, the cooperative effects of important parameters and multiple time delays are explored. Finally, the influence of noise on the system is researched.

3.1.   Influence of model parameters on the system

In this section, the model of M$\ddot{o}$nke et al. (let $\tau_{1} = \tau_{2} = \tau_{3} = \tau_{4} = 0$ in Eq (2.1)) is adopted to clarify the influence of different model parameters on the dynamic behavior of the P53 system. Here, we mainly study the influence of six parameters, including $P$, $d_{m}$, $T_{\omega}$, $C$, $T_{M}$ and $T_{W}$ that respectively refer to the dephosphorylation rate of ATM by Wip1, the basic degradation rate of mdm2, the maximum production rate of wip1, the production rate of P53, the maximum production rate of Mdm2 and the maximum production rate of Wip1. It should be noted that in the bifurcation diagram below, the red and black lines respectively represent stable and unstable steady states, green dots are the maxima and minima of the oscillation state corresponding to the stable limit cycles and the blue open circles denote the maxima and minima of the oscillation state corresponding to the unstable limit cycles.

Experiments have shown that the activation of ATM will induce sustained oscillation of P53. Next, we will study how the level of ATM affects the expression level of P53. For this reason, we draw the bifurcation diagram of P53 concentration versus $P$ which is the dephosphorylation rate of ATM by Wip1, as shown in Figure 2(A). When the dephosphorylation rate of ATM by Wip1 is very small, the concentration of P53 is at a high level, which promotes cell apoptosis. With the increase of the dephosphorylation rate of ATM by Wip1, the P53 system appears bistable approximately in the interval $(4.708, 13.13)$ and a subcritical Hopf bifurcation point appears at $P\approx13.13$, which leads to a series of sustained oscillations in the P53 system. These oscillations disappear when the dephosphorylation rate of ATM by Wip1 increases and exceeds the bistable region. Finally, the concentration of P53 enters a low state, which helps cells avoid premature senescence and apoptosis ^[20,28,39].

Figure 2.  Bifurcation diagrams of P53 with $P$ and $d_{m}$. (A) Bifurcation diagram of P53 with $P$. (B) Bifurcation diagram of P53 with $d_{m}$.

DownLoad: Full-Size Img PowerPoint

It can be seen in Figure 2(B) that the expression level of P53 is dependent on $d_{m}$ which is the basic degradation rate of mdm2. Specifically, when the basic degradation rate of mdm2 is small enough, there exists a unique stable steady state. With the basic degradation rate of mdm2 increasing gradually, the bistability of the system appears approximately in the interval $(0.1062, 0.2931)$ and a Hopf bifurcation point appears at $d_{m}\approx0.2931$, which leads to sustained oscillations of the P53 system. Interestingly, in the bistable region, the high steady state and the low steady state coexist. However, the concentration of P53 enters a low state when the basic degradation rate of mdm2 increases and exceeds the bistable region, which means that the cell is in a dormant state.

As can be seen in Figure 3(A), when $T_{\omega}$, i.e., the maximal production rate of wip1, is low, the concentration of p53 is in a high state. As the maximal production rate of wip1 increases gradually, bistability occurs in the range starting from $T_{\omega}\approx0.2609$ to $T_{\omega}\approx0.7303$ and a Hopf bifurcation point appears at $T_{\omega}\approx0.7303$, which leads to a series of sustained oscillations of the system. When the maximal production rate of wip1 increases to a certain extent, the concentration of P53 is in a low stable steady state, which means that the cells enter a normal state. A bifurcation point (HB) appears when the maximal production rate of wip1 takes an appropriate value, which indicates that the maximal production rate of wip1 is very important for the dynamics of P53.

Figure 3.  Bifurcation diagrams of P53 with $T_{\omega}$ and $C$. (A) Bifurcation diagram of P53 with $T_{\omega}$. (B) Bifurcation diagram of P53 with $C$.

DownLoad: Full-Size Img PowerPoint

As observed in Figure 3(B), the concentration of P53 is located at a high state if $C$ (the production rate of P53) is small enough. As the production rate of P53 increases gradually, bistability occurs approximately in the interval $(0.1549, 0.7123)$. In particular, a Hopf bifurcation point occurs at $C\approx0.7123$, which means that the system generates sustained oscillations. When the production rate of P53 increases to a certain extent, the concentration of P53 enters a low state, which corresponds to the normal cellular state.

$T_{M}$ and $T_{W}$ denoted the maximum productivity of Mdm2 and the maximum productivity of Wip1, respectively. Because the concentration of Wip1 and Mdm2 can affect the concentration of P53, we also studied the influence of $T_{M}$ and $T_{W}$ on the P53 dynamics, which may help to further understand the significance of the model parameters in the dynamic behavior of P53. To show this, the bifurcation diagrams of the P53 concentration versus $T_{M}$ and $T_{W}$ are given in Figure 4.

Figure 4.  Bifurcation diagrams of P53 with $T_{M}$ and $T_{W}$. (A) Bifurcation diagram of P53 with $T_{M}$. (B) Bifurcation diagram of P53 with $T_{W}$.

DownLoad: Full-Size Img PowerPoint

Figure 4(A) reveals that when the value of the maximum production rate of Mdm2 changes, the concentration of P53 also changed. It can be seen that there exists a unique stable steady state if the maximum production rate of Mdm2 is small enough. With the maximum production rate of Mdm2 increasing gradually, the system appears bistable approximately in the interval $(12.01, 37.79)$. In this range, the high steady state and the low steady state coexist. The concentration of P53 enters a high state when the maximum production rate of Mdm2 rises and exceeds the bistable region, which means that the cell undergoes apoptosis. Figure 4(B) shows that with the increase of $T_{W}$, there is a certain change in the concentration level of P53. When the maximum production rate of Wip1 is low, the bistability of the P53 concentration occurs. Then, as the increase of the maximum production rate of Wip1 further exceeds the bistable region, the concentration of P53 enters a low stable state, which indicates that the cells eventually enter a dormant state. Through numerical simulation, it is found that the proper parameter values can cause system oscillation. Therefore, we can change the dynamic properties of the system by regulating important parameters, and then we can control the dynamic characters of P53.

3.2.   Influence of time delay $\tau_{2}$ on system stability

Negative feedback is the most fundamental requirement for the generation of biological oscillations. However, for a network with one node, a negative feedback loop alone can not lead to biological oscillations unless it is coupled with other elements including time delays. The negative feedback loop among a network with two or more nodes can generate oscillation under appropriate parameter conditions, but if time delay is incorporated, it can enhance the robustness of the period and amplitude of the oscillation ^{[40,41,42,43]}. Here, we study the performance of the system by analyzing the effect of time delay. We have theoretically deduced that time delay is an important condition to induce Hopf bifurcation. Because this process is very tedious, for the convenience of reading, it has been placed in the supplementary material, and only the numerical analysis is given in the main text. In order to study the effect of time delay $\tau_{2}$ on the system, we assume that the values of time delay $\tau_{1}$, $\tau_{3}$ and $\tau_{4}$ are zero. It is easy to verify that the system has a unique positive equilibrium $E^*{\approx}(0.00164, 0.19733,$ $0.16481, 0.32428, 0.12678, 0.05512)$.

Figure 5 displays the time course of the P53 system. Clearly, the positive equilibrium point is stable if there are no time delays in the system, which can be seen in Figure 5(A). If there is a time delay in the system, its influence on P53 is shown as in Figure 5(B)–(D). The positive equilibrium point $E^{*}$ of the system is asymptotically stable when $\tau_{2} = 1$ and unstable for $\tau_{2} = 2$. This result indicates that as the time delay increases and crosses through the critical value, the equilibrium of the system loses its stability and the system undergoes a Hopf bifurcation. The critical value of time delay for the Hopf bifurcation was calculated to be $\tau_{2}^{0}\; {\approx}\; 1.20245$.

Figure 5.  Dynamic behavior of the system (2.1) under the influence of the time delay $\tau_{2}$ when $\tau_{1} = \tau_{3} = \tau_{4} = 0$. (A) The positive equilibrium point $E^{*}$ of the system is stable with the time delay $\tau_{2} = 0$. (B) The positive equilibrium point $E^{*}$ of the system is asymptotically stable when the time delay $\tau_{2} = 1$ $(\tau_{2} < \tau_{2}^{0})$. (C) and (D) The positive equilibrium point $E^{*}$ of the system is unstable when the time delay $\tau_{2} = 2$ $(\tau_{2} > \tau_{2}^{0})$.

DownLoad: Full-Size Img PowerPoint

Figure 6(A) is the time evolution diagram of the system under the action of time delay $\tau_{1}$. Obviously, the positive equilibrium point $E^{*}$ of the system is stable if $\tau_{1} = 0$. And as the time delay $\tau_{1}$ increases, the system is still in a stable state. In Figure 6(B), we can see that the system generates a series of sustained oscillations when the value of the delay $\tau_{2}$ is greater than the critical value, and that the period and amplitude of the oscillations increase with the increase of $\tau_{2}$. It can be seen in Figure 6(C) that the positive equilibrium point $E^{*}$ of the system is stable when the time delay $\tau_{3} = 0$. The positive equilibrium point of the system is asymptotically stable when the time delay $\tau_{3} = 1$. However, the positive equilibrium $E^{*}$ of the system (2.1) is unstable when the time delay $\tau_{3} = 1.5$ or $\tau_{3} = 2$. Figure 6(D) shows that the system is in a stable state no matter whether it is under the effect of the time delay $\tau_{4}$ or without the time delay $\tau_{4}$. This indicates that, for the parameter values shown in Table 1, the time delays $\tau_{2}$ and $\tau_{3}$ can induce system oscillation alone, while the time delays $\tau_{1}$ and $\tau_{4}$ can not do this.

Figure 6.  Dynamic behavior of P53 with the time delays $\tau_{1}$, $\tau_{2}$, $\tau_{3}$ and $\tau_{4}$. (A) The effect of the time delay $\tau_{1}$ on system dynamic behavior when $\tau_{2} = \tau_{3} = \tau_{4} = 0$. (B) The effect of the time delay $\tau_{2}$ on system dynamic behavior when $\tau_{1} = \tau_{3} = \tau_{4} = 0$. (C) The effect of the time delay $\tau_{3}$ on system dynamic behavior when $\tau_{1} = \tau_{2} = \tau_{4} = 0$. (D) The effect of the time delay $\tau_{4}$ on system dynamic behavior when $\tau_{1} = \tau_{2} = \tau_{3} = 0$.

DownLoad: Full-Size Img PowerPoint

3.3.   Influence of multiple time delays on system stability

In the previous subsection, we separately investigated the effect of the time delays $\tau_{1}$, $\tau_{2}$, $\tau_{3}$ and $\tau_{4}$ on the dynamic behavior of the system. We found that, with the selected parameter values, adjusting the time delays $\tau_{2}$ and $\tau_{3}$ not only induces system oscillations, but also changes the amplitude and period of the oscillations. However, adjusting the time delays $\tau_{1}$ and $\tau_{4}$ in the selected parameter values did not induce system oscillations. In order to better study the influence of time delay on the system, we will study the influence of the combination of multiple time delays on the dynamic behavior of the system.

As shown in Figure 7(A), (C), the system is in a sustained oscillatory state when only the action of $\tau_{2}$ is considered. However, when considering the influence of the time delays $\tau_{1}$ and $\tau_{4}$, the state of the system does not change significantly, which indicates that the system is not sensitive to the time delays $\tau_{1}$ and $\tau_{4}$. As shown in Figure 7(B), the system is in a low amplitude oscillation state when only the action of the time delay $\tau_{2}$ is considered. But the system is in an oscillation state with a higher amplitude when the effect of the time delay $\tau_{3}$ is considered, and the amplitude increases with the increase of $\tau_{3}$, which indicates that the time delay $\tau_{3}$ plays a role in promoting the oscillation of the system. Figure 8 can further verify the above conclusion. In general, the relative magnitude of these four time delays may have an effect on the stability, period and amplitude of the system.

Figure 7.  Dynamic behavior of the system (2.1) under the combinatorial effect of the time delay $\tau_{2}$ and one of other time delays when $\tau_{2} = 2$. (A) $\tau_{1}$ is taken as different values and $\tau_{3} = \tau_{4} = 0$. Changing the time delay $\tau_{1}$ does not change the state of the system. (B) $\tau_{3}$ is taken as different values and $\tau_{1} = \tau_{4} = 0$. The amplitude of the system oscillation increases with the time delay $\tau_{3}$. (C) $\tau_{4}$ is taken as different values and $\tau_{1} = \tau_{3} = 0$. Changing the time delay $\tau_{4}$ does not change the state of the system.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Comparison of the dynamic behavior of the system (2.1) under the influence of a single time delay and multiple time delays. (A) The effects of the time delay $\tau_1$ and the interaction between the time delay $\tau_1$ and the other three time delays on the system. The system is in a stable state when only time delay $\tau_1$ acts, while oscillation appears when the rest time delays are considered. (B) The effects of the time delay $\tau_2$ and the interaction between the time delay $\tau_2$ and the remaining three time delays on the system. When there is only the time delay $\tau_2$ acting on the system, the oscillation amplitude is low, while the amplitude becomes larger when the other time delays are considered. (C) The effects of the time delay $\tau_3$ and the interaction between time delay $\tau_3$ and the rest of the three time delays on the system. The system is in a low amplitude oscillatory state when only the delay $\tau_3$ acts, while it is enlarged after considering the action of the remaining time delays. (D) The effects of the time delay $\tau_4$ and the interaction between the time delay $\tau_4$ and the remaining of three time delays on the system. The system is in a stable state when only the time delay $\tau_4$ acts, while oscillation appears when the rest of the time delays are considered.

DownLoad: Full-Size Img PowerPoint

3.4.   Combinatorial regulation of time delay and parameter interaction on system stability

In the previous study, we found that the dephosphorylation rate of ATM by Wip1, the basic degradation rate of mdm2, the maximum production rate of wip1, the production rate of P53, the maximum production rate of Mdm2 and the maximum production rate of Wip1 could cause changes in the P53 concentration. At the same time, it is also mentioned in literature^[16,17] that the parameters in the mathematical model have a great influence on system dynamics. In order to further understand the dynamics of the P53 system, we will study the dynamic behavior under the combinatorial effects of parameters and time delays. To achieve this, we fixed the time delay and selected different parameter values for numerical simulation.

As shown in Figure 9(A), with the increase of $C$, i.e., the production rate of P53, the system is always in a stable state. As can be seen in Figure 9(B), the system changes from a stable state to an oscillating state with the increase of the production rate of P53, which means that the system state can be changed by changing the value of the production rate of P53. Figure 9(C) reveals that when $C = 0.2$, the system is in a stable state. With the increase of $C$, the state of the system changes from a stable state to a state of sustained oscillation, and its period and amplitude increase with the production rate of P53. In Figure 9(D), it can be seen that the system is always in a stable state when the production rate of P53 increases.

Figure 9.  Time evolution diagram of P53 in the system (2.1) under the combinatorial effect of the production rate of P53 $C$ and time delays. (A) The time history diagram of P53 with the combinatorial effect of $C$ and time delay $\tau_{1}$ ($\tau_{1}$ = 3) when $\tau_{2} = \tau_{3} = \tau_{4} = 0$. (B) The time history diagram of P53 with the combinatorial effect of $C$ and the time delay $\tau_{2}$ ($\tau_{2}$ = 2) when $\tau_{1} = \tau_{3} = \tau_{4} = 0$. (C) The time history diagram of P53 with the combinatorial effect of $C$ and the time delay $\tau_{3}$ ($\tau_{3}$ = 2) when $\tau_{1} = \tau_{2} = \tau_{4} = 0$. (D) The time history diagram of P53 with the combinatorial effect of $C$ and the time delay $\tau_{4}$ ($\tau_{4}$ = 3) when $\tau_{1} = \tau_{2} = \tau_{3} = 0$.

DownLoad: Full-Size Img PowerPoint

In Figure 10(A), (D), we can see that with the increase of $d_{m}$, i.e., the basic degradation rate of mdm2, the system lies in a stable state. However, Figure 10(B), (C) indicate that as the basic degradation rate of mdm2 increases, the state of the system transitions from a steady state to an oscillating state. It can be seen in Figure 10 that the four time delays have some differences in regulating the oscillation dynamics of the system.

Figure 10.  Time evolution diagram of P53 concentration under the combinatorial regulation of the basic degradation rate of mdm2 $d_{m}$ and time delays. (A) The time history diagram of P53 concentration under the combinatorial regulation of the time delay $\tau_{1}$ ($\tau_{1}$ = 3) and $d_{m}$ when $\tau_{2} = \tau_{3} = \tau_{4} = 0$. (B) The time history diagram of P53 concentration under the combinatorial regulation of time delay $\tau_{2}$ ($\tau_{2}$ = 2) and $d_{m}$ when $\tau_{1} = \tau_{3} = \tau_{4} = 0$. (C) The time history diagram of P53 concentration under the combinatorial regulation of the time delay $\tau_{3}$ ($\tau_{3}$ = 2) and $d_{m}$ when $\tau_{1} = \tau_{2} = \tau_{4} = 0$. (D) The time history diagram of P53 concentration under the combinatorial regulation of the time delay $\tau_{4}$ ($\tau_{4}$ = 3) and $d_{m}$ when $\tau_{1} = \tau_{2} = \tau_{3} = 0$.

DownLoad: Full-Size Img PowerPoint

Figure 11 plots the time evolution diagram of P53 concentration when the time delay is combined with $T_{M}$, i.e., the maximum production rate of Mdm2. As shown in Figure 11(A), (D), the state of the system did not change significantly with the increase of $T_{M}$. However, in Figure 11(B), (C), we can see that as the maximum production rate of Mdm2 increases, the system changes from an oscillation state to a stable state. Therefore, we can regulate the state of the system by using parameters and time delay.

Figure 11.  Time evolution diagram of P53 concentration under the combinatorial effect of time delays and $T_{M}$, which that is the maximum production rate of Mdm2. (A) The time history diagram of P53 concentration with the combinatorial effect of $T_{M}$ and $\tau_{1}$ ($\tau_{1}$ = 3) when $\tau_{2} = \tau_{3} = \tau_{4} = 0$. (B) The time history diagram of P53 concentration with the combinatorial effect of $T_{M}$ and $\tau_{2}$ ($\tau_{2}$ = 2) when $\tau_{1} = \tau_{3} = \tau_{4} = 0$. (C) The time history diagram of P53 concentration with the combinatorial effect of $T_{M}$ and $\tau_{3}$ ($\tau_{3}$ = 2) when $\tau_{1} = \tau_{2} = \tau_{4} = 0$. (D) The time history diagram of P53 concentration with the combinatorial effect of $T_{M}$ and $\tau_{4}$ ($\tau_{4}$ = 3) when $\tau_{1} = \tau_{2} = \tau_{3} = 0$.

DownLoad: Full-Size Img PowerPoint

3.5.   Influence of noise in the system

Physical, chemical and biological systems may be disturbed by the noise generated by the random birth and death of a single molecule or the fluctuation of biochemical reactions ^{[23,38,44,45]}. It can not only cause bistability but also oscillation and bifurcation ^[34,46]. These behaviors play an important role in biological systems ^[24,34,46]. In order to better understand the dynamics of the P53 model, we introduce random perturbations into the deterministic model to reveal richer and more complex dynamics. Various researchers have used a binomial $\tau$-leap algorithm, Monte Carlo and other methods to study the properties of noise; it was found that the noise will affect the state of the system^{[35,47,48,49,50]}. Here, we established stochastic models of P53-Wip1 and P53-Mdm2 networks and use the Wiener process to deal with the noise term, which focuses on the impact of noise and time delay on the system.

According to the randomness of DSB dynamics^[7,51], we consider the influence of external noise on the interaction between ATM and P53; then, the system becomes as follows:

where $\xi(t)$ is Gaussian white noise, which has the properties $E[\xi(t)] = 0$ and $E[\xi(t)\xi(t+\iota)] = 2\sigma\delta(\iota)$ and $\sigma$ is the noise intensity. Here, the influence of noise on the system is represented by the Wiener process. For the noise term, we mainly focus on noise intensity.

From the previous numerical simulation, it is known that all states of the molecules in the studied system are synchronous, which refers to the qualitative properties including oscillation, and that the steady state of the six molecules are concordant. In this paper, we mainly focus on the switch between oscillation and steady state. Therefore, to study if oscillations of the system occur, one only needs to discuss the oscillation dynamics of one of six molecules. In the following study, we only choose the variable P53 for numerical simulation. The time evolution of P53 under three different noise intensities is shown in Figure 12(A). We have found that the noise can not only cause the oscillation of the system but also induce the state switching without time delay. When the noise intensity $\sigma = 0$, the concentration gradually approaches a certain stable value ^[24]. At a noise intensity $\sigma = 0.02$, the system has a slight oscillation. Interestingly, oscillation with a large amplitude occurs when the noise intensity $\sigma = 0.12$.

Figure 12.  Influence of noise intensity on the system when $\tau_{1} = \tau_{2} = \tau_{3} = \tau_{4} = 0$. (A) Time history diagram of P53 concentration in the system (3.1) under different noise intensities. (B) Phase diagram of the variables P53 and Mdm2 in the system (3.1) under the influence of different noise intensities.

DownLoad: Full-Size Img PowerPoint

Figure 12(B) shows the phase plane diagram of the system (3.1) without time delay. When the noise intensity $\sigma = 0$, the phase trajectory is represented by the magenta curve and the system is in a steady state. When the noise intensity $\sigma = 0.02$, the motion track of the system is represented by a black curve. At this time, the track has obvious deviation, which was obviously caused by noise. When the noise intensity $\sigma = 0.12$, the phase trajectory of the system is represented by the cyan curve, which shows that the system is in an oscillation state. Thus, Figure 12(B) further visually captured the phenomenon in Figure 12(A).

In order to further understand the dynamics of the P53 model, we will study the dynamics of P53 under the combined action of time delay and noise to reveal richer and more complex dynamics of P53. The time history diagram of P53 under the influence of a single time delay and noise is shown in Figure 13. Figure 13(A), (D) shows that when the time delay fails to induce system oscillation, adding noise with positive intensity to the system can cause the system to oscillate. In Figure 13(B), (C), we can see that the system is in a state of low amplitude oscillation under the action of a single time delay, while the system is in a state of high amplitude oscillation when noise with a positive intensity is added to the system. This indicates that noise can promote the oscillation of the system.

Figure 13.  Time history diagram of P53 concentration in the system (3.1) under the action of a single time delay and noise. (A) $\tau_1$ = 3, $\tau_2 = \tau_3 = \tau_4 = 0$. (B) $\tau_2$ = 2, $\tau_1 = \tau_3 = \tau_4 = 0$. (C) $\tau_3$ = 2, $\tau_1 = \tau_2 = \tau_4 = 0$. (D) $\tau_4$ = 3, $\tau_1 = \tau_2 = \tau_3 = 0$.

DownLoad: Full-Size Img PowerPoint

The time evolution process of P53 under the influence of multiple time delays and noise is shown in Figures 14–16. It can be seen in Figures 14(A), (B) and (D) that when $\sigma = 0$, the system is in a low-amplitude sustained oscillation state. While $\sigma = 0.12$, the system is in a state of sustained oscillation with high amplitude. This means that noise can enhance the oscillation of the system. Figure 14(C) shows that the system is in a stable state when $\sigma = 0$, while the system is in an oscillation state when $\sigma = 0.12$, which indicates that the noise can cause the system to oscillate.

Figure 14.  Time history diagram of P53 concentration in the system (3.1) under the action of noise and multiple time delays. (A) $\tau_1$ = 3, $\tau_2$ = 2, $\tau_3 = \tau_4 = 0$. (B) $\tau_1$ = 3, $\tau_3$ = 2, $\tau_2 = \tau_4 = 0$. (C) $\tau_1$ = 3, $\tau_4$ = 3, $\tau_2 = \tau_3 = 0$. (D) $\tau_2$ = 2, $\tau_3$ = 2, $\tau_1 = \tau_4 = 0$.

DownLoad: Full-Size Img PowerPoint

Figure 15.  Time history diagram of P53 concentration in the system (3.1) under the action of noise and multiple time delays. (A) $\tau_2$ = 2, $\tau_4$ = 3, $\tau_1 = \tau_3 = 0$. (B) $\tau_1$ = 3, $\tau_2$ = 2, $\tau_3$ = 2, $\tau_4 = 0$. (C) $\tau_1$ = 3, $\tau_2$ = 2, $\tau_4$ = 3, $\tau_3 = 0$. (D) $\tau_1$ = 3, $\tau_3$ = 2, $\tau_4$ = 3, $\tau_2 = 0$.

DownLoad: Full-Size Img PowerPoint

Figure 16.  Time history diagram of P53 concentration in the system (3.1) under the action of noise and multiple time delays. (A) $\tau_2$ = 2, $\tau_3$ = 2, $\tau_4$ = 3, $\tau_1 = 0$. (B) $\tau_1$ = 3, $\tau_2$ = 2, $\tau_3$ = 2, $\tau_4$ = 3.

DownLoad: Full-Size Img PowerPoint

From Figures 15 and 16, it is found that when $\sigma = 0$, the system is in a low-amplitude oscillation state, and when the noise intensity increases to $\sigma = 0.12$, the system is in a high-amplitude oscillation state, which indicates that the noise plays a role in promoting the oscillation under the action of multiple time delays. The above numerical results show that the introduction of noise into the system can not only cause system oscillation, but it can also enlarge its amplitude. The effect of noise on the dynamics of the P53 system has been studied in several studies ^[47,48,49]. In ^[47], the random behavior of the P53-Mdm2 network is described by using the Langevin equation; also, both the intrinsic noise and the extrinsic noise are considered. The Monte Carlo simulation results show that the two kinds of noises can cause system oscillation. In ^[48], random simulations using the binomial $\tau$-leap algorithm showed that noise is critical to system oscillation. In ^[49], it is found that the noise will induce sustained oscillations with random excursions in the larger cycle corresponding to spikes. In ^[50], it is pointed out that the mechanism of oscillation induced by noise maybe occurs via a slowly varying Ornstein–Uhlenbeck process and three numerical examples are showed to explain it. The above articles introduce the methods of dealing with noise from different perspectives, and the focus of this paper is to pay attention to the impact of noise and time delay on the system. More specifically, time delay and noise can induce oscillation alone under certain parameter ranges. Although the emphasis is different, the results of these papers are consistent with that of this paper.

4.   Conclusions

In recent years, a large number of studies have pointed out that P53 is mutated in cancer cells, which affects the fate of cancer cells. The P53 in normal cells can induce the repair or apoptosis of damaged cells, thereby inhibiting cancer, while the mutated P53 does not have such a function, which will lead to the occurrence of cancer. The important and interesting function of P53 is achieved by regulating the cell cycle process ^[1]. Due to the inherent complexity of biological systems and the importance of P53 in cell fate decisions, it is particularly important to consider the dynamic behavior of the transcription factor P53 under the influence of time delay and noise. In particular, M$\ddot{o}$nke et al. ^[7] combined quantitative single cell data with methods from dynamical systems theory and proposed a network architecture and the corresponding mathematical model based on positive and negative feedback loops that reproduced experimentally measured P53 dynamics over a wide range of conditions.

In this paper, four time delays were introduced to the model proposed by M$\ddot{o}$nke et al. ^[7]. Hopf bifurcation and numerical simulation of the model were used to study the dynamic characteristics of the P53 system under the influence of time delay and noise, including stability and bifurcation. First, we performed bifurcation analysis on the parameters $P$, $d_{m}$, $T_{\omega}$, $C$, $T_{W}$ and $T_{M}$ in the original model proposed by M$\ddot{o}$nke et al. ^[7], and the results showed that important parameters could induce P53 oscillation in an appropriate range. Second, the time delay was found to be ubiquitous in all kinds of systems. We studied the performance of the P53 system by analyzing the influence of time delay. By employing Hopf bifurcation theory, the stability of the system and the existence of Hopf bifurcation were studied by using time delay as the bifurcation parameter. It was found that time delay plays a key role in inducing Hopf bifurcation and adjusting the period and amplitude of oscillation. At the same time, it was also found that the combination of time delays can not only promote the oscillation of the system, but it also offers better robustness.

Then, we studied the influence of important parameters and time delay on the dynamic behavior of the system. The results showed that the P53 system changed from a stable state to an oscillating state with the increase of the production rate of P53 and the basal degradation rate of mdm2 under a certain time delay, which leads to cell cycle arrest and the repair of DNA damage. However, with the increase of the maximal production rate of Mdm2, the P53 system changes from an oscillating state to a stable state, which indicates that appropriately changing the parameter values can change the system state. In addition, due to the low molecular copy number and the fluctuation of the environment, we also considered the influence of noise on the system. Through numerical simulation, we found that the noise can change the state of the system without time delay. Under the action of a single or multiple time delays, the noise can not only change the state of the system, but it can also enlarge its amplitude.

In summary, $\tau_2$ and $\tau_3$ can induce P53 oscillations alone, while $\tau_1$ and $\tau_4$ cannot, but they can fine-tune the oscillations. $\tau_2$ refers to the time delay required for P53 degradation that is dependent on the transport of Mdm2, $\tau_3$ refers to the time delay required for Mdm2 protein synthesis, $\tau_1$ refers to the time delay required for dephosphorylation of the ATM$^*$ that is dependent on the transport of Wip1 and $\tau_4$ refers to the time delay required for Wip1 protein synthesis. This suggests that the P53-mdm2 core negative feedback loop plays an essential role in the induction of P53 oscillations, while the P53-Wip1-ATM negative feedback loop plays an auxiliary role ^[52,53,54]. Second, in the absence of time delay, noise can also cause P53 oscillations, which reveals that, when the DNA suffers a DSB, intracellular molecular movement is violent and noise is strengthened, thereby activating the P53 signaling pathway and initiating the program of repairing DNA damage ^[55,56]. In addition, important parameters involved in the system, such as the basic production rate of P53, the basic degradation rate of mdm2 RNA, and the maximum transcription rate of mdm2 RNA, can change the critical value of the Hopf bifurcation induced by time delay. This suggests that important parameters, time delay and noise, cooperate to influence the oscillatory behavior of the P53 signaling pathway. Biological facts show that Pp53 oscillations contribute to the repair of damaged cells, while high levels of P53 promote the apoptosis of irreparable cells or cancer cells, avoiding the transmission of damage information to the next generation of daughter cells. Therefore, we can modulate the dynamics of the P53 system by modulating these important parameters, time delays and noise, resulting in cell cycle arrest or apoptosis and thereby preventing tumor development.

Acknowledgments

The authors express gratitude to the anonymous referee for his/her helpful suggestions and the partial supports of the National Natural Science Foundation of China (12162032/12062027), the Scholarship for Academic Leader (2019HB015), Ten Thousand Talent Plans for Young Top-notch Talents, the Science and Technology Innovation Team about Modern Applied Mathematics and Systems Biology of Yunnan Education Department, and Applied Basic Research Foundation (LS21005) of Yunnan Province.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix

The characteristic equation of System (2.1) expanded at the equilibrium point $E = (x_1^{*}, x_2^{*}, x_3^{*}, x_4^{*}, x_5^{*}, x_6^{*})$ is

where

and

Assuming that $\pm \omega i$ is the root of Eq (A.1), we can get

Separating the real and imaginary parts, we can get

Through calculation, the following equations are obtained

where

Noticing the formula $\cos({\omega}{\tau_2})^2+\sin({\omega}{\tau_2})^2 = 1$, we can get

where

Assume $h_{5} > 0$, $h_{4} > 0$, $h_{3} > 0$, $h_{2} > 0$ and $h_{1} > 0$. If $h_{0} < 0$, then Eq (7.5) has a unique positive root $\omega_{0}$, and the corresponding critical value of time delay is

Therefore, when $\tau_{2} = \tau_{2}^{0}$, Eq (A.1) has a pair of pure imaginary roots $i\omega_{0}$, and all other roots have negative real parts.

In addition, let $\lambda(\tau_2) = \upsilon(\tau_2)+i\omega(\tau_2)$ be the root of Eq (A.1), and let it satisfy $\upsilon(\tau_{2}^{0}) = 0$ and $\omega(\tau_{2}^{0}) = \omega(0)$, we can get

By Eq (A.1), the derivative of $\lambda$ with respect to $\tau_2$ can be obtained

So there is

where

Given that $h_{5} > 0$, $h_{4} > 0$, $h_{3} > 0$, $h_{2} > 0$ and $h_{1} > 0$, $H^{'}(\omega_{0}^{2}) > 0$. At the same time, we have

The results show that the P53-Mdm2-Wip1 network is asymptotically stable for $0{\leq}\tau_2 < \tau_{2}^{0}$, but unstable for $\tau_2 > \tau_{2}^{0}$ at the equilibrium point $E^{*}$. So, the system undergoes Hopf bifurcation when the time delay $\tau_2 = \tau_{2}^{0}$. Therefore, it is very important to judge the critical value $\tau_{2}^{0}$ of the P53-Mdm2-Wip1 network.