Understanding the mechanism of cellular response to ionizing radiation (IR) damage is important from the perspectives of both radiotherapy and mitigation of radiation damage. Cell response to IR involves several protein-DNA and protein-protein interactions, as well as the formation of free radicals that alter cellular biochemistry[@b1]. Cell death usually takes place several hours after radiation injury. Even if the exposure to radiation is brief, its effect on cellular biochemistry may be long-lived depending on the strength of IR[@b1]. Moreover, several proteins that are expressed transiently after radiation damage may trigger downstream responses that are manifested long after the original insult. The responses to treatments that aim at alleviating radiation damage (or decreasing the susceptibility to apoptosis in damaged cells) depend on the dosage and duration of exposure, the treatment timing, and the dynamics of the proteins that regulate apoptotic events.

The tumor suppressor protein p53 is a main mediator of cell response to genotoxic stress. p53 regulates apoptosis via both transcription-dependent and -independent pathways[@b2][@b3], in addition to regulating cell/tissue-specific response to radiation by apoptosis-independent mechanisms[@b4]. The transcription-independent effect of p53 is mediated by its translocation to the mitochondria, although the mechanism is still debated (see review[@b5]). Previous efforts to model cell response to radiation have been in part stimulated by the observed oscillatory dynamics, or repeated pulses, of p53 in response to radiation damage[@b6][@b7][@b8]. To this end, deterministic methods[@b6][@b7][@b8][@b9] and, to a significantly lower extent, stochastic simulations[@b10] have been adopted. Apoptosis itself has been mathematically modeled independent of p53 response to radiation, using deterministic[@b11][@b12][@b13] as well as probabilistic methods[@b14][@b15]. Likewise, there have been efforts to establish the link between p53 activities to DNA damage and cell fate using deterministic simulations[@b16][@b17] and methods of limited stochasticity[@b18].

With accumulating experimental data, we are now in a better position to construct more detailed models for p53-mediated signal transduction in response to IR and use them as a platform for evaluating new polypharmacological strategies. Here, we focus on the biochemical network associated with IR-induced apoptosis and examine the time-dependence of p53-mediated apoptotic events. Our approach incorporates cell heterogeneity and subcellular localization, and aims at estimating the response to targeted therapies following IR. It applies to conditions where certain types of molecules are very small in number yet are major determinants of system behavior.

We consider several outstanding issues: (i) the significance of the oscillatory behavior of p53 in determining the onset of downstream apoptotic events in response to IR, (ii) the role of its transcription-dependent and -independent activities in regulating cell susceptibility to apoptosis, (iii) the effect of the interactions involving anti-apoptotic Bcl-2 and pro-apoptotic Bax on cell fate, (iv) the impact of the positive feedback loop mediated by Bid/caspase-3, and (v) the efficacy of various treatment strategies, e.g. how particular combination therapies may elicit anti-apoptotic responses to mitigate IR-induced damage.

Our results indicate that p53 oscillations are insufficient to induce apoptosis *per se*. Activation of Bax on the outer mitochondrial membrane (OMM) plays a key role in driving apoptosis, and Bid and Bax emerge as targets most relevant for regulating apoptosis. In particular, truncated Bid (tBid)/caspase-3 feedback loop is a determinant of cell fate or treatment efficacy: if the feedback is sufficiently strong, inhibition of pro-apoptotic proteins like PUMA mitigates damage; but the effect weakens with delay in inhibitor administration. These results offer new insights into novel polypharmacological strategies for alleviating IR damage.

Results
=======

Mathematical modeling of IR-induced apoptosis
---------------------------------------------

We first constructed a mathematical model for the biochemical network associated with IR-induced apoptosis. Our model consists of three closely interconnected subnets ([Figure 1](#f1){ref-type="fig"}): an upstream p53 module, a mitochondrial module and caspase activation events triggered upon mitochondrial cytochrome *c* (cyt *c*) release. The interactions are distributed across three subcellular locations: nucleus (N), cytosol (C) and mitochondria (M), indicated by superscripts.

### p53 module

The p53 module consists of the translocation of cytoplasmic p53 (p53^(C)^) to the nucleus, especially when the cell is exposed to stress[@b19], p53^(N)^ tetramerization into (p53^(N)^)~4~[@b20], and the transcriptional activation of Mdm2 by p53. The latter involves four steps: (i) generation of messenger RNA, mRNA~Mdm2~^(N)^, (ii) its translocation to the cytoplasm, (iii) translation of mRNA~Mdm2~^(C)^ into Mdm2^(C)^, and (iv) translocation of Mdm2^(C)^ to the nucleus, which serves as a negative regulator of p53 by promoting its unbinding and ubiquitination[@b21]. p53~ub~^(N)^ can translocate to the cytoplasm, and p53~ub~^(C)^ to the mitochondria[@b22].

Under severe genotoxic stress, p53^(N)^ gains transcriptional activity (p53~A~) for pro-apoptotic BH3-only proteins[@b23], including PUMA and mitochondrial outer membrane permeability (MOMP) inducers represented by Bax[@b24]. In parallel, it retains its ability to tetramerize and activate mRNA~Mdm2~, and to translocate to the mitochondria.

The p53-Mdm2 interaction is mediated by upstream DNA-damage-signaling proteins such as ATM, Chk2 or Wip1[@b8]. DNA double-strand breaks resulting from exposure to IR induce ATM autophosphorylation[@b25], which then disrupt the p53-Mdm2 interaction upon phosphorylating both proteins. We represent ATM-associated upstream events as follows: DNA-damage activates ATM, which, in turn, stabilizes p53.

### p53-modulated mitochondrial events and release of cytochrome c

Mitochondrial events modulated by p53 proceed via two paths: transcription of PUMA and Bax, and translocation of p53 to the mitochondria followed by binding of p53^(M)^ to Bcl-2 to inhibit its anti-apoptotic action[@b26]. PUMA also binds Bcl-2 on the mitochondrial membrane[@b27] with a higher affinity than does p53^(M)^ and can displace p53^(M)^ from its complex with Bcl-2[@b3], thus, freeing up p53 for further activity[@b28].

Bax translocates between the cytoplasm and the mitochondria[@b29], its retranslocation into the cytosol being mediated by Bcl-x~L~[@b29]. Bax^(M)^ interacts with Bcl-2 to form a complex[@b30] that prevents its activation (into Bax\*) and ensuing oligomerization to form a MOMP pore. However, p53^(M)^ can bind to Bcl-2 stronger than does Bax and displace Bax, thus countering/alleviating this effect[@b28]. Activation of Bax^(M)^ is facilitated by the localization of tBid to the OMM[@b31], which induces the insertion of Bax into the OMM[@b32]. Bax is also activated by p53^(M)3^ and by PUMA[@b33]. Bax\* oligomerizes on the OMM[@b34], to form a MOMP pore, which, in turn, promotes the release of cyt *c*^(M)^ into the cytoplasm[@b11]. Cyt *c* release is usually considered as the point of no return in mitochondria-mediated apoptosis. MOMP pore also enables the release of Smac/Diablo^(M)^ that inactivates the inhibitors of apoptosis (XIAPs), further promoting apoptosis[@b35].

### Events triggered by cyt c release

Cyt *c*^(C)^ forms a complex with Apaf-1 in an ATP-dependent manner, which assembles into the apoptosome complex[@b36] upon heptamerization, and recruits inactive procaspase-9 molecules to activate them into caspase-9 (C9) and catalyze the cleavage of procaspase-3 (proC3) to form active C3[@b13]. C3 and C8 activated by external death signals[@b37] truncate Bid, resulting in a positive feedback loop that amplifies cyt *c* release[@b38], while XIAP inhibits the apoptosome[@b39] and promotes the proteasomal degradation of C3[@b40].

### Synthesis, degradation and inhibition of components

Synthesis and degradation of monomeric species (not shown in [Figure 1](#f1){ref-type="fig"}; see [Tables S1--S2](#s1){ref-type="supplementary-material"}) help establish and maintain steady state conditions in the absence of stimuli. We consider four drug targets: PUMA, Bid, C3 and Bax. Their inhibitors are designated as I~PUMA~, I~Bid~, I~C3~ and I~Bax~, respectively.

Model simulation, calibration and validation
--------------------------------------------

We adopted stochastic simulations for two reasons. First, the quantity of some proteins such as caspase-3 are expected to be extremely low (or non-existent) under homeostatic conditions. Second, as [Figure 2A](#f2){ref-type="fig"} shows, stochastic simulations reproduce the sustained oscillations of p53^(N)^ and Mdm2^(N)^ in accord with experiments[@b6][@b7][@b41], while deterministic simulations result in damped oscillations.

To establish initial concentrations of proteins that take into account of cell-to-cell variability, we first performed a run of 300 h to allow the system to reach steady state conditions in the absence of IR. Then, we changed the system parameters to account for IR-induced perturbations, and allowed the system to evolve. Return to unstressed state occurs upon restoring the parameters after Δ*t*. [Figure 2B](#f2){ref-type="fig"} shows the p53^(N)^ levels as a function of time for transient (Δ*t* = 12 h) and sustained (Δ*t* = 56 h) exposure to radiation. Note that Δ*t* *does not* refer to the duration of radiation but to the length of time the kinetic parameters are altered from unstressed (normal) to radiation (WT) values (see [Supplementary Table S3](#s1){ref-type="supplementary-material"}). The transient and sustained exposures are analogous to 0.3 Gy and 10 Gy of γ-irradiation, respectively. In what follows, simulations were performed under sustained exposure condition.

We obtained kinetic parameters from experiments whenever available and deduced others from the ranges adopted in previous computational studies[@b11][@b13]. Details on the sources and choices of model parameters are presented in the [Supplementary Text](#s1){ref-type="supplementary-material"}. These were further calibrated to accurately describe previous experimental observations[@b6][@b7][@b41].

[Figure 3A](#f3){ref-type="fig"} shows the comparison of the model-predicted time profiles of p53^(N)^ and Mdm2^(N)^ (*blue lines*) with time series western blot data[@b41] (*red dots*). 1,000 stochastic trajectories were generated (under the same initial conditions but with different random seeds) and averaged to mimic the heterogeneous cell population behavior. Good agreement is achieved between predicted time profiles and experimental data.

We next validated the model using single-cell-based experimental data[@b7][@b42]. We compute the histogram of the period and amplitude values of Mdm2^(N)^ oscillations using the 1,000 trajectories simulated under IR. As shown in [Figure 3B and C](#f3){ref-type="fig"}, the simulation results match the experimentally detected period[@b7] and amplitude[@b42] of oscillations, indicating that the model captures well cell-to-cell variability property.

Our model also reproduced key observations downstream the signaling cascade ([Figure 3D](#f3){ref-type="fig"}). Note that two stochastic simulations may lead to different dynamics even with identical initial conditions, and conclusions drawn from a single simulation may not be reliable. We applied a statistical model checking (SMC) technique[@b43] in order to make predictions with given confidence level. Briefly, we encode dynamical properties of the system using a bounded linear temporal logic and devised a SMC procedure to check if the system that satisfies a given property passes a statistical test with predefined type-I (false positive) and type-II (false negative) errors. A detailed description can be found in the Methods. [Figure 3D](#f3){ref-type="fig"} summarizes the dynamical properties predicted with high confidence. These properties are consistent with existing experimental observations[@b7][@b12][@b41][@b44][@b45][@b46][@b47]. This analysis further validates the model and supports its use for predictive purposes, presented next.

Sensitivity analysis
--------------------

In order to evaluate the sensitivity of the model to initial concentrations and kinetic parameters we performed a sensitivity analysis based on the level of caspase-3 activation. The results show that the system is robust to perturbations in the initial concentrations of major species, in support of the approach adopted for defining initial distributions.

The calculated sensitivities to model parameters are presented in [Supplementary Figure S2](#s1){ref-type="supplementary-material"}. Strong control over the system response is distributed among the p53-Mdm2 interactions (*k~45~*), activation of the nuclear p53 (*k~11~*) and ensuing transcriptional activation of PUMA (*k~20~*), translocation/retranslocation of Bax between the mitochondria and the cytosol (*k~24~*^±^), Bax activation (*k~30~*) and caspase-3 regulation (*k~38~^+^* and *k~42~*). The reactions with high influences are colored red in [Figure 1](#f1){ref-type="fig"}.

The analysis highlights the significant role of p53-mediated transcription-dependent pathway: . The p53-transcription-independent events (involving Mdm2 and p53^(M)^) have moderate effects and become influential when they feed into the Bax activation pathway. These results suggest that the strength of coupling of the nuclear p53 module to the mitochondrial apoptotic machinery will play an important role in determining the cell fate. In the next sections, we focus on the effects of the reactions identified here to be most influential. We begin with the examination of the relative contributions of the p53-transcription-dependent and p53-translocation-dependent activities of p53 to eliciting mitochondrial apoptotic responses.

Coupling of the p53 module to mitochondrial pore opening
--------------------------------------------------------

p53~A~ and p53^(M)^ initiate the p53-transcription-dependent and -independent apoptotic events, respectively, by either upregulating pro-apoptotic proteins (represented by PUMA and Bax), or interacting with antiapoptotic Bcl-2 family proteins. These may lead to MOMP opening, cyt *c* release and caspase activation depending on the coupling between the two pathways, which in turn depend on the relative rates of (i) transcriptional activation of pro-apoptotic proteins by p53^(N)^ (*k~11~*) and (ii) translocation of p53^(C)^ to the mitochondria (*k~10~*). We varied *k~11~* and *k~10~*, and simulated the time profiles of caspase-3 activation (commonly used as an indicator of apoptosis) toward elucidating the relative sensitivity of caspase-3 activation to these two p53-mediated mechanisms - higher *k~11~* and *k~10~* values representing the dominance of transcription-dependent and -independent activities of p53, respectively. The results are presented in [Figure 4A and 4B](#f4){ref-type="fig"}.

[Figure 4](#f4){ref-type="fig"} shows that the onset of apoptosis is elicited either by enhancing the transcriptional activity of p53 (increasing *k~11~*) (**A**), or accelerating the p53^(C)^ translocation to the mitochondria (increasing *k~10~*) (**B**), while low efficiency of either process is not sufficient to cause apoptosis. To make an assessment of the efficacy of either pathway on determining the cell fate, we further predicted the response curves of caspase-3 activation to *k~11~* and *k~10~*. As shown in [Figure 4C](#f4){ref-type="fig"}, the apoptotic response to *k~11~* is clearly greater than that to *k~10~* indicating that the p53-transcription-dependent pathway is playing a dominant role in determining cell fate.

Next, we examined the role of p53 oscillation in mediating apoptotic response. The p53 level oscillate with a \~5.5 h period due to radiation exposure ([Figure 2B](#f2){ref-type="fig"}). It is unclear if the apoptotic response is mainly determined, or predominantly affected, by the p53 oscillation frequency. To address this question, we perturbed the parameters ([Supplementary Table S3](#s1){ref-type="supplementary-material"}) to produce a condition, under which the p53 level oscillates with a lower frequency (LF) (\~11 h period) in response to IR ([Supplementary Figure S1](#s1){ref-type="supplementary-material"}). Under such conditions, the apoptotic response previously observed for *k~11~* = 2 × 10^−5^ s^−1^ remained unaffected ([Figure 4D](#f4){ref-type="fig"}), suggesting that slowing down p53 oscillations does not dampen the p53-transcription-dependent apoptotic response. Interestingly, an apoptotic response could be triggered even with a moderate transcriptional activity of p53 (*k~11~* = 10^−5^ s^−1^, [Figure 4D](#f4){ref-type="fig"}), which would fall short of inducing apoptosis under WT conditions (see panel **A**). This might be explained by the increase of p53^(C)^ level ([Supplementary Figure S3](#s1){ref-type="supplementary-material"}) due to LF condition, which, in turn enhances the p53-transcription-independent pathway. These results imply a dependency of the p53-mediated apoptosis on the p53 oscillation frequency, and also suggest that p53-transcription-independent pathway might be capable of defining cell fate when p53-transcription-dependent pathway is suppressed.

A crucial insight that emerges from these findings is that IR-induced cell death depends on the coupling between the regulatory p53 module and downstream apoptotic cascade of events, and among the two regulatory activities of the p53 module, transcription of pro-apoptotic proteins appears to be more effective that those driven by the translocation of p53 to the mitochondria.

Next, we examine the response of proteins downstream of the p53 module toward shedding further light on the significance of particular interactions/reactions in determining cell fate.

Bax activation
--------------

The above analysis highlights the importance of transcriptional regulation by p53, while also drawing attention to the complementary role of transcription-independent mechanisms. Bax is actually affected by both mechanisms. The sensitivity analysis above indeed highlighted the importance of Bax activation (*k~30~, k~51~, k~74~*) for initiating apoptosis ([Supplementary Figure S2](#s1){ref-type="supplementary-material"}). As shown in [Figure 1](#f1){ref-type="fig"}, in the p53-transcription-dependent pathway, p53~A~ prompts Bax activation by upregulating the expressions of Bax and PUMA which activates Bax either directly or indirectly (via binding Bcl-2 and thereby preventing its anti-apoptotic effect on Bax). Similarly, p53^(M)^ in the p53-transcription-independent pathway also activates Bax either directly or indirectly via binding Bcl-2. Thus, the affinity of Bax^(M)^ for Bcl-2, accounted for by the ratio *k~25~^−^*/*k~25~^+^*, controls the efficacy of the indirect Bax activation (via either p53-transcription-dependent or independent pathways).

*k~25~^−^*/*k~25~^+^* widely varies with experimental conditions and cell type[@b26][@b27]. Values spanning several orders of magnitude, from micromolar to nanomolar, have been reported for this dissociation constant by different groups (see the [Supplementary Text](#s1){ref-type="supplementary-material"}). [Figure 5](#f5){ref-type="fig"} shows caspase-3 and Bax^(M)^.Bcl-2 complex concentrations for two different *k~25~^−^*/*k~25~^+^* values. Comparison of [Figure 5A and C](#f5){ref-type="fig"} shows that the caspase-3 levels (*blue* curve, *right* ordinate) are unaffected by variation in the dissociation constants over two orders of magnitude. In contrast, vastly differing amounts of Bax^(M)^. Bcl-2 complex may be formed (*black* curve, *left* ordinate) consistent with the change in the dissociation constant. Similarly, the different dissociation constants had an insignificant effect under anti-apoptotic conditions (panels 5B and D).

The model is therefore robust to variations in Bax^(M)^.Bcl-2 association strength. This suggests that cells exhibiting significantly different propensities for Bax^(M)^.Bcl-2 complex formation or dissociation can still exhibit the same cellular fate. It also implies that the indirect activation of Bax (via its dissociation from the complex Bax^(M)^.Bcl-2) has only a secondary effect on the overall cell response. This conclusion is further supported by results presented in [Supplementary Figure S4](#s1){ref-type="supplementary-material"}, which shows that knocking down the interactions between PUMA and Bcl-2 does not significantly decrease caspase-3 levels ([Supplementary Figure S4](#s1){ref-type="supplementary-material"}).

Caspase/tBid positive feedback loop
-----------------------------------

Bax activation leads to the release of cyt *c* and causes caspase-3 activation. Caspase-3 cleaves Bid and the resulting tBid^(M)^ further activates Bax. In this section, we examine the role of this positive feedback loop, Bax\* → C3 → tBid → Bax\*.

[Figure 6](#f6){ref-type="fig"} displays the time evolution of PUMA, tBid^(M)^ and caspase-3, after exposure to IR. PUMA (*black*) mirrors the p53 oscillations, presumably due to its direct transcriptional regulation by p53. Interestingly, PUMA levels decay after termination of IR exposure, but there is sustained tBid fluctuations and newly elicited caspase-3 production even after the decay of PUMA. Sustained activation becomes more prominent with increasing coupling to p53 transcriptional machinery (high *k~11~*, [Figure 6B](#f6){ref-type="fig"}). Caspase-3 activation depends on whether sufficient tBid is activated by the time PUMA decays. The increase in \[tBid\] (and associated positive feedback to Bax) is critically important for sustained caspase activity. Otherwise, PUMA and Bax upregulation may fall short of triggering efficient cyt *c* release. Under these circumstances, the extent of Bax activation by p53^(M)^ (*k~10~, k~29~*) or PUMA (*k~30~*) may become critical.

The observed importance of the caspase/tBid feedback loop in modulating apoptosis warrants further study. Two of the possible ways that this strength can differ in cells is through truncation of Bid by caspase-3 and through activation of Bax by tBid (controlled by *k~28~*). We varied *k~28~* to study the effect of this feedback on cell fate (high *k~11~*, [Figure 6C-D](#f6){ref-type="fig"}). Unless *k~28~* is sufficiently high (*k~28~* = 0.1 μM^−1^s^−1^), no activation of Bax^(M)^ by tBid takes place, i.e., the positive feedback loop that promotes apoptosis is decoupled from the apoptotic machinery, and Bid has no influence on apoptosis. In the intermediate values of *k~28~* = 0.3 μM^−1^s^−1^, caspase-3 reaches a level comparable to that attained with *k~28~* = 0.5 μM^−1^s^−1^, although the response is slower. Therefore, in contrast to the strength of Bax/Bcl-2 association, the intensity of the positive feedback loop significantly affects the time evolution of \[C3\], suggesting that targeting the caspase-3/tBid loop might be an effective therapeutic strategy for mitigating radiation damage.

Efficacies of polypharmacological strategies
--------------------------------------------

In the sections above, we have analyzed the dynamics of apoptotic mediators in response to IR. We now focus on identifying potential targets and polypharmacological strategies suitable for mitigating radiation damage. Of special interest are treatment strategies that are effective even if not administered immediately after IR exposure.

[Figure 7](#f7){ref-type="fig"} illustrates the efficacy of simulated drug treatments administered at two different times after exposure to IR: immediately (with a time lag of 15 minutes succeeding IR exposure, [Figure 7A](#f7){ref-type="fig"}) and after a delay of 12 hours ([Figure 7B](#f7){ref-type="fig"}). The treatments involve administration of individual inhibitors for four targets, PUMA, Bid, caspase-3 and Bax, or combined therapies targeting pairs of these proteins. All (virtual) inhibitors are assumed to have nanomolar binding affinity to their target proteins. In the absence of treatment, the given radiation dose leads to sustained caspase-3 activity.

Individual inhibitions of PUMA, Bid or Bax are very effective in mitigating radiation damage, provided that treatment is available immediately after radiation damage. On the other hand, if there is a longer delay before the treatment is administered ([Figure 7B](#f7){ref-type="fig"}), individual inhibitions of Bax and Bid still help in abrogating radiation damage but PUMA inhibition is no longer an effective treatment. The ineffectiveness of PUMA inhibition is due to the fact that a 12-hour PUMA activity is sufficient for the activation of the tBid/Caspase-3 positive feedback loop and activation of Bax via tBid dominates the apoptotic response at this late stage.

In contrast to inhibiting PUMA, Bid and Bax, inhibition of caspase-3 itself does not prevent apoptosis irrespective of the timing of inhibition. This is an intriguing observation since complete abrogation of caspase-3 would be expected to inhibit apoptosis (although not caspase-independent cell death). However, we have a continual production of procaspase-3, and the activation of caspase-3 upon cleavage of procaspase-3, which is enabled by Bid and PUMA, dominates the outcome. This observation suggests that targeting caspase-3, itself, may not be an efficient strategy for preventing apoptotic response.

All the combination therapies considered here abrogate apoptosis completely (e.g., inhibition of PUMA + Bax, Bid + Bax, PUMA + Bid), or significantly (Bid + C3) when administered immediately after radiation exposure. This might be expected as I~PUMA~, I~Bid~ and I~Bax~ are already effective alone and their joint administration apparently does not have an unexpected effect. When the inhibitors are administered after 12 hours, the combination I~PUMA~ + I~C3~ is ineffective, consistent with the failure of these two inhibitors to prevent caspase-3 accumulation individually ([Figure 7B](#f7){ref-type="fig"}). However, the time dependence of \[C3\] is interesting: it shows a short-term response suggestive of mitigation of apoptosis, followed by an apoptotic response after a day or two. This combination therapy is therefore not effective -- again, due to the attainment of sufficiently high level of tBid to sustain the pro-apoptotic positive feedback loop that is enhancing Bax activation. We then simulated the attenuation of the strength of this feedback loop by reducing *k~28~* (from 0.5 to 0.3 μM^−1^s^−1^). The results confirmed that I~PUMA~ administered after 12 hours would be effective in mitigating radiation-induced apoptosis provided that this positive feedback loop is weakened (*blue* curve, [Figure 7D](#f7){ref-type="fig"}).

These results highlight the importance of examining the long-term behavior before making an assessment on the efficacy of a treatment, the role of timing in ensuring an effective treatment, and the significance of identifying most susceptible targets for efficacy of the treatment.

Discussion
==========

Recently, Batchelor *et al* reported that stimulus-dependent temporal dynamics of p53 is an important determinant of cell regulation[@b48]. A further study suggests different p53 dynamics (repeated pulses vs sustained response) may alter cell fate[@b47]. The present study focused on a detailed examination of the dynamics of p53-mediated and mitochondrial-dependent modules for regulating apoptosis in response to IR damage. It also highlights the significance of the timing between genotoxic stress and therapeutic intervention in the context of the cell stochastics.

Our model analysis showed that the p53 oscillatory behavior in response to IR is not, *per se*, sufficient to explain a cell\'s susceptibility to apoptosis; rather, the dynamic coupling between the p53 module and mitochondrial machinery is important. Among transcription-dependent and -independent roles of p53, the former plays a major role in driving apoptosis via transcriptional activation of Bax. The present study also sheds light to the strength of caspase/tBid feedback loop as a determinant of apoptotic response as well as treatment efficacy.

Several experimentally testable results emerge from the analysis: (i) oscillatory behavior of PUMA in response to radiation damage ([Figure 6A](#f6){ref-type="fig"}), (ii) elevated tBid activity even after the decay of PUMA when increased caspase-3 levels occur ([Figure 6B](#f6){ref-type="fig"}), (iii) inhibition of Bid and/or Bax as more effective strategies for mitigating radiation damage, compared to inhibition of PUMA ([Figure 7](#f7){ref-type="fig"}), and (iv) the need to have more potent inhibitors of caspase-3, compared to those of PUMA, Bid or Bax, in order to effectively mitigate radiation-induced cell death.

The last two points are especially relevant to developing efficacious therapies to regulate apoptosis. Our earlier experiments and pharmacophore modeling studies showed that inhibition of PUMA can be an effective strategy for mitigating radiation-induced cell death[@b49]. However, in the presence of a strong caspase/tBid positive feedback mechanism that amplifies Bax activation, even a potent inhibitor of PUMA would fail if not administered immediately after IR-exposure. In contrast, for cells with a reduced strength of caspase/tBid feedback, PUMA inhibition could be effective even if administered after substantial delay ([Figure 7D](#f7){ref-type="fig"}). The above considerations do not account for the structural promiscuity of these targets: PUMA is a BH3 domain-only protein and Bax also contains a BH3 domain. An inhibitor of PUMA may also bind and inhibit Bax, thus resulting in the effective therapy of combined I~PUMA~ + I~Bax~ ([Figure 7B](#f7){ref-type="fig"}).

It is possible to test the hypothesis (iii) above using the following two approaches. Firstly, a knockdown of PUMA when performed several hours *after* radiation damage should be ineffective in mitigating apoptosis (although when administered immediately, it would) for cells with a strong caspase/Bid feedback. Secondly, knocking down Bid mRNA several hours after radiation damage should be more effective in mitigating apoptosis at both shorter and longer time delays. Inhibitors for Bax[@b50], Bid[@b51] and PUMA[@b49] have been synthesized, that could be used for these testing purposes.

We have also performed an extensive sensitivity analysis of the network presented in [Figure 1](#f1){ref-type="fig"}, and identified robust reactions whose parameter changes do not significantly affect the cell behavior. For instance, the binding/dissociation constant between pro-apoptotic Bax and anti-apoptotic Bcl-2 (*k~25~*^±^), on the mitochondria, has a relatively small effect ([Supplementary Figure S2](#s1){ref-type="supplementary-material"}). On the other hand, our study highlights the significance of the translocation of Bax from the cytosol to the mitochondria (succeeding the transcriptional activation of Bax and PUMA by p53) for eliciting apoptotic response. The translocation unambiguously emerged from our sensitivity analysis to be a key initiator of the highly influential interactions that lead to apoptosis, indicated by *red* arrows in [Figure 1](#f1){ref-type="fig"}. This finding is in support of the recently proposed retrotranslocation of Bax from the mitochondria to the cytosol, mediated by Bcl-x~L~, as an effective mechanism for regulating pro-apoptotic activity[@b29][@b52]. Enhancing retrotranslocation may indeed be an essential strategy for decreasing susceptibility to cell death induced by IR.

An important concept that has recently emerged in cancer therapy is mitochondrial priming -- whereby cells that are exposed to higher concentrations of BH3-only proteins such as PUMA, NOXA and Bim are designated as 'well-primed\'[@b53]. Well-primed cancer cells are closer to the apoptotic threshold than are normal cells, leading to more efficient responses to therapies that induce apoptosis. Although, mitochondrial priming is used in conjunction with chemotherapy, the results in [Figure 4](#f4){ref-type="fig"} show some parallels to the priming process. Increased *k~11~* beyond a threshold value (panel **D**) presumably leads to overexpression of BH3-domain-only proteins (e.g. PUMA), promoting apoptosis. The same panel further shows that a switch from anti- to pro-apoptotic state (indicated by an abrupt change in \[C3\]) is alternatively elicited when the rate for the translocation p53^(C)^ → p53^(M)^ exceeds a threshold value (*k~10~* \> 2 × 10^−4^s^−1^), indicating that transcription-independent activity of p53 alone may also induce apoptosis. Inhibition of p53 binding to mitochondria has indeed been reported to protect mice from gamma radiation, consistent with the predicted behavior[@b54].

The development and application of systems models (whether stochastic or deterministic) to understanding cellular fate is an evolving field. We note that the sustained oscillation of p53 and Mdm2 can be maintained under deterministic conditions using certain parameter settings[@b18]. However, it is often difficult for deterministic models to reproduce the observed statistical features (e.g. distribution of oscillation periods and amplitudes) of a cell population ([Figure 3B and C](#f3){ref-type="fig"}), because these quantities tend to be constant under deterministic condition. The current model and parameters under stochastic conditions capture the cell-to-cell variability observed in single-cell experiments[@b41][@b42] while the damped oscillation behavior under deterministic condition is also consistent with the work of Alon and coworkers[@b42]. Further experimental data on changes in particular protein levels in response to radiation exposure will help refine the models, as well as establish rate parameters (via direct implementation of experimental data, or statistical inferences from expression/activity data). In the present setting, we included in our model the key events relevant to p53-mediated apoptotic response to IR. The basal activity and effects of other mechanisms such as mitochondrial fusion/fission events[@b55][@b56] are implicitly captured by the kinetic parameters. A more comprehensive study of the role of p53 dynamics would require to expand the model to include p53-regulated cell cycle arrest and DNA repair events as well mitochondria-targeted inhibitors of cyt *c* peroxidase[@b57] and caspase-independent apoptotic interactions[@b58][@b59]. These further extensions may help design more efficacious polypharmacological strategies for controlling cell susceptibility to apoptosis under different disease states and environmental challenges.

Methods
=======

Extended Gillespie algorithm
----------------------------

We adopt an extension of Direct Reaction version of the Gillespie algorithm[@b60], which takes account of the stochasticity of interactions and heterogeneity of the cellular environment[@b61]. In this approach, one reaction is chosen at a time, and the time is advanced based on the overall reaction propensity at that time. Consider the following schematic reaction,

The stochastic rate, *c*, of the reaction is related to the macroscopic kinetic rate constant, *k*, as *c = kV*, where *V* is the reaction volume; and the instantaneous propensity for this reaction is[@b60]:

where *N~X~* and *N~Y~* are the number of molecules *X* and *Y* in the reaction volume and the subscript *α* denotes that the reaction index in the system of *M* chemical reactions (1 ≤ *a* ≤ *M*). The underlying assumption is that the system is in thermal equilibrium (the number of reactive collisions are much less than that of non-reactive collisions) and well-mixed. However, cellular systems are spatially heterogeneous: for example, the mitochondria and the cytoplasm provide vastly different environments; and, there may be heterogeneities within these localizations themselves. We adopted a protocol[@b61] developed to address this issue, where the cell is subdivided into subvolumes (or compartments) and each reactant is treated as a different molecule in each subvolume, *e.g*., *X^(i)^* ≠ *X^(j)^* for the same compound *X*, where *i* and *j* refer to different subvolumes and the reaction propensities are altered accordingly. Each subvolume is assumed to be well-mixed - thus allowing for the use of the Gillespie algorithm. This protocol results in an increase in the number of reactions because (i) each reaction is treated independently in different subvolumes and (ii) additional "reactions" relating to the conversion of *X^(i)^* to *X^(j)^* (due to diffusion or translocation) are included, as described in detail by Bernstein[@b61]. The nucleus (N), cytosol (C) and the mitochondrial membrane (M) are taken as three distinctive subvolumes here. No superscripts are appended for designating compounds in the cytosol, except for clarifying ambiguous cases.

Statistical model checking
--------------------------

In order to test whether our stochastic model satisfies a given dynamical property with guaranteed confidence levels, we adopted a statistical model checking (SMC) based framework developed recently[@b43]. SMC is a highly scalable simulation-based verification approach that can check system properties encoded as logic formulas. We use a strengthened bounded linear time temporal logic (BLTL) to specify qualitative or quantitative properties of interest. Specifically, let *S* be a finite set of real-valued variables and *T* be a positive integer. An atomic proposition (*AP*) in our logic is of the form *x*\#*y*, where *x* and *y* are arithmetic expressions over real-valued variables in *S*, and \# ∈ {\>, \<, =, ≥, ≤}. The logic operators in our BLTL consist of (and), (or), (negation), and time bounded **U** (until), **G** (global), and **F** (future). The formulas of BLTL are defined as: (i) every *AP* as well as the constants *true* and *false* are BLTL formulas; (ii) if *ψ*, *ψ*\' are BLTL formulas then and are BLTL formulas; (iii) if *ψ*, *ψ*\' are BLTL formulas and *t* ≤ *T* is a positive integer then *ψ***U**^≤*t*^*ψ*\' and *ψ***U***^t^ψ*\' are BLTL formulas. The derived operators such as , **G**^≤*t*^ and **F**^≤*t*^ are defined in the usual way.

A trajectory in our model is a series of time-dependent states of the form *σ* = (*s*~0~, *t*~0~), (*s*~1~, *t*~1~),..., meaning that the system jumps to state *s~i~*~+1~ after staying in state *s~i~* for *t~i~*. We interpret the formulas of our logic at the finite set of time points **T** = {0,1,...,*T*}. A trajectory that *σ*satisfies a BLTL specified property *φ* at time *t* ∈ **T** is written as *σ*, *t*\| = *φ*. The semantics of the logic is defined as follows: (i) *σ*, *t*\| = *AP* iff *AP* holds *true* in state *s~t~*; (ii) and are interpreted in the usually way; (iii) *σ*, *t*\| = *ψ***U**^≤*k*^*ψ*\' iff there exists *k*\' such that *k*\' ≤ *k*, *t* = *k*\' ≤ *T* and *σ*, *t* + *k*\'\| = *ψ*\'. *σ*, *t* + *k*″\| = *ψ* for every 0 ≤ *k*″ \< *k*\'; (iv) *σ*, *t*\| = *ψ***U***^k^ψ*\' iff *t* + *k* ≤ *T* and *σ*, *t* + *k*\| = *ψ*\'. *σ*, *t* + *k*\'\| = *ψ* for every 0 ≤ *k*\' \< *k*.

The statements we make are in the form of *M*\| = Pr~≥*θ*~(*ψ*), meaning that the probability that the system *M* satisfies a property *ψ* is at least *θ*. For example, we express the property "caspase-3 level sustains once it reaches certain threshold" in [Figure 3D](#f3){ref-type="fig"} as follows:

The BLTL expressions of other properties can be found in the [Supplementary Text](#s1){ref-type="supplementary-material"}. Using SMC, the verification of such properties can be carried out approximately but with guaranteed confidence levels and error bounds. Specifically, we verify each property using a sequential hypothesis test between the null hypothesis H0: *p* ≥ *θ* + *δ* and the alternative hypothesis H1: *p* ≤ *θ* − *δ*, where *p* is the probability of *M* satisfying *ψ* and *δ* specifies the indifference region. The strength of the test is determined by parameter *α* and *β* which bound the type-I and type II errors, respectively. The test proceeds by generating a sequence of sample trajectories, *σ*~1~, *σ*~2~,... A corresponding sequence of Bernoulli random variables *x*~1~, *x*~2~,... are assumed, where *x~k~* = 1 if *σ*, 0\| = *ψ*, otherwise *x~k~* = 0. For each generated sample, we update the score *ω~n~* by the following function:

where *n* is the number of generated samples. We accept hypothesis H0 if *ω~n~* ≥ (1 − *β*)/*α*, and hypothesis H1 if *ω~n~* ≤ *β*/(1−*α*); otherwise, we draw another sample.

Sensitivity analysis
--------------------

Global sensitivity analysis is performed using a SMC-based multi-parametric sensitivity analysis (MPSA) method[@b43]. We encode the steady state of caspase-3 (the model output) as a BLTL formula. The MPSA procedure involves drawing a representative set of samples from the parameter space. For each sampled combination of parameter values, we compute the objective value with respect to the BLTL property. The sampled parameter sets are classified into two classes using a threshold objective value. The sensitivities are then computed as the Kolmogorov-Smirnov statistics of cumulative frequency curves associated with the two classes.

Author Contributions
====================

B.L., D.B., Z.N.O. and I.B. conceived and designed research; B.L., D.B. performed research; B.L., D.B. and I.B. analyzed the results with input from J.S.G., B.L., D.B., Z.N.O., J.S.G. and I.B. wrote the paper. All authors approved the final version of the paper.

Supplementary Material {#s1}
======================

###### Supplementary Information

Supplementary Information

Support from NIH awards P41-GM103712 (IB) and U19-AI068021 (JSG and IB) is gratefully acknowledged.

![Simplified reaction network diagram of the mathematical model.\
The diagram highlights the major reactions in the model. Basal protein synthesis and degradation reactions are included in the model but not shown. The full list of components, reactions and kinetic equations, and parameters are presented in the [Supplementary Tables S1 and S2](#s1){ref-type="supplementary-material"}. Complexes are denoted by the names of their components, separated by a dot. Single-headed solid arrows characterize irreversible reactions and double-headed arrows, reversible reactions. Dotted arrows represent enzymatic reactions. The reactions computed by sensitivity analysis ([Supplementary Figure S2](#s1){ref-type="supplementary-material"}) to play a significant role are shown by *red* arrows. Among them, the kinetic steps 20 and 10 (or associated rate constants *k~20~* and *k~10~*) lead to the respective transcription-dependent and -independent activities of p53.](srep06245-f1){#f1}

![Simulation of p53 and Mdm2 dynamics.\
(A) Comparison of stochastic (*left*) and deterministic (*right*) simulations. The time profiles of p53^(N)^ and Mdm2^(N)^ were simulated using the stochastic approach (*left* *panel*) and a deterministic approach with the same kinetic parameters (*right panel*). The stochastic simulation shows sustained oscillations while the deterministic simulation results in damped oscillations. (B) Radiation exposure is initiated at *t* = 0 (upon alteration of kinetic parameters, which applies for a duration of 12 h (*left*) or 56 h (*right*)), and leads to p53 oscillations. The amount of the radiation doses for 12 h and 56 h cases are comparable to 10 seconds and 5.6 min γ-irradiation (^60^Co, 1.8 Gy min^−1^) treatments, respectively[@b7].](srep06245-f2){#f2}

![Model predictions and validation.\
The time profiles of p53^(N)^ and Mdm2^(N)^ (A), the histograms of the period (B) and amplitude (C) of oscillations for Mdm2^(N)^ under IR are simulated and compared against previous experimental observations. *Blue* solid lines and bars depict the simulation results and *red* dots and bars indicate experimental data. The data in *blue* lines and *red* dots were normalized so that their maximum value was equal to 1. The experimental data in (A), (B) and (C) were extracted from Lahav et al, 2004[@b41], Geva-Zatorsky et al, 2006[@b7] and Geva-Zatorsky et al, 2010[@b42], respectively. (D) The dynamical properties predicted by the model are supported by experimental evidence[@b7][@b12][@b41][@b44][@b45][@b46][@b47] (type I error = 0.05, type II error = 0.05, see further details in [Supplementary Material](#s1){ref-type="supplementary-material"}).](srep06245-f3){#f3}

![Simulated apoptotic response with different strengths of coupling between p53 dynamics and mitochondrial events.\
(A) Caspase-3 time profile for different rates of p53 transcriptional activation of pro-apoptotic proteins (represented by *k~11~*). (B) Caspase-3 time profile for different p53 translocation rate from the cytoplasm to the mitochondria (represented by *k~10~*). (C) The simulated coupling strength-caspase activation response curves. (D) Caspase-3 time profile under low frequency (of p53 pulses) conditions. This condition is denoted as LF, while normal frequency conditions are denoted as WT. The associated parameters can be found in [Supplementary Table S3](#s1){ref-type="supplementary-material"}.](srep06245-f4){#f4}

![Robustness of caspase-3 activity with respect to association strength of Bax^(M)^ and Bcl-2 in the mitochondria.\
Panels A and C show onset of apoptosis (*blue* curve, for C3; *right* ordinate) despite vastly differing dissociation constants, *k~25~^−^/k~25~^+^*, for the Bax^(M)^.Bcl-2 complex formed in the mitochondria (*k~11~* = 2 × 10^−5^ s^−1^ in these two panels). Note the different scales for the \[Bax^(M)^.Bcl-2\] profiles in panels A and C (*black* curve; *left* ordinate). Likewise, panels B and D show the insensitivity of the system behavior (low C3 levels) to *k~25~^−^/k~25~^+^*, under conditions conducive to cell survival (*k~11~* = 10^−5^ s^−1^).](srep06245-f5){#f5}

![Role of tBid activation for sustained caspase-3 activity and effect of the strength of caspase/Bid positive feedback loop in mediating apoptosis and its inhibition.\
(A--B) Time evolution of the concentrations of PUMA, caspase-3, and tBid at in response to long time IR exposure predicted for two different coupling strengths of p53 transcriptional activation to mitochondrial pore opening events (represented by *k~11~*). PUMA exhibits oscillations echoing the behavior of p53. tBid production remains limited (*red* curve) on the *left* panel A (low *k~11~)* such that caspase-3 levels are negligibly small. Enhanced coupling (*right* panel) yields a sustained increase in tBid levels, which in turn induces a significant increase in caspase-3 levels, via a positive feedback loop. Note the ordinate scale difference between the two panels. (C) The caspase-3 time profiles were simulated for three *k~28~*, which correspond to different strengths of caspase/Bid feedback loop. Upon decreasing *k~28~*, caspase-3 concentration progressively decreases. (D) Feedback loop strength-caspase activation response curves.](srep06245-f6){#f6}

![Potential mitigation of radiation damage via individual and combination therapies.\
Panels (A) and (B) show the effects of therapies based on single targets, PUMA, Bid, caspase-3, Bax, or their combination when inhibitors are administered at 15 min (A) and 12 h (B) after radiation exposure. Bid and Bax appear to be the most effective targets even when inhibited with a time delay. The control (no inhibition) is shown in both panels. The efficacy of PUMA inhibition strongly depends on the timing of the drug treatment. (C) Time evolution of caspase-3 in the absence of inhibitors (*black*) and in the presence inhibitors of PUMA and C3 (*pink*), and PUMA and Bax (*blue*). (D) A reduction in caspase/tBid feedback strength (*k~28~*) (from 0.5 to 0.3 μM^−1^s^−1^) rescues the efficacy of PUMA inhibitor administrated after a delay of 12 h.](srep06245-f7){#f7}
