Mathematical model identifies effective P53 accumulation with target gene binding affinity in DNA damage response for cell fate decision

ABSTRACT

Functional p53 signaling is essential for appropriate responses to diverse stimuli. P53 dynamics employs the information from the stimulus leading to selective gene expression and cell fate decision. However, the decoding mechanism of p53 dynamics under DNA damage challenge remains poorly understood. Here we mathematically modeled the recently dual-phase p53 dynamics under doxorubicin treatment. We found that p53 could perform sequential pulses followed by a high-amplitude terminal pulse at relatively low doxorubicin treatment, whereas p53 became steadily accumulated when damage level was high. The effective p53 integral above a threshold but not the absolute accumulation of p53 precisely discriminated survival and death. Silencing negative regulators in p53 network might promote the occurrence of terminal pulse. Furthermore, lower binding affinity and degradation rate of p53 target genes could favorably discriminate high and low dose doxorubicin treatment. Grouping by temporal profiles suggested that the p53 dynamics rather than the doxorubicin doses could better discriminate cellular outcomes and confer less variation for effective p53 integral reemphasizing the importance of p53 level regulation. Our model has established a theoretical framework that p53 dynamics can work cooperatively with its binding affinity to target genes leading to cell fate choice, providing new clues of optimized clinical design by manipulating p53 dynamics.

KEYWORDS:

• P53 terminal pulse
• DNA damage
• cell fate decision

Introduction

Many transcription factors can respond to external or internal signals and initiate specific downstream programs. The specificity between the input and output may be obtained through multiple strategies such as regulations at transcriptional, translational and posttranslational levels and spatial locations [1–3]. Another layer of regulation derives from the temporal dynamics of transcription factors in which signaling molecules may decode response specificity into rich dynamic features [1]. A notable example is the tumor suppressor p53, which can manifest itself in pulsatile or sustained activation to dictate cell fate [4]. The encoding mechanisms of the complex dynamic features can be ascribed to intricate coordination of feedback loops [5]. Understanding the underlying mechanisms which detect the time evolution of transcription factors and decode them into differential gene expression patterns and cell fate remains a daunting task.

The tumor suppressor p53 acts as a central node in sensing various stimuli and initiates appropriate responses ranging from cell cycle arrest, senescence and apoptosis [6]. The DNA double strand breaks (DSBs) elicited by either irradiation or DNA damaging drugs (e.g. doxorubicin, etoposide) represent a lethal type of signals which can faithfully activate p53 through several checkpoint kinases (e.g. ataxia-telangiectasia-mutated, ATM) [7]. There are a number of positive and negative auto-regulatory feedback loops which actively participate in the regulation of p53 network [8–10]. Consequently, the p53 dynamics at single cell levels under DNA damage exhibit complex patterns ranging from spontaneous or synchronized pulses to sustained elevation [4,11–14]. In a recently report, Wu et al. identified a novel p53 dynamics in response to DNA damage in which a prolonged and mild challenge using genotoxic drugs can evoke sequential p53 pulses followed by a terminal pulse with high amplitude to enact apoptosis [13]. Moreover, an acute and high-dose treatment will trigger a monotonic increase of p53 with similar rate of apoptosis and decreased lifespan [13]. They argued that p53 should exceed a threshold with sufficient accumulation to trigger apoptotic programs and they termed this novel pattern as effective cumulative p53 [13]. However, the exact mechanism underlying how effective cumulative p53 is dynamically generated and manipulates cell fate specifically through theoretical perspectives still remains elusive.

In current work, we investigated how p53 dynamics determines cell fate by constructing a simplified p53 model using ordinary differential equations (ODEs). We first faithfully reproduced the complex patterns of p53 responses with doxorubicin (Dox) treatment under either wild type or genetically manipulated conditions using stochastic simulations. Generally, p53 can either display a finite number of pulses with or without the occurrence of a terminal pulse or a sustained activation depending on the doses of genetic drugs. Then, we identified a potential mechanism which unified p53 dynamics, binding affinity and turnover rates of target genes to dictate cell fate. We further postulated that dynamic features of p53 signaling rather than the amount of applied genetic drugs can encode phenotypic responses (e.g. survival and death). In addition, p53 dynamic features per se can dictate cell fate with markedly reduced variation. Our findings characterized the intricate interplay between effective accumulative p53 and cell fate using mechanistic models and may provide clues to optimized clinical strategies for chemotherapeutic drugs.

Results

Construction of a dynamic p53 model identifies the terminal pulses

To investigate the dynamic p53 responses under doxorubicin (Dox) treatment, we constructed a mathematical model which primarily incorporated two negative p53 feedback loops and one implicit positive feedback loop [15] (Figure 1(a), Table S1 and S2). To match our model to experimental data, we varied the stimulus strength by changing the parameter Dox in silico, we found that cells performed sequential pulses and followed by a terminal pulse (as defined in Wu et al.’s work [13]) with markedly high amplitude (Figure 1(b)). The elevation of total p53 was increased with higher Dox levels (Figure 1(b)). Notably, the amplitude and duration of pulses did not vary significantly under relatively low Dox treatment (Figure 1(b), we in silico set Dox = 0.01, 0.05 or 0.1, the same doses as exemplified in Wu et al.’s work [13]). Note that in current work, “Dox treatment” corresponds to the situation where we in silico changed the Dox parameter. For example, Dox = 0.01 in our model corresponds to 0.01 μM doxorubicin treatment). However, the maximal p53 level and accumulation rate of effective p53 integral (defined as effective p53 integral divided by time duration) increased monotonically with all Dox doses during the simulation within 40 hrs (Figure 1(c,d)). We next performed local sensitivity analysis and found that all kinetic parameters affected effective p53 integral more profoundly (Figure 1(e)). Moreover, pulse period was more robust to parametric change than amplitude (Figure 1(e), compare green and red bars). Specifically, p53 phosphorylation by activated ATM (k₁₂) and translation rate of MDM2p-bound p53 mRNA (k_complex) positively influenced signaling amplitude and effective p53 accumulation, whereas MDM2-mediated p53 degradation (k₁₁) and Wip1-induced ATM* deactivation (k₂₄) may efficiently terminate p53 activation (Figure 1(e), only relatively sensitive parameters were labeled). We compared the correlation among all sensitivity measures and found that only the local sensitivity of effective p53 integral and pulse height was highly correlated (Figure 1(f)). These deterministic data identified the sensitive parameters and terminal pulses in p53 signaling.

Figure 1. A mathematical p53 model identifies p53 terminal pulses. (a) Schematic representation of p53 regulatory network for our model. The numbers correspond to kinetic parameters in Table S2. The parameters k_f and k_r represent association and dissociation respectively between p53 mRNA and phosphorylated MDM2. (b) Concentration of total p53 level at various doxorubicin doses (Dox, from 0.01 to 1). (c) Maximal p53 levels under increasing Dox levels. (d) Time-averaged effective p53 integral (i.e. rate of effective p53 increase in units of μM) against increasing Dox levels was plotted. The end of the simulation is 40 h. (e) (Relative) local sensitivity coefficient of kinetic parameters for effective p53 integral (violet), p53 pulse height (red) and period (green). The most sensitive parameters were marked. The local sensitivity is in units of 1. (f) Correlation between sensitivity measures was shown. R denotes Spearman correlation coefficient. ‘Ep53ʹ denotes effective p53 integral.

Simulation of dynamic behavior mediated by effective p53 integral

The deterministic model did not consider the process of doxorubicin clearance as well as the biological noises, therefore we performed stochastic simulations. Simulation results showed that a fraction of cells only showed discrete number of pulses suggesting that Dox might be totally removed or cleared below the threshold when Dox was set to 0.1 (Figure 2(a), left). However, cells (‘cells’ in in silico simulation indicated a stochastic representation during simulation thereafter) may display sequential pulses accompanied by a high-amplitude terminal pulse at low or medium Dox challenge (Figure 2(a), left and middle, FigureS1). Acute high dose Dox could even lead to sustained activation of p53 (Figure 2(a), right, FigureS1). The fraction of cells with terminal pulses also increased with elevated Dox stimulation (Figure 2(b)). By sorting cells by cell fate, we noticed that the effective p53 integral in apoptotic cells was significantly higher than that in survival cells (Figure 2(c) and Fig. S2). Instead, the lifespan and terminal pulse height were not varied markedly across different doses of Dox (Figure 2(d,e)). The effective p53 integral of apoptotic cells, however, was not significantly changed with varying Dox values (Figure 2(f)). The pulse height and duration were conserved when cells were challenged with relatively low doses of Dox (Figure 2(g)). Furthermore, the dynamic properties above (i.e. statistical comparison) were not affected by locally changing the minimum effective level (data not shown). These results were qualitatively consistent with experimental findings [13] and validated our stochastic model.

Figure 2. Dynamic features of dual-phase p53 under stimulation by different Dox dosage. (a) Selected total p53 trajectories when Dox = 0.01, 0.1 and 10. Note that the unit of ‘p53 levels’ in stochastic simulation is 1. (b) Fraction of terminal pulses at different Dox levels. (c) Comparison of effective p53 integral between survival and apoptotic cells. The unit for effective p53 integral during stochastic simulation is ‘min’ thereafter. (d-e) Lifespan (d) and terminal pulse height (e) treated by relatively low levels (0.01, 0.05 and 0.1) of Dox. Data were represented as mean±SD. (f) Effective p53 integral of apoptotic cells under different Dox levels. (g) Pulse height (top) and duration (bottom) for cells with Dox = 0.01, 0.05 and 0.1.

Silencing negative regulators in p53 pathway enhances apoptosis without altering effective p53 integral

To further validate the concept of effective p53 integral, we simulated the effect of manipulating MDM2 expression. We first knocked down MDM2 by in silico lowering MDM2 expression. We found that MDM2 silence promoted steady accumulation of p53 and generally induced an earlier onset of terminal pulses (Fig. S3A). In addition, the average lifespan was consistently shortened by decreasing MDM2 (Fig. S3B). However, the effective p53 integral did not vary significantly in apoptotic cells (Fig. S3C). We also simulated the effect of Nutlin-3 treatment at 0, 6 and 16 h with Dox = 0.1. Simulating Nutlin-3 treatment in silico accelerated the accumulation of p53 (Figure 3(a)). The earlier addition of Nutlin-3 induced significantly shortened lifespan in cells (Figure 3(b)). However, the terminal pulses held constant across different regimes (Figure 3(c)). The total accumulative p53 (i.e. total p53 integral) steadily increased with different treatments (Figure 3(d)), the effective p53 integral was conserved irrespective of the time for Nutlin-3 addition (Figure 3(e)). We also silenced another negative regulator Wip1 in p53 signaling and found that lowering Wip1 in silico profoundly increased p53 accumulation and earlier onset of terminal pulses (FigureS4A). The lifespan was significantly lowered with reduced Wip1 expression (Fig. S4B). Not surprisingly, the effective p53 integral was not changed significantly with or without in silico si-Wip1 (FigureS4C). These results supported that repression of negative regulators in p53 pathway could promote apoptosis without changing the effective p53 integral. Altogether, our mathematical model faithfully recapitulated key experimental findings and could be used to further unravel the dynamic implications of dual-phase p53 patterns.

Figure 3. Simulation of effect Nutlin-3 treatment. (a) Representative traces of p53 treated with Dox = 0.1 combined with Nutlin-3 added immediately or at 6 and 16 h, following Dox treatment. To simulated the effect of Nutlin-3 treatment, k₁₁ and k₁₄ were simultaneously lowed by 5 fold. (b-e) Lifespan (b), terminal pulse height (c), total p53 integral (d) and effective p53 integral in apoptotic cells (e) for 500 simulations with nutlin-3 addministration immediately or at 6 and 16 h (Dox = 0.1), respectively. ANOVA was used to detect the statistical significance. **: P < 0.01.

Gene expression specificity determined by effective p53 integral

Wu et al. proposed that effective p53 integral cooperates with binding affinity to regulate downstream gene expression [13]. We then modeled this scenario by sampling two parameters (K_m and δ) associated with gene expression (see notes in Table S1 and Table S2). Then, stochastic simulation was performed for each parametric combination for three Dox doses (Dox = 0.00005 as control, Dox = 0.01 and 0.1 as stressed condition, Figure 4(a)). The integrated relative fold change of downstream effector with Dox = 0.1 compared with that with Dox = 0.01 was shown (Figure 4(a), colormap was normalized to the maximum expression under control condition). We found that higher K_m (i.e. lower binding affinity) and lower degradation rate could facilitate the fold change in gene expression between Dox = 0.1 and Dox = 0.01 (Figure 4(a)). The temporal dynamics of downstream effector induction further showed that slowly degraded effector with moderate or high K_m could discriminate between relatively high and low Dox (Figure 4(b)). Moreover, downstream effector with higher binding affinity (low K_m) failed to differentiate between Dox = 0.1 and Dox = 0.01 especially at early time points (Figure 4(b)). The scatter plot further supported that genes with lower binding affinity and degradation could enlarge the differences in gene expression between Dox = 0.1 and Dox = 0.01 (n = 50 for each dose, Figure 4(c) and Fig. S5). The effective p53 integral was also significantly larger in cells under 0.1 μM Dox treatment (Figure 4(d)). These results were qualitatively consistent with Wu et al.’s data [13] and further argued that gene expression specificity might be conferred via binding affinity of transcription factor and effector degradation rates.

Figure 4. Gene expression specificity by p53 binding affinity and target gene stability. (a) Relative fold change of target gene expression under 0.1 and 0.01 Dox treatment with varying binding affinity (K_m in μM, equivalent to k₁₁ in Table S1) and target gene degradation rates (δ, min⁻¹). Arrows represent exemplified situations in (B). (b) Stochastic simulations of target gene expression with varying transcription factor binding affinity and target gene stability. The mean trajectories (solid) for control (Dox = 0.00005, red), 0.01 (blue) and 0.1 (black) Dox were shown. Shaded areas represent 95% confidence intervals. (c) The integral of putative downstream effector (i.e. p53 target gene) under control, 0.01 and 0.1 Dox treatment was depicted with δ = 0.0001 (min⁻¹). N = 50 cases of stochastic simulations in (B) were shown. K_m values (from left to right: 0.01, 0.07, 0.15 and 0.3 μM) were displayed on top. (d) Effective p53 integral under 0.01 and 0.1 Dox treatment, respectively. Mann-Whitney test was used.

P53 overexpression facilitates terminal pulses and promotes apoptosis

We have simulated the threshold effect in p53 activation, we next argued that overexpressing p53 may also increase the effective p53 integral and favor apoptosis induction. By gradually increasing the production rate of p53 mRNA, we found that the system could ease the occurrence of p53 terminal pulses even at Dox = 0.1 (Fig. S6A). The fraction of cells with apoptosis was also increased with enhanced p53 overexpression (Fig. S6B). These data supported that manipulating p53 levels per se can also favor p53 terminal pulses and promote apoptosis.

Dynamic features of p53 signaling determines phenotypic responses

In stochastic simulation, we have found pronounced cell-to-cell viability with varying Dox levels (Fig. S7A). For example, dynamic responses when Dox = 1 either showed transient responses with sequential number of pulses or terminal pulses resembling those typically observed upon stimulation with Dox = 0.1 (Figure 5(a) and S7A). Therefore, dynamic responses were overlapping markedly between input levels and only partially determined by stimulus strength. To identify the encoded information despite the strong heterogeneity, we hypothesized that phenotypic responses were determined by p53 activation patterns rather than the amount of Dox applied to the system. We used k-mean clustering to group all signaling p53 dynamics into five classes (Figure 5(b)). The median time courses of five signaling clusters showed distinct dynamic patterns compared with corresponding cells defined by Dox levels (Figure 5(a,b) and S7. Note that to allow clustering, all temporal courses of p53 within 60 h were retained irrespective of apoptotic status). Cluster 1 showed on average a pulsatile behavior (low responder) while the other clusters (clusters 2–5) revealed significantly elevated p53 responses with varying duration (high responder, Figure 5(b)). However, grouping dynamics by Dox levels displayed a gradual evolution from pulsatile to high-amplitude elevation of p53 (Figure 5(a)). As expected, increasing Dox concentrations may induce a shift from low responder to high responders (Figure 5(c)). The shift is gradual implying that cells in a given dose could contain cells in several signaling clusters. Meanwhile, cells with the same Dox stimulation were more different in their dynamic behavior than cells grouped into a signaling cluster (Figure 5(d,e)). There were more positive silhouette scores in signaling classification by dynamic clusters than those in classification by Dox levels (Figure 5(d,e)). The phenotypic responses, such as lifespan and effective p53 integral showed in general a gradual dependence on Dox levels (Figure 5(f,g), top panels). However, these phenotypic responses exhibited a clear-cut pattern if p53 dynamics were clustered by signaling dynamics (Figure 5(f,g), bottom panels). Therefore, signaling clusters may provide a better separation of phenotypic responses compared with Dox doses. These results suggested that signaling dynamics rather than the stimulus doses encode for cellular behavior and could determine phenotypic outcomes.

Figure 5. P53 dynamics decompose into distinct signaling classes. (a) Time-resolved analysis of p53 dynamics for varying stimulus levels as indicated (Dox = 0.01–10). 20 representative cases were shown as well as the population median (red). (b) P53 dynamics were clustered into 5 signaling classes according to their time-resolved data. The population median values were shown (red). (c) Distributions of signaling classes depending on Dox levels. (D-E) Silhouette plots of cells sorted according to signaling classes (d) and Dox levels (e). Positive silhouette scores indicate that P53 dynamics are more similar to the own group, whereas negative scores indicate that the trajectory is closer to any of the other groups. (f) The lifespan distribution of cells treated with different Dox levels (top) or dynamic clusters (top). Note that if cell survived during the whole simulation period (60 h), the lifespan was recorded as 60 h. (g) Effective p53 integral for cells with varying Dox (top) and different clusters (bottom).

Decomposition into dynamic clusters confers sharp transition and less variability to effective p53 integral

Since the effective p53 integral can determine cell fate [13], we then investigated how variability of effective p53 integral evolves with distinct Dox doses or dynamic clusters. The distribution under low dose Dox treatment was skewed while it became more concentrated with increasing Dox doses (Figure 6(a)). However, there was a sharp transition of effective p53 integral from high to low variability between cluster 1 and 2–5 (Figure 6(b,c)). Notably, the variation gradually decayed with increasing Dox (Figure 6(c)). Instead, dynamic clusters clearly differentiate low (cluster 2–5) from high (cluster 1) variation in effective p53 integral (Figure 6(c)). Meanwhile, the distribution and variability were comparable among clusters 2–5 (Figure 6(b,c)). Next, we used information theoretical metrics to calculate the mutual dependence between effective p53 integral and dynamic processes (the rates v₁-v₁₈ were defined in Table S3). In most cases, the mutual dependence gradually decreased with increasing Dox (Figure 6(d), violet) whereas in principle a swift decline occurred from cluster 1 to clusters 2–5 (Figure 6(d), rose carmine). This was primarily due to the reduced variability (i.e. enhanced robustness) in effective p53 integral at higher Dox or cluster 2–5 (Figure 6(c)). Notably, the reduced mutual information for basal p53 (v₇) and Wip1 turnover (v₄ and v₁₅) with Dox = 0.01 or at cluster 1 indicated that Wip1 and p53 were kept at low levels with relatively less variation (Figure 6(d)). These data suggested that decomposition into dynamic clusters induced sharp transition of effective p53 integral.

Figure 6. Classification into dynamic clusters produces switch-like transition. (a) Distribution of effective p53 integral for varying Dox levels (0.01–10). (b) Distribution of effective p53 integral for five dynamic clusters. (c) The coefficient of variation for cells stimulated with various Dox (violet) or in distinct clusters (rose carmine). (d) Mutual information index with increasing Dox (violet) and different clusters (rose carmine). The details for kinetic processes (v₁-v₁₈) were listed in Table S3.

Discussion

Efficient information processing by the p53 pathway is crucial during stress response. Our current model qualitatively replicated the dual-phase dynamics of p53 in response to chemotherapeutic drugs. We modeled the observation that the occurrence of p53 terminal pulses correlated with cell death and manipulating p53 regulators could affect the terminal pulses without significant changes in effective p53 integral of apoptotic cells. Furthermore, we argued that it was the dynamic patterns could determine effective p53 accumulation. Dynamic clusters may provide less variability and more faithful discrimination for effective p53 integral and cell fate.

The occurrence of terminal pulses relied heavily on positive feedback loops. In current study, we considered the effect of MDM2-p53 mRNA binding to enhance p53 mRNA translation [15]. Notably, the p53 system is pervaded with positive feedback loops [8–10]. Incorporating any positive feedback loop in our model may support the p53 terminal pulses although the detailed mechanisms and kinetic parameters may differ (data not shown). Previous studies have shown that p53 may be activated in either pulsatile or sustained mode in response to DNA damage depending on stimulus type [4,11,12,16,17]. An intriguing hypothesis would be that multiple feedback loops in p53 network may display cell-specific activities and contribute to the diversity of p53 dynamics similar to the observation in MAPK signaling [18,19]. Furthermore, Purvis et al. found that sustained p53 activation induces senescence whereas pulsed p53 allows recovery and growth although pulsed and sustained p53 signaling show equivalent cumulative p53 levels [4]. Since the effective p53 integral under high doxorubicin treatment or sustained p53 activation might be significantly higher than that under pulsatile condition, the conceptual framework of effective p53 integral can also explain the observation in Purvis et al.’s work [4].

An earlier model developed by Zhang et al. also described a ‘two-phase’ p53 dynamics in response to DNA damage [20]. In their model, p53 only performs sequential pulses if low levels of DNA damage are applied. Moreover, when DNA damage is high, p53 always shows fixed number of pulses before the onset of sustained p53 activation [20]. However, a mild challenge by doxorubicin triggers sequential p53 pulses and occasionally results in terminal pulses. Applying a high dose doxorubicin treatment will abruptly lead to sustained p53 activation without the appearance of sequential pulses in Wu et al.’s experiments [13]. Therefore, Zhang et al.’s model at best is only valid under low dose doxorubicin treatment and cannot be used to characterize the sustained elevation of p53 in response to high doxorubicin doses.

The p53 dynamics can work cooperatively with binding affinity to differentiate the transactivation of target genes in Wu et al.’s work [13]. It is the effective p53 integral above the threshold rather than the absolute level of p53 integral that orchestrates gene expression and cell fate. Under mild doxorubicin challenge, although p53 pulses consecutively exceed the threshold, the cumulative integral above the threshold is insufficient to ignite apoptosis. When the DNA damage persists at high level, the accumulation of effective p53 then becomes competent to markedly transactivate proapoptotic genes and trigger apoptosis. The competence could be ascribed to the disruption of the balance between positive and negative feedback loops in p53 signaling to induce high-amplitude terminal pulses. Notably, we only focus on the dynamic aspects of p53-mediated cell fate decision. Other mechanisms such as post-translational modification, molecular interaction or subcellular translocation may provide additional complexity.

Our model suggested that lower binding affinity and degradation could better differentiate gene expression between high and low DNA damage (Figure 4). We noted that some key downstream factors such as MDM2 and PPM1D can mirror p53 dynamics indicative of short half-life [4,21]. However, factors involved in terminal fate (e.g. Bax, APAF-1 and PUMA) are relatively stable and show steady accumulation when p53 network becomes highly activated [22,23]. Therefore, the increased stability for p53 downstream regulators involved in apoptosis (at least for those mentioned in Wu et al.’s work [13]) could help differentiate the transactivation of different target genes and promote cell fate decision. Furthermore, previous study indicates that common p53 ‘binding affinity’ (K_m) for downstream targets is primarily distributed within 0.04 ~ 0.58 μM similar to the settings in our model (Figure 4) whereas apoptotic genes generally showed lower binding affinity [24]. Therefore, individual or combinatorial effect of stability and binding affinity may facilitate gene expression specificity and promote cell fate decision.

The signaling clusters represent mathematically identifiable classes for time series and provide more homogeneous grouping compared to stimulus concentrations. The definition of five clusters remains a heuristic choice [25]. Notably, p53 dynamics vary abruptly from cluster 1 to cluster 2–5 and provide clear-cut response patterns (e.g. lifespan, effective p53 integral), whereas grouping by doxorubicin doses represents a continuum of dynamic patterns. Noticeably, p53 temporal profiles will terminate when cells undergo apoptosis and therefore, most temporal p53 trajectories in cluster 2–5 will not persist for the entire 60 h (i.e. profiles in cluster 2–5 become more homogeneous). As a result, grouping by temporal profiles may better discriminate cell fate. In cluster 1, we found that on average cells performed consecutive pulses with occasionally terminal pulses. Profiles in cluster 1 indicated low effective p53 integral and longer lifespan, implying that these dynamic patterns favored survival (survival cluster). Instead, nearly all profiles in cluster 2–5 ended in p53 terminal pulses with significantly reduced lifespan and high effective p53 integral (dead cluster). Therefore, the dynamic clustering suggested that pulsed p53 may favor survival until the occurrence of high-amplitude terminal pulses. The stimulus strength will determine the fraction of cells in each cluster and the fraction in dead cluster was gradually elevated with increasing stimulus. Collectively, manipulating wild type p53 levels (i.e. effective p53 integral) might be an effective strategy to tip the balance between survival and death.

There are also limitations in our current work. We did not consider the subcellular localization of specific proteins although this may possibly confer enhanced ultrasensitivity/nonlinearity in p53 signaling. ATM-independent effect on DSB repair was also not considered as it only affected the repairing rate but not the feedback architecture in p53 signaling [26]. Previous studies showed that MDMX and c-Abl also affect p53 levels [27,28]. However, these regulatory mechanisms do not function through feedback mechanisms and therefore do not exert significant influence on dual-phase nature of p53 under doxorubicin treatment. We only emphasized the effect of interlocked negative and positive feedback loops in p53 network and other detailed regulatory effects were not considered for simplicity.

In conclusion, we have modeled the dual-phase p53 dynamics under doxorubicin treatment and elucidated the putative determinant for p53-mediated cell fate decision. As aberrant p53 signaling is usually linked to tumor progression, our current model may help ignite new clues to optimized clinical design for chemotherapeutic drugs.

Materials and methods

Model construction

We modeled the p53 regulatory network containing 11 species and 30 reactions (Figure 1). For p53, MDM2 and Wip1, we incorporated both mRNA and protein representation. The conventional genotoxic drug doxorubicin (Dox) was shown to induce DNA double strand breaks (DSBs) with a half-life around 11.7 h [29] corresponding to an approximate clearance rate of 0.001 min⁻¹. Therefore, we assumed that Dox clearance underwent an exponential decay with a clearance rate of 0.001 min⁻¹ in stochastic simulation. The DNA damage can phosphorylate and activate the checkpoint kinase Ataxia-telangiectasia mutated (ATM) [30]. Total ATM levels were considered as a constant because data have shown that total ATM does not change significantly within 72 hours post stimulation [31]. Activated ATM (ATM*) can phosphorylate p53 and MDM2, which become stabilized or destabilized, respectively [32,33]. We incorporated two negative feedback loops (P53-MDM2 and P53-Wip1-ATM) as suggested by Lahav et al [34]. In response to DNA damage, p53 can be phosphorylated, activated and induce MDM2 expression, which in turn targets p53 for proteasome degradation. The second negative feedback loop involves p53 induced Wip1 upregulation and ATM deactivation [35]. The levels of Dox denote model input and affect the DNA damage induced ATM activation as previously [36] (Table S1 and S2).

Gajjar et al. have found that under doxorubicin treatment, ATM can phosphorylate MDM2 at Ser395 and allosterically altered MDM2 can bind p53 mRNA [15]. The phosphorylated MDM2 at Ser 395 can further enhance p53 mRNA translation [15]. We incorporated this interaction in our model and assumed that p53 mRNA can bind phosphorylated MDM2 (MDM2p) and form a p53 mRNA-MDM2p complex. We also assume that this complex can promote p53 protein translation with higher efficiency (Table S1 and S2). We did not consider the subcellular localization not only for simplicity but also owing to the fact that the detailed mechanism by which MDM2-dependent nucleolar targeting of p53 mRNA enhances its translation and its detailed subcellular trafficking events remain elusive [15]. All the dynamic reactions were formulated using ordinary differential equations (ODEs) (Table S1) with parameter description listed in Table S2.

Model parameterization

p53 mRNA has a long half-life from 10 to 22–23hours in different cells [37]. We assigned 0.001 min⁻¹ (k₂) corresponding to a half-life ~11.6 h. The half-life of MDM2 mRNA is reported to be within 1 to 2 hours [38,39] implying a degradation rate from 0.005775 to 0.01155 min⁻¹. Therefore, we set this parameter to be 0.006 (k₄). Basal mdm2 induction rate (k₃), MDM2 dependent P53p degradation rate (k₁₄), mdm2 translation rate (k₁₅), basal MDM2 degradation rate (k₁₆), the threshold concentration of MDM2 dependent p53 (K₃) and P53p degradation (K₆) were set according to Ma et al’s model [40]. The turnover rate of p53 (k₁₀) and MDM2 dependent p53 degradation (k₁₁) were assigned according to Cai et al [37]. MDM2 undergoes phosphorylation and accelerated degradation in DNA damage response and we therefore assigned a roughly 10-fold increase in MDM2p degradation rate (k₁₉) [41]. Wip1 dependent dephosphorylation rate (k₁₃) of phosphorylated p53 (P53p), basal Wip1 degradation rate (k₂₁) and Wip1 mediated ATM deactivation rates (k₂₄ and K₁₀) were derived from Zhang et al within the same order of magnitude [20]. The dissociation constant k_r/k_f was set to be 100 nM qualitatively similar to Gajjar et al.’s measurement [15]. The basal p53 production rate (k₁), ATM* induced MDM2 activation (k₁₇), the threshold concentration of p53 induced mdm2 (K₁) and wip1 (K₂) transcription were determined according to Chong et al’s data [42]. All other parameters were estimated from our previous work [43] to reproduce the dual-phase dynamics of p53.

Local parameter sensitivity

Local parameter sensitivity analysis provides dynamic responses to an infinitesimal disturbance in kinetic parameters. A dynamic system can be formulated by x’ = F(x, p), where x and p donate state vector and parameter vector for the system, respectively. Pulse amplitude and width sensitivity capture the variations of relative amplitude and duration in response to parametric changes. Relative amplitude sensitivity S_A and width sensitivity S_W are defined as

(1) $S_A=\frac{\mathrm{\partial }A/A}{\mathrm{\partial }p/p}=\frac{\mathrm{\partial }ln(A)}{\mathrm{\partial }ln(p)}$(1)

(2) $S_W=\frac{\mathrm{\partial }W/W}{\mathrm{\partial }p/p}=\frac{\mathrm{\partial }ln(W)}{\mathrm{\partial }ln(p)}$(2)

where A and W denote relative pulse amplitude and width, respectively. Note that normalized sensitivities are only locally valid in parameter space.

Mutual information index (MII)

A general measure of input-output associations can be provided by mutual information regardless of the shape of underlying distributions. Randomized perturbations of input parameters (X₁, X₂,…, X_N) in parameter space within predefined range create a random output Y with a probability density p(y) . The entropy associated with Y is defined as

$H(Y)=-\sum _{y}p(y)\cdot {log}_{2}p(y)$

Conditional entropy can further be defined with one fixed X_i

$H(Y|X_i)=\sum _{x}p(x)\cdot H(Y|X_i=x)$

The mutual information is defined as the difference between total entropy H(Y) and one preferred conditional entropy H(Y | X_i) to characterize the impact X_i exerts on Y.

$I(X_i;Y)=H(Y)-H(Y|X_i)$

Mutual information index (MII) can be defined as [44]:

$m_i=\frac{I(X_i;Y)}{H(Y)}$

Metrics for terminal pulse, lifespan and effective p53 integral

Prolonged low-dose or acute high-dose doxorubicin treatment initiated a series of p53 pulses followed by an abrupt increase to high p53 amplitude (defined as ‘terminal pulse’) [13]. The lifespan denotes the time duration from initial Dox treatment to the time point when apoptotic effector (e.g. Bax) reaches a pre-defined threshold (see below). According to Wu et al’s reports, p53 accumulation below a certain threshold may not activate proapoptotic genes [13]. This threshold is defined as ‘minimal effective level’ (MEL) and the p53 accumulations (i.e. temporal integral) above the threshold are termed as ‘effective cumulative p53ʹ.

$Ep53=\underset{0}{\overset{T}{\int }}(P53-MEL)dt$

T denotes the time when apoptosis occurs. Note that in our simulation, the ‘MEL’ was assigned 1.8 × 10⁴). The time-averaged Ep53 is defined as follows:

$Ep{53}_{AverT}=\frac{1}{T}\underset{0}{\overset{T}{\int }}(P53-MEL)dt$

In current study, ‘effective p53 integral’ refers to ‘effective cumulative p53ʹ and is sometimes abbreviated as Ep53.

Stochastic simulation

Since the maximum level of dynamic species was >10⁵, we performed stochastic simulations by binomial τ-leap method according to Chatterjee et al and Leier et al’s work [45,46]. The kinetic delays were varied by 10% around the reference values during the simulation. We further assumed that the transcription is burst-like and the burst size positively correlated with the amount of transcription factors [47]. In stochastic simulation, we assumed that Dox obeys an exponential decay with clearance rate 0.001. Meanwhile, p53 can transactivate the expression of an apoptotic effector (e.g. Bax in our model). If the Bax level reaches the threshold (2.1 × 10⁴), cells will immediately trigger apoptosis owing to the ‘Bax activation switch’ described in previous work [48]. Therefore, the p53 trajectories may terminate at different time.

Model simulation

The ordinary differential equations were numerically integrated using the ode23s solver. Both stochastic and deterministic simulations were performed using MATLAB (MathWork, Version 8.1, R2013a).

Acknowledgments

This work is supported by National Natural Science Foundation of China (31400714, 31800316), Key Project of Natural Science Foundation of Anhui Provincial Education Department (KJ2017A359).

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the Young Scientists Fund [31800316]; Young Scientists Fund [31400714].