Modeling the roles of 14-3-3 σ and Wip1 in p53 dynamics and programmed cell death*

Nan Liu(劉楠),Hongli Yang(楊紅麗)and Liangui Yang(楊聯(lián)貴)

School of Mathematical Sciences,Inner Mongolia University,Hohhot 010021,China

Abstract The transcription factor p53,also known as the‘guardian of the genome’and the‘gatekeeper’,encodes programmed cell death(apoptosis)information in its dynamics.The anti-cancer protein 14-3-3 σ and oncogenic protein Wild-type p53-induced phosphatase 1(Wip1)are both in the p53-murine double minute 2(Mdm2)gene network.However,how these two proteins modulate the dynamics of p53-Mdm2 module and further cell apoptosis remains elusive.Therefore,the authors constructed a relatively comprehensive ordinary differential equation model,taking into account the roles of 14-3-3 σ and Wip1 in p53 dynamics and cell apoptosis.By using numerical methods,we found that p53 exhibits abundant dynamics under the common control of 14-3-3 σ and Wip1,and explained the anti-apoptotic effect of Wip1 and the pro-apoptotic effect of 14-3-3 σ.By comparing experimental phenomena,we proposed that the p53 oscillation caused by the participation of Wip1 and 14-3-3 σ is closer to the real scene.This paper provides a new insight for p53 dynamics mechanism and apoptosis outcome mediated by Wip1 and 14-3-3 σ,and may be helpful for clinical application.

Keywords:p53 dynamics,14-3-3 σ,Wip1,apoptosis

1.Introduction

In recent years,both experimental approaches and mathematical models have become the mainstream of studying biological/medical systems[1,2].In general,experimental phenomena and biological knowledge provide the‘raw materials’for the construction of mathematical models,and in turn,mathematical model-generated interpretation and prediction are helpful for biological understanding and experimental design,forming a cycle work flow[3].There are many well-studied biological systems,and one of them is the p53 protein molecular network.As far as we know,the application of the mathematical model about p53 to explain the experimental phenomena can be traced back at least to the beginning of this century,that is,Lev et al used Western blotting to observe the attenuated pulse of p53 and established the corresponding p53-murine double minute 2(Mdm2)vibrator model[4].Subsequently,a wave of research on the p53 oscillation mechanism started.There are many ways to make p53 oscillate,including linear or nonlinear feedback,and regulation with or without time delay[5].Moreover,the types of bifurcations leading to the initiation of p53 oscillations are also diverse[6].Which type of bifurcation causes p53 oscillation is still an open question.The p53 oscillation produced by the supcritical Hopf bifurcation(sup-HB)can explain the experimental findings that the period is more robust than the amplitude[6].However,the p53 oscillation caused by subcritical Hopf bifurcation(sub-HB)or saddle node invariant circle bifurcation(SNIC)also makes sense in some scenarios[7].In addition,the oscillation of p53 tends to induce cell cycle arrest experimentally[8,9].Therefore,whether the oscillation behavior of p53 can be related to the downstream cell cycle arrest gene is also worth considering.

Here,a class of recognized cell cycle arrest proteins 14-3-3 σ attracts our attention[10].Mdm2 is the main negative regulator of p53 in most scenarios.As a positive regulator of p53,14-3-3 σ outstandingly shortens the half-life period of Mdm2 on the one hand and thereby prolongs the half-life period of p53 as expected on the other hand[11].In addition to Mdm2,wild-type p53-induced phosphatase 1(Wip1),another negative regulator of p53,was experimentally found to be a key factor in p53 oscillation[12].As a phosphatase,Wip1 can catalyze the dephosphorylation of p53 and Mdm2,as well as ataxia telangiectasia mutated(ATM),an upstream protein of p53-Mdm2 module[13].From the perspective of carcinogenesis,p53 dysfunction has been discovered in many cancers[14],and the situation in 14-3-3 σ is similar.It has reported that 14-3-3 σ is down-regulated in several cancers,including breast cancer[15],ovary cancer[16],salivary glands adenoid cystic carcinoma[17],hepatocellular carcinoma[18],prostate cancer[19],and lung cancer[20].On the contrary,Wip1 overexpression has been also reported in a variety of human cancers,such as breast cancer[21,22],and neuroblastoma[23].However,the specific mechanisms of Wip1 cancer promotion and 14-3-3 σ cancer suppression are still poorly understood.

Motivated by the above considerations,we developed a mathematical model to qualitatively illustrate the dynamics of p53 and the mechanisms underlying programmed cell death initiation.Firstly,we gave the gene network framework that is concerned by us and the corresponding ordinary differential equation dynamic system,and found a parameter set that can induce this dynamic system to vibrate repeatedly.Then,the pulse dynamics of p53 were explained under two types of cell pressure and the p53 steady levels under the mediation of Wip1 or 14-3-3 σ were tested.The robustness of this p53-Mdm2 vibrator was studied.We emphasized the timeliness of cell stress,and pointed out that the influence of Wip1 on Mdm2 may have a weak effect on p53 oscillations.Subsequently,the oscillating or bistable p53 kinetic modes were revealed under the joint regulation of Wip1 and 14-3-3 σ.We highlight that the p53 vibrator most likely needs both positive and negative feedback and the key players are Wip1 and 14-3-3 σ.Finally,the principle of 14-3-3 σ deletion or Wip1 enrichment inducing cancer cells to evade death fate was demonstrated.Our work can have guiding value for future p53 dynamics research and medical applications.

2.Model and methods

The core of the model is the negative feedback between p53 and Mdm2,that is,p53 promotes the transactivation of Mdm2 gene,but Mdm2 accelerates the degradation rate of p53 via the ubiquitin-proteasome pathway[4].After DNA doublestrand break,ATM quickly senses DNA damage and is activated to a phosphorylated form[12].The activated ATM can further promote self-activation and the phosphorylation of p53 and Mdm2[24–26].For p53,ATM-induced phosphorylation enables it to acquire transcription ability and enhance ubiquitination resistance[24],i.e.the activity and stability of p53 are improved;For Mdm2,ATM-induced phosphorylation makes it more susceptible to self-degradation,and its E3 ubiquitin ligases function is lost[25].Activated p53 can induce the expression of Wip1,which stands on the opposite side of ATM:the proteins phosphorylated by ATM are restored by Wip1 in our model[13].Therefore,Wip1 participates in at least two pathways to reduce the activity of p53:(i)promoting dephosphorylation of p53 is a direct pathway,and(ii)promoting dephosphorylation of ATM is an indirect pathway.In addition,Wip1 is also involved in reducing the stability of p53:promoting the dephosphorylation of Mdm2 leads to the stabilization of Mdm2 and the restoration of Mdm2?s E3 ubiquitinase function,those are not conducive to the stable existence of p53[13,27].14-3-3 σ,another target protein of p53 that is antagonistic to Wip1 by promoting the degradation of Mdm2[11].These proteins form the protein signal network shown in figure 1.It is difficult to intuitively understand the dynamics of p53 without a mathematical model.

Without considering the interaction between different substrates,the basic equations(1)–(8)can be written by using the law of mass action,i.e.the reaction P→P*can be expressed as a function term:k[P],where[P]represents the dimensionless concentration of protein P;P*stands for the protein P with a phosphate group;k is the reaction rate.

?The form of ATM switches between active and inactive,but their total amount is almost unchanged[26].These facts make the kinetic terms of ATM’s generation and degradation are not included in the model,which is the main difference between ATM’s equations and other substances

?Under resting conditions,p53 is rapidly degraded after production[4].In other words,inactive p53 has poor stability.After DNA damage,p53 is activated into a stable phosphorylated form[24].Therefore,the equation of p53 should consist of birth,death,activation,and deactivation dynamic terms The degradation rate of p53 should be greater than that of p53*,so we set

?Interestingly,Mdm2,the main negative regulator of p53,changes exactly the‘opposite’to that of p53 after DNA damage.The‘opposite’here refers to stability and enzyme activity,i.e.the phosphorylated Mdm2 induced by DNA damage is rapidly degraded on the one hand,and the effect of Mdm2*ubiquitin ligase is weakened on the other hand[25].In the kinetic equations,the degradation rate of Mdm2 should be much smaller than that of Mdm2*,i.e.Besides,the kinetic terms in the Mdm2 equations are similar to those of p53.

?Without involving complex situations,the equations of Wip1 and 14-3-3 σ only contain a synthesis term and a degradation term

When the mutual regulation relationship within the protein components in the network shown in figure 1 is introduced into the model,some reaction rate parameters should be replaced by a generalized Hill function[28].That is,the reaction(P catalyzes the conversion of P1into P2)is expressed as a dynamical termin the equation.Here,v is the maximum reaction rate,ρ(0≤ρ≤1)is a proportional constant characterizing the dependence degree of the reaction on the enzyme,j represents the halfsaturated concentration of the enzyme,and n is the Hill coefficient determined by the number of enzyme molecules in polymers.Since enzymes can significantly improve the reaction efficiency,we set ρ≥0.9.Active monomer ATM plays a catalytic role[26],so its Hill coefficient is 1.P53 acts as a transcription factor in a tetramer form[29],thus the Hill coefficient of p53*is 4.It is assumed that the Hill index of Wip1 and 14-3-3 σ are 1 and 2,respectively.The subfunctions of the equations(1)–(8)are as follows equations(9)–(19):

Indeed,there are few valid experimental data that can evaluate parameter values.Many mathematical models are for investigating the qualitative mechanism of a certain topological network rather than the quantitative data fitting[30,31].Our work is also like this.In order to reproduce the experimental phenomenon,the unknown parameters are obtained by the‘trial and error’method.The default parameters are in table 1.The unit of concentration is dimensionless‘C’,and the unit of time is minutes.

The initial values in table 2 are selected as the steady state level of various proteins without DNA damage(in the case of kpatm=0).Numerical analyses are performed by using the softwares OSCILL8 and XPPAUT.In the analysis,the parameter ranges in the one-dimensional bifurcation diagrams were selected the cases that can clearly show the limit cycle.Assume that the minimum rate of a biochemical reaction is 0,and the maximum is about 2 times of the default value.

Table 1.Simulation parameters.

Table 2.Initial condition.

3.Results

3.1.Reproduce the dynamic behaviors of p53

In our theoretical model,DNA damage directly causes ATM activation,which is reflected in the parameter kpatm.The kinetics of p53 is dependent on the DNA damage dose in an environment of permanent stimulation[9].When a chemotherapy drug,Etoposide,is used as the DNA-damage stimulus,only an appropriate drug dose can trigger p53 oscillations.Thereby,we carried out the experiments in silico as shown in figures 2(a)and(b),and found that excessive DNA damage can totally destroy the dynamic behavior of p53 pulses,which is consistent with biological experiments[9].Ionizing irradiation is another stimulus type that can trigger p53 oscillations.Different from the treatment with Etoposide solution,ionizing irradiation is a transient pressure on cells,which results in countable p53 pulses[32].In a hybrid mathematical model of p53,the DNA repair process after radiation is characterized by using Monte Carlo methods:kpatmis almost fixed during DNA repair and the time spent by DNA repair is roughly proportional to the radiation dose[33].Therefore,we assume that the cell is exposed to radiation at t=500 min,and it takes 900 min to complete DNA repair.Simulation results are shown in figure 2(c),the p53-Mdm2 module generates three pulses in this case,which is in agreement with the role of p53-Mdm2 feedback loop as a‘digital’clock[32].Furthermore,the experimental phenomenons that DNA damage induces Mdm2?s destabilization and the pulse of Mdm2 follows the pulse of p53 also appeared in our numerical simulations[32,34].

Additionally,it has been experimentally found that the level of p53 induction decreased in Wip1 overexpression cells;in contrast,p53 can reach a high level in the cells treated with Wip1 RNA inhibitor during the DNA damage period[12].In figure 3,we reproduced these two points by adjusting the Wip1 synthesis rate parameter swip1.14-3-3 σ positively regulates p53 level[11].In agreement with its role,figure 3 indicates that the absence of 14-3-3 σ may lead p53 to approach the basal level,while the enrichment of 14-3-3 σ can create a high steady state of p53.One might be curious,which of the three arms of Wip1 inhibiting p53 activity and stability is dominant in mediating p53 levels?The answer given by figure 3 is that,the lower level of the p53 steady state is determined by the Wip1-ATM arm or the Wip1-p53 arm,while the higher level of the p53 steady state is determined by the Wip1-Mdm2 arm.When the Wip1-regulated dephosphorylation rate of ATM,p53,or Mdm2 doubles,the lower level of p53 appears in the first two scenarios.And when the biochemical reaction of Wip1-regulated dephosphorylation of ATM,p53,or Mdm2 is removed,the highest level of p53 occurs in the last one scenario.Note that the increase in kqmdm2does not affect the p53 oscillation,we will study this issue in the next section.

We have explained the feasibility of the model.To our knowledge,oscillation is dependent on the topology of the gene network,the form of regulation function,and parameter selection.Hence,it is necessary to test the sensitivity of oscillation to parameters.In our theoretical model,the oscillation is confined to the parameter interval between two HBs in figure 4.In the second column of table 3,the value range of the parameter that can make p53 oscillate is given as HB1–HB2.It should be noted that when HB is not in the range of physical significance,take the maximum/minimum value with physical significance,that is,for ρ,HB2 on the right side of the oscillation interval is not greater than 1,and for all parameters,HB1 on the left side of the oscillation interval is not less than 0.If only the size of the parameter range where the oscillation occurs is compared,the robustness conclusion obtained is unreliable.For example,the acceptable interval length of ATMt is much larger than kpatm.However,compared with kpatm,ATMt is inherently large,so it is not accurate to claim that ATMt is more robust than kpatmfor p53 oscillations.It seems more reasonable to compare the percentage ranges in which the default parameters can be

3.2.Robustness of oscillation to parameters variation

Table 3.Robustness analysis of p53 oscillations.

increased(Id%)or decreased(Dd%).From this view,ATMt is a sensitive parameter,compared with kpatm.In the third column of table 3,the Dd% of kpatmis bigger than that of ATMt.Interestingly,the Dd%values corresponding to ρ6and j6are 100%,which means that the oscillation of p53 is completely independent of Wip1?s regulation on Mdm2.Indeed,Wip1?s role in Mdm2 dephosphorylation is ignored in many p53 oscillator models[35,36].

In addition,the ranges in which some parameters can be increased or decreased is extremely unbalanced.For example,the range of ρ6can be increased or decreased to 2.04% or 100%,respectively.It is still not very convenient to study the global robustness.Here,we define S=(HB2?HB1)/(HB2+HB1)to describe the sensitivity of oscillation behavior to parameters.Obviously,the ideal(most robust)parameter value is(HB1+HB2)/2,and half the length of the oscillation interval is(HB2-HB1)/2.The larger the S value corresponding to a certain parameter is,the less sensitive the oscillation is to the relative variation of this parameter under the ideal condition.Our calculation results are shown in table 3 fifth column.Approximately 66% of the parameters have an S value of more than 10%,indicating the robustness of this oscillator to the parameter set.

3.3.Dynamics of p53-Mdm2 module regulated by 14-3-3 σ and Wip1

Figure 5 exhibits the dynamic behaviors of the system under the common regulation of Wip1 and 14-3-3 σ.The parameter plane stretched by the production rates of these two proteins is roughly divided into three types of regions by the two types of bifurcation curves in figure 5(a).Here,roughness mainly means that only the Fold bifurcation and the Hopf bifurcation of the steady-state branch are displayed,while the bifurcations of the limit cycle oscillation branch are ignored,such as the Fold bifurcation of the limit cycle oscillation.We do this omission because the bifurcation curves of the limit cycle oscillation branch are too close to the Hopf bifurcation curve.In the white parameter area,the system will eventually converge to a unique steady state,regardless of the initial state.As shown in the pseudo-phase plane of figure 5(b),there is only one stable fixed point.In the blue parameter area,the system can finally converge to two steady states,and the initial state determines the final steady state of the system.As shown in the pseudo-phase plane of figure 5(c),there are two stable fixed points.In the yellow parameter area,the system will eventually have lasting and stable limit cycle oscillation dynamics,regardless of the initial state.As shown in the pseudo-phase plane of figure 5(d),there is a stable limit cycle.Interestingly,the Fold bifurcation represented by the black line seems to have no contribution to the division of the dynamic regions.In fact,the Fold bifurcation is the fundamental reason why the system appears bistable,but the Hopf bifurcation within the two Fold bifurcation points reduces the bistable region.Rigorously,this dynamic system also has the excitable state in the white parameter area,but it beyond our research and is not considered.Clearly,in this theoretical model,the p53 oscillation is generated by the interconnection between Wip1 mediated negative feedback and 14-3-3 σ mediated positive feedback.This is a typical characteristic of relaxation oscillations.

Figure 1.Schematic diagram of the network model.Promotion and state transition is represented by colored and black arrow-headed lines respectively.The promotion of degradation is indicated by cycle-headed lines.

Figure 2.Time evolution diagram of total p53(blue line)and total Mdm2(red line).(a)kpatm=0.25;(b)kpatm=0.5;(c)kpatm=0.25 if t?[500,1400]min,else,kpatm=0.

Figure 3.Time evolution diagrams of total p53([p53]+[p53*])in some special situations.

Figure 4.Codimensional 1 bifurcation diagram.HB represents the Hopf bifurcation point;the red line is the stable equilibrium state;the black line is the unstable equilibrium state.The green dots are the upper and lower bounds of the stable limit cycles;the blue cycles are the upper and lower bounds of the unstable limit cycles.

Figure 5.(a)Codimension 2 bifurcation diagram.The green line represents the Hopf bifurcation,and the black line represents the Fold bifurcation.The blue area indicates the bistable state;the white area indicates the monostable state;the yellow area indicates the oscillation.(b)–(d)are pseudo-phase plane,black dots represent stable fixed points,and black closed orbit represents stable limit cycle.

Figure 6.Pseudo-phase plane analysis.The blue line and red line represent the pseudo-collar slope lines.The gray line indicates orbit.Black dot represents unstable fixed point.

Figure 7.Codimension 1 bifurcation diagram.sup-or sub-HB represents the supcritical or subcritical Hopf bifurcation point respectively;the red line is the stable equilibrium state;the black line is the unstable equilibrium state.The green dots are the upper and lower bounds of the stable limit cycles;the blue cycles are the upper and lower bounds of the unstable limit cycles.

Figure 8.The schematic diagram of ECP,the blue solid line is the time history of p53*,the blue dashed line is the threshold Th,and the gray area is the accumulation of p53 that is conducive to apoptosis.

Figure 9.The functional relationship between ECP and kqatm(green inverted triangles),kqp53(cyan upright triangles),kqmdm2(black diamonds),swip1(red squares),and s14-3-3σ(pink five-pointed stars).

Figure 10.ECP as a function of swip1 and s14-3-3σ.The colors indicate the value of ECP,as the color bar on the right.The red arrow points to the direction of canceration.

In order to further explain the mechanism of relaxation oscillation,we adopt the pseudo-phase plane analysis method.Firstly,let the right side of equation(7)be zero,that is,[Wip1]degenerates into a parameter;by making a bifurcation curve of[14-3-3 σ]versus[Wip1],a Z-shaped pseudo-nullcline(zero slope line)can be obtained as shown in figure 6 blue solid line.Secondly,let the right side of equation(8)be zero to get another pseudo-nullcline.The two nullclines intersect in the middle branch of the Z-shaped line,so the fixed point is unstable.After DNA damage,the concentration of 14-3-3 σ rise rapidly,and the level of Wip1 enhance follow the increase of the concentration of 14-3-3 σ;when the concentration of Wip1 exceed the right knee of the Z-shaped solid line,the level of 14-3-3 σ begin to decrease rapidly until nears the lower branch of the Z-shaped solid line;meanwhile,the concentration of Wip1 will decrease to a level below the left knee of the Z-shaped solid line,resulting in the 14-3-3 σ level rapidly enhance once again.The increase and decrease process of this cycle can create the oscillations of p53.Mathematically,the Folded(Z-shaped)bifurcation curve cannot be separated from the positive feedback.Thereby,besides Wip1-mediated negative feedback,a positive feedback is also required for the relaxation oscillations of p53.Indeed,there are experimental evidences that preliminarily prove that the oscillation of p53 requires both positive and negative feedback.As a negative regulator of p53,Wip1 can significantly affect the occurrence of p53 oscillations has been confirmed experimentally by Batchelor et al[12];recently,exciting experiments have found that a type of negative regulator of Mdm2(miR-192,a transactivation target of p53,similar in function to 14-3-3 σ,can reduce Mdm2 level)also has a constructive effect on the p53 vibrator[37].As aforementioned,there are many types of bifurcations that can give rise to p53 oscillations[6,7].In this model,p53 oscillations can be born from both sup-HB and sub-HB in figure 7.Therefore,there is a great possibility that the oscillation of p53 is a relaxation oscillation,just like our model.

3.4.Impact of Wip1 and 14-3-3 σ expression levels on the programmed cell death

There has been much theoretical work on revealing the relationship between p53 dynamics and cell apoptosis.In the original hypothesis,it was believed that p53 needs to reach a high steady state to initiate apoptosis[38],which is based on the biological knowledge that the binding affinity of p53 to pro-apoptotic genes is weaker than that of pro-survival genes[39].Latter,two research groups from Nanjing University in China proposed the hypothesis that p53 pulses induce apoptosis:in the model of Sun et al,pro-apoptotic proteins slowly accumulate in p53 pulses until the threshold is reached[40];in the model of Zhang et al,the post-translational modification of p53 gradually changed from pro-survival to proapoptotic[41,42].In these models,apoptosis seems to require p53 oscillations peak to be high enough,and the number of p53 pulses to be enough[40–42].In general,the high steady level of p53 induces faster apoptosis,while the excessive number of p53 pulses causes slower apoptosis.Recently,Wu et al proposed an apoptosis hypothesis,which can explain both the rapid and the slow induction of apoptosis by p53 high steady state and excessive pulses,respectively[43].For simplicity,this paper adopts the hypothesis of Wu et al.

Regardless the dynamic patterns of p53,in the default initial condition,as the level of Wip1 increases,the value of ECP will decrease(red squares in figure 9),and when swip1increases from 1.3 times of the default value to 1.4 times,there is a big gap in ECP,like a switch.It should be noted that the‘switch’behavior in here may also change continuously,or equivalent to the so-called‘ultrasensitive’[44].As mentioned above,Wip1 regulates the p53-Mdm2 module through three pathways.Among them,the pathways that directly or indirectly inactivate p53 have a more obvious impact on ECP(figure 9 green and cyan triangles).Our model shows that,like the kinetics of p53(the indexes of robustness of kqmdm2are all relatively high),ECP is not sensitive to the parameter kqmdm2(figure 9 black diamonds),which means that Wip1?s promotion on Mdm2 stability is not very effective in blocking apoptosis.Contrary to Wip1,as the level of 14-3-3 σ increases,the value of ECP will increase(figure 9 pink fivepointed stars),and when s14-3-3σincreases from 0.7 times of the default value to 0.8 times,the ECP‘switch’is‘toggled open’,and the apoptosis probability is greatly increased.In order to show the mechanism of canceration more intuitively,we make a three-dimensional map as shown in figure 10.The red arrow indicates the direction of canceration.It can be seen that Wip1 level is too high or 14-3-3 σ level is too low will not be conducive to cell apoptosis.Our results are well in line with the tumor suppressor function of 14-3-3 σ and the tumor engine function of Wip1.

4.Conclusion and discussion

In summary,this article uses the ordinary differential equation dynamics model to study the physiological phenomena of dynamics in the p53 protein signal network and apoptosis underlying p53 dynamics.Our innovation lies in exploring the simultaneous effects of 14-3-3 σ and Wip1 on the p53-Mdm2 module.We explained the reason of p53 oscillation,that is,the relaxation oscillation caused by the interconnection of 14-3-3 σ mediated positive and Wip1 mediated negative feedbacks,and explained the roles of Wip1 and 14-3-3 σ in apoptosis,that is,the accumulation of active p53 exceeding the threshold within a certain period of time will be downregulated by Wip1 and up-regulated by 14-3-3 σ.Tumor cells often have strong proliferation ability but weak apoptosis ability.We explained the tumor-promoting function of Wip1 and the tumor-inhibiting function of 14-3-3 σ from the perspective of apoptosis.As for cancer suppression strategies,targeting Wip1 axis and 14-3-3 σ axis at the same time may be a more optimized approach.

Note that in the experiments of Wu et al[43],Bax is activated when apoptosis occurs,which is consistent with the model prediction of Sun et al[40];the concentration of p53 increases significantly when apoptosis is initiated,which is consistent with the model prediction of Zhang et al[42].From the aspect that the p53 concentration burst during the apoptotic phase,the model of Wee et al is also reasonable[38].In the mitochondrial apoptotic pathway,Bax inserts into the outer mitochondrial membrane,enhances the permeability of the mitochondrial membrane,releases cytochrome C,and then forms a positive feedback with caspase 3 to amplify the apoptotic signal cascade[45,46].In the simplified mathematical models,a certain step in these processes can be selected as an apoptotic marker[40–42].In our view,‘ECP’can be a key indicator in these models of crosstalk between p53 module and apoptotic module if Th is selected to a proper value.We think that it is also meaningful work to establish more computational models for the network motifs downstream of p53,similar to[31].Moreover,whether the similarswitch behaviors in the figure 9 are continuous is an open question for the readers.

Undoubtedly,this paper has some significance for understanding the roles of Wip1 and 14-3-3 σ in the p53 protein signaling network,but there are also some shortcomings at the same time.For model building,a partial differential equation model with spatial distribution can be considered[47],a stochastic differential equation model with internal or external noise can be considered[35],and so on.For cellular results,there are actually more options besides apoptosis,such as DNA repair guided by p53 dynamics[40,41],cell cycle arrest guided by p53 dynamics[41,42],senescence guided by p53 dynamics[8],etc.In addition,due to the high dimensionality of the model,only numerical calculations have been applied.In the future,we will reduce the dimensionality of the model to facilitate the application of more mathematical theoretical methods,such as the linear stability analysis and bifurcation theory[48],to better understand the relaxation oscillation of p53.We hope that this work can provide some enlightenment to researchers in related disciplines.

Acknowledgments

The authors are very grateful for the professional and detailed suggestions of the anonymous reviewers.