STAT1–p53 Dynamics Program Cell Fate through p21 and PUMA

STAT1–p53 Dynamics Program Cell Fate through p21 and PUMA | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results STAT1–p53 Dynamics Program Cell Fate through p21 and PUMA View ORCID Profile M.K. Francis Didier Tshianyi , Rebecca Walo Omana , Apollinaire Ndondo Mboma , Didier Kumwimba , View ORCID Profile Didier Gonze doi: https://doi.org/10.1101/2025.10.28.685172 M.K. Francis Didier Tshianyi a Unit of Theoretical Chronobiology, Faculty of Sciences, Université libre de Bruxelles (ULB) , Brussels, Belgium c Department of Mathematics and Computer Science, Faculty of Sciences, University of Lubumbashi , Lubumbashi, DR Congo Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for M.K. Francis Didier Tshianyi For correspondence: francis.tshianyi.mwana.kalala{at}ulb.be tshiakal{at}gmail.com Rebecca Walo Omana b Department of Mathematics and Computer Science, Faculty of Science and Technology, University of Kinshasa , Kinshasa, DR Congo Find this author on Google Scholar Find this author on PubMed Search for this author on this site Apollinaire Ndondo Mboma c Department of Mathematics and Computer Science, Faculty of Sciences, University of Lubumbashi , Lubumbashi, DR Congo d Department of Science and Technologies, Faculty of Sciences and Technologies, New Horizon University , Lubumbashi, DR Congo Find this author on Google Scholar Find this author on PubMed Search for this author on this site Didier Kumwimba c Department of Mathematics and Computer Science, Faculty of Sciences, University of Lubumbashi , Lubumbashi, DR Congo Find this author on Google Scholar Find this author on PubMed Search for this author on this site Didier Gonze a Unit of Theoretical Chronobiology, Faculty of Sciences, Université libre de Bruxelles (ULB) , Brussels, Belgium Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Didier Gonze Abstract Full Text Info/History Metrics Preview PDF Abstract Signalling pathways tightly regulate stress-induced cell fate decisions, with p 53 stabilization via attenuation of the p 53– Mdm 2 feedback being central to effective responses. Signal Transducer and Activator of Transcription 1 (STAT1) both represses Mdm2 transcription and co-activates p53 targets ( p21 and PUMA ), yet these dual roles have not been quantified within a unified framework. We build a mechanistic model that couples a biologically calibrated stress scale ( S ) and STAT1 activity (Σ) to the canonical p 53– Mdm 2 core and to downstream p21 and PUMA modules. Key processes are not described by fixed parameters but by explicit functions of S , allowing the model to self-adjust across stress levels without manual re-fitting. This formulation ensures that oscillatory, damped, or plateau responses naturally emerge by varying S alone. STAT1 acts as a dynamical gain on the core (controlling mean level, pulse amplitude, and duration) and as a transcriptional co-activator scaling promoter strengths. Bifurcation analysis reveals a two-dimensional Hopf region in the ( S , Σ) plane; increasing Σ shifts this region to higher S and progressively narrows it, ultimately quenching oscillations at large Σ . Simulations of a generic cell with time-varying S and Σ, combined with a two-stage decision rule (early transient detection followed by stationary readout), map p 53 dynamics to fate: sustained moderate oscillations align with arrest, damped intermediate responses with senescence, and strongly damped high plateaus with apoptosis. The model reproduces cell fate distributions reported in literature for different cell lines (MCF-7, HCT116, U2OS) without kinetic parameter re-fitting, and highlights cell-type-specific sensitivity to p21 versus PUMA . Our framework identifies STAT1 as a tunable amplifier and oscillation quencher of the stress-responsive p53 network, providing testable dynamics-based predictions for fate control. Introduction Cellular stress can profoundly reshape protein dynamics. It arises from diverse agents (UV, γ -irradiation, viral infection…) and perturbs genome integrity and homeostasis. To cope with these adverse conditions, cells deploy signalling networks in which the tumour suppressor p53 is central. Under basal conditions, p 53 is kept low by its principal regulator Mdm2 , which binds and ubiquitinates p 53, promoting nuclear export or proteasomal degradation; this forms a negative feedback in which p 53 induces Mdm2 , and Mdm2 accelerates p 53 turnover ( Bar-Or et al., 2000 ; Lahav et al., 2004 ; Ciliberto et al., 2005 ; Proctor & Gray, 2008 ). A hallmark of this module is its ability to generate oscillations in response to DNA damage, observed experimentally and captured by prior models ( Abou-Jaoudé et al., 2009 ; Ouattara et al., 2010 ; Batchelor et al., 2011 ). Beyond the canonical DNA-damage axis, we consider the transcription factor STAT1. Biologically, STAT1 is activated by cytokines (notably type I/II interferons) via JAK/STAT signalling; upon activation it dimerizes, translocates to the nucleus and modulates gene expression ( Levy & Darnell, 2002 ; Giles et al., 2017 ). Importantly, studies indicate that STAT1 can stabilise p 53 and cooperate to regulate the expression of p 53 target genes ( Brzostek-Racine et al., 2011 ; Yu et al., 2015; Heitmeier et al., 1999 ). In practical terms, this gives two conceptually distinct inputs: a stress/DNA-damage signal (e.g. ATM/Chk-mediated) and a cytokine/STAT1 signal (e.g. interferon bursts). These inputs can combine in time—for instance, after irradiation the effective DNA damage may decrease as lesions are repaired while paracrine interferon induces a transient STAT1 peak ; conversely, inflammatory contexts can sustain STAT1 activity with minimal ongoing damage. It is not straightforward to predict how different mixtures of stress and STAT1 activity, acting through p 53, will shape the dynamics of downstream protein effectors, namely p21 and PUMA , nor how cells interpret these signals to choose between cell-cycle arrest, senescence or apoptosis . Although biology suggests that STAT1 represses Mdm2 and enhances p 53 target programmes ( Brzostek-Racine et al., 2011 ; Feng et al., 2020 );, to our knowledge, no mathematical studies have explicitly integrated STAT1 with the p 53– Mdm 2 oscillator and its downstream effectors. This gap underlies the need for a kinetic model. We develop a mechanistic model of the canonical p 53– Mdm 2 core that explicitly takes the stress level S and STAT1 activity Σ as inputs and propagates their effects to the two key effectors p21 and PUMA . In the model, STAT1 both represses Mdm2 transcription (stabilising p 53) and co-activates p 53 targets ( p21, PUMA ). We introduce a two-stage fate criterion: (i) an early transient detector marking commitment once p21 / PUMA exceed thresholds for a minimal dwell time, and (ii) a late stationary assessment within 60–72 h using the same threshold-based rules. A bifurcation analysis reveals a two-dimensional oscillatory region delimited by Hopf bifurcations in the ( S , Σ) plane that organises transitions between oscillatory and plateau regimes of p 53. In a single “standard” cell under time-constant S , increasing Σ precipitates the early detector and shortens decision times, biasing outcomes toward the PUMA -dominated (pro-apoptotic) branch; low–moderate Σ favours p21 -driven arrest/senescence. At the population level, we apply the same rule—without re-fitting biochemical parameters—to archetypes cell lines (MCF-7, HCT116, U2OS), obtaining literature-consistent fate splits and revealing cell-type-specific sensitivity to p21 vs PUMA ( Abbas & Dutta, 2009 ; Jeffers et al., 2003 ). Throughout, qualitative p 53 motifs map to fate tendencies (sustained moderate-mean oscillations are linked to cell cycle arrest; damped intermediate responses are linked to senescence; strongly damped high plateaux are linked to apoptosis), providing testable predictions under time-varying S and Σ. Most mathematical models of the p 53– Mdm 2 module assume fixed parameters. To reproduce different behaviours under low, moderate, or high stress, one typically has to re-tune rate constants or re-fit parameter sets for each condition (see e.g. ( Bar-Or et al., 2000 ; Proctor & Gray, 2008 ; Ciliberto et al., 2005 ; Eliaš et al., 2014 )). By contrast, our framework is intrinsically stress-modulated : the input S explicitly shapes several key processes through functional dependencies (Hill-type or saturating terms). As a consequence, no additional manual adjustment is required to generate the repertoire of dynamical responses. Our model allows to predict cell fate under various profiles of stress and Stat1. We integrate a biologically calibrated stress scale with a dual-role STAT1 modulation of the p 53– Mdm 2 network, introduce a decision rule that scales from single cells to heterogeneous populations, and reveal a two-dimensional Hopf region in ( S , Σ) that highlights oscillation–plateau transitions. Our main contributions are: Integrating the dual-role STAT1 in one mechanistic framework: repression of Mdm2 (stabilizing p 53) and co-activation of p21 / PUMA . Identifying the two-dimensional Hopf region in ( S , Σ)that shifts rightward and narrows with Σ, reshaping oscillatory domains of p 53. Implementation of Stress-modulated parameters that embed S directly into the system’s kinetics, eliminating the need for manual re-tuning across stress levels. Defining a two-stage fate rule (early detector + late stationary readout) mapping p 53 dynamics to arrest/senescence/apoptosis, with thresholds and minimal dwell times. A population-level validation across archetypes (MCF-7, HCT116, U2OS) without biochemical re-fitting , reproducing literature distributions and revealing sensitivity to p21 versus PUMA . Section 1 introduces the model and the assumptions used; Section 2 presents the bifurcation analysis and p 53 dynamical patterns; Section 3 formalises the two-stage fate rule and illustrates their application on single-cells; Section 4.1 introduces the archetypes, their calibration targets, and demonstrates the two-stage fate rule on them. Section 4.4 discusses biological implications, limitations, and testable predictions. 1. The model In this section we present the hypothesis and the different processes we used to construct our model. We consider a minimal p 53– Mdm2 module driven by two external inputs : the DNA-damage/stress signal S ( t ) and the STAT1 activity Σ( t ). In this module they are treated as prescribed forcings (constants or time profiles), i.e. as parameters rather than state variables. p53 balance. p 53 is produced at a basal rate β 0 and is further induced by stress through a sigmoidal Hill term of order 4, . Removal of p 53 by Mdm2 is described by a Michaelis–Menten term , multiplied by so that increasing S weakens the effective ubiquitylation capacity and stabilises p 53. A first–order term δ p [ p 53] accounts for basal turnover. Mdm2 RNA and protein . Transcription of Mdm2 mRNA is activated by p 53 with apparent cooperativity (quadratic dependence) and repressed by STAT1, where the quadratic [ p 53] 2 term follows ( Hunziker et al., 2010 ), who used it to capture the effective cooperativity/tetramerisation of p 53 in driving Mdm2 expression. The mRNA decays at rate δ m . Mdm2 protein is translated at rate k tl [ Mdm2 RNA ] and degraded with rate δ M [ Mdm2 ]. The model reads: Throughout the manuscript, square brackets denote concentrations . For any molecular species X , we write [ X ] for its concentration (typically in nM or µ M), e.g. [ p 21], [ PUMA ], [ p 53]. When time dependence is relevant we use [ X ]( t ). For readability, the ODE state variables may be written without brackets (e.g., p 21) but are understood to represent the same concentration quantity with the same units unless explicitly stated otherwise. 1.1 Stress scale Stress ( S ) is a single, dimensionless control parameter that modulates all stress–sensitive processes in the model, including the p 53 half-life via ubiquitination, as well as the production of p 53, p 21, and PUMA. To enable comparison across experiments, we calibrated a simple scale S ∈ [0, 25] from published responses to UV and γ -irradiation (1). These results were calibrated primarily based on: (i) non-transformed murine cells (UV); (ii) human skin in vivo (UV); and (iii) the human cell lines MCF-7, MCF10A, and RPE-hTERT for fine measurements of dynamics and half-life (UV and γ ). In the manuscript, S = 0 denotes no insult and larger values encode stronger damage. We will use the term generic cell or standard cell to refer to any cell that reproduces the behaviour described in Table 1 . View this table: View inline View popup Download powerpoint Table 1. Calibrated stress scale S : correspondence with benchmark γ -irradiation doses, qualitative biological effects, and p 53 dynamics. UV exposures were mapped to S by matching p 53 stabilization metrics (see text). View this table: View inline View popup Download powerpoint Table 2. Parameters for the p53–Mdm2–p21–PUMA model. In Ref , estimated indicates no published numeric interval was found that contains the table value. The parameter values given here account for a generic cell. 1.2 Stat1 and stress In this article Σ( t ) denotes an interferon (IFN)–driven JAK–STAT signal that activates STAT1 (Y701 phosphorylation, dimerization, nuclear translocation). Σ( t ) captures the cytokine context (paracrine/autocrine IFN- α/β or IFN- γ ) that co-activates parts of the p 53 program, including PUMA and p 21 ( Levy & Darnell, 2002 ; Giles et al., 2017 ; Kalliara et al., 2022 ; Yu et al., 2015; Abbas & Dutta, 2009 ; Nakano & Vousden, 2001 ). A single insult (irradiation, UV, some chemotherapies) can simultaneously activate the p 53 arm via DNA damage and, through transient interferon bursts, the JAK–STAT pathway and STAT1; in such coupled cases, S and Σ rise together but with distinct kinetics (often a delayed Σ) ( Bar-Or et al., 2000 ; Lahav et al., 2004 ; Proctor & Gray, 2008 ; Batchelor et al., 2011 ). Conversely, cytokine exposure or infection can elevate Σ with little DNA damage, while sterile genotoxic stress can raise S with minimal interferon signalling ( Vousden & Ryan, 2009 ). These scenarios motivate treating S ( t ) and Σ( t ) as independent, externally controlled inputs and exploring their combinations/timings to reveal how they reshape p 53 dynamics and the downstream balance between p21 and PUMA. 1.3 p21 and PUMA p21 ( CDKN1A ; p21 Cip1 / Waf1 ) is the cyclin–dependent kinase inhibitor CDKN1A , historically referred to as p 21 Cip1 / Waf1 (also SDI1). It promotes cell–cycle arrest and, under sustained activation, senescence. In the p 53 network, p 21 acts as a relatively slow , integrative branch: its mRNA/protein turnover is slower than immediate–early responders, so its level reflects the time–integrated history of p 53 activity ( El-Deiry et al., 1993 ; Harper et al., 1993 ; Xiong et al., 1993; Abbas & Dutta, 2009 ; Barr et al., 2017 ; Jin et al., 2003 ; Coleman et al., 2003 ) . PUMA stands for p53 Upregulated Modulator of Apoptosis . It is a p 53–induced BH3–only protein that antagonizes anti–apoptotic BCL–2 family members. Functionally, PUMA is a fast , switch–like effector: it can rise quickly following strong p 53 activity and directly commits cells to apoptosis when its threshold is crossed ( Nakano & Vousden, 2001 ; Jeffers et al., 2003 ; Purvis et al., 2012 ; Hanson et al., 2019 ; Yu et al., 2015) . We extend the core p 53– Mdm 2 module in (1.1) with a branch for p 21 and PUMA. p21 (cell-cycle arrest/senescence) . To capture p 21’s early engagement under mild damage, stress S ( t ) enhances transcription through a low-threshold sigmoidal term, so that half-activation occurs at small S ( Abbas & Dutta, 2009 ; Stewart-Ornstein & Lahav, 2017 ; Barr et al., 2017 ). At stronger stress, multiple post-transcriptional/post-translational checks are known to limit p 21 (e.g., altered stability and proteasomal turnover) ( Jin et al., 2003 ; Coleman et al., 2003 ; Hanson et al., 2019 ; Schwanhäusser et al., 2011 ). We therefore include a stress-dependent translational factor which is near-maximal at intermediate S but attenuates p 21 production at very high S , reflecting the reported dampening of p 21 under strong stimuli. Regarding STAT1, we allow only a modest co-activation of p 21, consistent with literature indicating that p 21 is primarily governed by p 53 at low–moderate stress while STATs can provide auxiliary input ( Levy & Darnell, 2002 ; Abbas & Dutta, 2009 ; Yu et al., 2015). PUMA (apoptosis) . PUMA is a pro-apoptotic p 53 target that is robustly induced under genotoxic stress ( Nakano & Vousden, 2001 ; Jeffers et al., 2003 ). In contrast to p 21, we posit a higher transcriptional threshold for stress, so that PUMA half-activation occurs at stronger S , aligning with its role in apoptosis ( Purvis et al., 2012 ; Stewart-Ornstein & Lahav, 2017 ). To model stress resilience of translation, we include a gain motivated by stress-sensitive translational control observed in the p 53 network and proteome-wide measurements ( Schwanhäusser et al., 2011 ; Hanson et al., 2019 ). Unlike p 21, PUMA is assumed to be more sensitive to STAT1: we multiply its promoter gain by a steeper co-activation factor of Σ (e.g. a higher Hill sensitivity), consistent with interferon/STAT1 signalling intersecting pro-apoptotic p 53 programmes ( Levy & Darnell, 2002 ; Giles et al., 2017 ; Brzostek-Racine et al., 2011 ). The extended model is given by: Where: The following table gives the parameters values and their references: 2. Stability and bifurcation In this section, we characterise both the dynamical behaviour and the mathematical properties of the model that couples stress ( S ) and STAT1 activity (Σ) to the canonical p 53– Mdm 2 core and its downstream effectors p21 and PUMA . 2.1. Stress–dependent dynamics of the p53–Mdm2 core We first examined how the core module responds when only the stress input S is varied (at fixed STAT1 activity Σ = 100). Figure 2 illustrates the time courses of p 53, Mdm 2 RNA and Mdm 2 protein for representative levels S ∈ {2, 10, 18, 25} . All simulations are initialized close to the low-stress steady state of the model. Specifically, each species starts from a small but non-zero basal concentration ( p 53 = 10 − 3 , Mdm 2_ RNA = 10 − 7 , Mdm 2 = 1.5 × 10 − 7 , PUMA = 10 − 6 , p 21 = 10 − 6 ). These values approximate the quasi-equilibrium reached in the absence of stress ( S ≈ 0), where p 53 is maintained at a low level by efficient Mdm2-mediated degradation, and both PUMA and p21 remain transcriptionally silent. They prevent numerical artifacts associated with strictly zero initial conditions while preserving the biological interpretation of a basal, unstressed cell. Unless otherwise stated, these initial conditions are used for all single-cell simulations and results presented in this work. Download figure Open in new tab Figure 1. Scheme of the model: Red dashed labels (“Stress”) denote the stress input S ( t ), which stabilises p 53 and modulates downstream processes; purple dashed labels (“Stat1”) denote interferon/STAT1 activity Σ( t ). p 53 activates Mdm2 mRNA transcription, and Mdm2 protein induces p 53 degradation. STAT1 has a dual role: it represses Mdm2 transcription (bar placed on the p 53 → Mdm2 mRNA link) and co-activates the p 53 programme toward p21 (cell-cycle arrest) and PUMA (apoptosis). Arrows indicate activation; blunt arrows indicate inhibition; dashed connections denote regulatory influences (not material flow). The triangular wedge on the p21 mRNA → p 21 arrow encodes a stress-dependent translation gate: at low S , the effective translation rate decreases as S increases, and it shuts off for extreme stress levels. Download figure Open in new tab Figure 2. Stress scan at fixed Σ = 100. Time courses of p 53, Mdm 2 RNA and M dm 2 protein for S ∈ {2, 10, 18, 25} . Low S : basal fixed point; intermediate S : sustained oscillations (period ∼ 6–8 h); high S : damped transients to an elevated plateau. We have three dynamical regimes For Low stress (e.g. S ≲ 4; S = 2 shown). The system remains close to basal: p 53 is quickly cleared by Mdm 2 and trajectories converge to a low steady state with negligible oscillations. At Intermediate stress (oscillatory band; S = 10). We observe sustained p 53 oscillations, mirrored by delayed oscillations of Mdm 2 RNA and protein. The mean p 53 level rises with S , the peak–to–peak amplitude is maximal in this band, and the inter-peak interval (IPI) is about 6.1 h. At High stress (damped/plateau; S = 18–25). A large overshoot is followed by damped oscillations that settle to an elevated plateau; for very high stress ( S = 25), the response is essentially monotone. Here, the rising Mdm2 synthesis rapidly restores control, suppressing sustained oscillations. Effect of S on p53 period, amplitude, and peak count The global sensitivity ( SRC , standardized regression coefficient, final–mean), presented in Section S2.2 of the Supplement, quantifies the normalized linear influence of each parameter on a given output variable across simulations. The SRC analysis is consistent with these trends: the period increases with S (positive SRC), the peak count decreases (negative SRC), and the stationary amplitude shows a weaker, overall negative dependence on S —in line with the large initial transient followed by diminished sustained oscillations at high stress. A bifurcation analysis will pinpoint the Hopf threshold S * that separates the pulsatile regime from the damped, high-plateau regime. Mechanistic interpretation Stress acts on two levers of the core loop: (i) it increases p 53 production through the sigmoidal term ; and (ii) it weakens Mdm 2-mediated removal via the factor . The core module exhibits behaviours like the one summarized in Table 1 . 2.2. Existence, positivity, and boundedness of solutions We consider the full five–dimensional system defined on , with strictly positive parameters and initial conditions. Here g 21 and g pu denote bounded Hill-type functions satisfying: 0 ≤ g 21 , g pu ≤ 1. Theorem: Positivity and boundedness of the full model Theorem: Positivityl model For any nonnegative measurable inputs S ( t ), Stat 1( t ) and strictly positive initial conditions, the system above admits a unique global solution ( p 53, Mdm 2 RNA , Mdm 2, p 21, PUMA )( t ) such that: The positive cone is invariant (all variables remain nonnegative). The trajectories are uniformly bounded in time: there exist explicit constants C p 53 , C mRNA , C Mdm 2 , C p 21 , C PUMA depending only on the parameters such that p 53( t ) ≤ C p 53 , Mdm 2 RNA ( t ) ≤ C mRNA , Mdm 2( t ) ≤ C Mdm 2 , p 21( t ) ≤ C p 21 , PUMA ( t ) ≤ C PUMA , ∀ t ≥ 0. Sketch of proof. Positivity At each boundary face (e.g. p 53 = 0, Mdm 2 RNA = 0, p 21 = 0, PUMA = 0), the corresponding right-hand side is nonnegative, so that the positive cone Λ is invariant. Bound for p53 (by the Grönwall–Bellman inequality) Since and , we obtain By the differential Grönwall inequality, Bound for p21 (Grönwall) The Hill terms satisfy and , while g 21 ( S ) ≤ 1. Hence By Grönwall, Other variables For Mdm 2 RNA , Mdm 2 and PUMA , one uses the previously obtained bound on p 53 together with the fact that the Hill functions are bounded by 1. Each equation is of the form u ′ ( t ) ≤ A − Bu ( t ) with A, B > 0, so the Grönwall–Bellman inequality yields uniform bounds (full details and explicit constants C mRNA , C Mdm 2 , C PUMA are given in the Supplement). Moreover, since the right-hand side is 𝒞 1 in the state variables on Λ, the Picard–Lindelöf theorem ensures existence and uniqueness of a local solution. This establishes that the problem is well posed both mathematically and biologically. Conclusion The entire five-dimensional model is well-posed: all species concentrations remain nonnegative and bounded for all time. The Grönwall–Bellman inequality’s provides explicit bounds for p 53 and p 21; the other cases are analogous and deferred to the Supplement for clarity. 2.3. Stability No bistability For physiological consistency and a history–independent response, we must ensure that the system does not admit bistability. This is guaranteed by the following results: Equilibria (characterization) .. Any equilibrium ( p 53 * , Mdm 2 * RNA , Mdm 2 * , p 21 * , PUMA * ) satisfies with p 53 * solving Theorem: Uniqueness of the equilibrium point Consider the subsystem defined on with strictly positive parameters and initial data. Then, for any S > 0, the system admits a unique strictly positive equilibrium, obtained as the unique solution of where (a positive constant), and F is strictly decreasing. Hence, there is neither bistability nor multiplicity of equilibria for any S . Sketch of proof At equilibrium one has and , so the p 53 steady-state equation reduces to F ( p 53 * ) = 0 with Thus F is strictly decreasing on (0, ∞). Since and F ( p 53) → −∞ as p 53 → ∞, the intermediate value theorem yields a unique root p 53 * > 0. Back–substitution gives a unique . □ Why this suffices for the full model . In the full five-dimensional model, every equilibrium takes the form i.e., all components are functions of p 53 * . In particular, there can be no equilibrium multiplicity and hence no classical bistability of equilibria. 3D phase portrait and stability We fix STAT1 at Σ = 100 and track the Jacobian spectrum at the global equilibrium as stress S varies. The real part of the dominant complex pair crosses zero at S * ≈ 4.32, indicating a supercritical Hopf bifurcation ( Fig. S2 ). For S S * , a stable limit cycle emerges. At Σ = 100 ( Fig. 3 ), each column corresponds to a fixed value of S , while the rows display (Mdm2_RNA, Mdm2, p 53), (Mdm2, p 53, p 21), and (Mdm2, p 53, PUMA). Tight inward spirals indicate convergence toward a stable fixed point, whereas tubular coils correspond to stable limit cycles characterized by the expected phase lags, with p 53 leading Mdm2 and p 21/PUMA responding with additional delay. Download figure Open in new tab Figure 3. 3D phase portraits for p 21 and PUMA at S = 2, 4.32, 12, 24 with Σ = 100. At S = 2 (low stress), the trajectories rapidly relax to a low fixed point: p 53 is efficiently cleared by Mdm2 (reflecting high k cat values), and both p 21 and PUMA remain far below their activation thresholds. This regime corresponds to a cycling, unstressed state. For S = 4.32 (near the Hopf bifurcation), a small-amplitude limit cycle emerges. p 21 exhibits weak, delayed pulses, while PUMA remains negligible. The system shows transient oscillations consistent with a reversible arrest without full commitment. For S = 12 (robust oscillations), the trajectories form larger and slower cycles. p 21 pulses are strong and regular, whereas PUMA remains below its threshold. This configuration reflects a sustained arrest state, which may progressively evolve toward senescence if p 21 exceeds higher cutoffs over time. At S = 24 (extreme stress), the dynamics collapse onto a high fixed point as the feedback weakens (due to the reduction of k cat at large S ). In this regime, p 21 expression is suppressed, while PUMA rises toward a high plateau, consistent with an early apoptotic commitment. 2.4. Oscillatory window and bifurcation landscape in the (S, Σ) plane Figure 4 summarizes the dynamical structure of the p 53– Mdm 2– p 21–PUMA network across the ( S , Σ) plane. The maps combine steady and oscillatory readouts: the first row reports maximal concentrations of p 53, p 21, and PUMA, while the second row shows their one–sided amplitudes evaluated on a late-time window (60–72 h). Superimposed dotted and dash–dotted curves indicate the numerically determined Hopf loci and , which delimit the region of sustained oscillations. Download figure Open in new tab Figure 4. Oscillatory window and bifurcation landscape in the ( S , Σ) plane. (A1–A3) Maximal concentrations and (B1–B3) one–sided amplitudes of p 53, p 21, and PUMA. Dotted and dash–dotted curves show the Hopf loci and . Increasing Σ shifts the oscillatory domain rightward and compresses it, eventually collapsing the pulses into a steady regime. As Σ increases, the onset (dotted) shifts rightward and the offset (dash–dotted) moves leftward, progressively narrowing the oscillatory band. For Σ ≳ 3.2 × 10 2 , the two boundaries nearly coincide and the system becomes aperiodic for all S . This behaviour reflects a decreasing effective loop gain in the p 53→ Mdm 2 feedback: STAT1 enters through the regulatory factor in the Mdm 2 RNA equation, so increasing Σ reduces Mdm 2 production and damps the feedback. At low Σ, moderate stress destabilises the fixed point (Hopf onset), whereas at higher S the stress–dependent term k cat ( S ) decreases the Mdm 2-mediated p53 degradation rate, restoring stability (Hopf offset). Together, these mechanisms explain the closed oscillatory window bounded by and . Figure S2 in the supplement illustrates Σ effect on bifurcation values. Maximal concentrations (A1–A3). p 53 reaches its highest levels along the upper-right flank of the Hopf region, reflecting the cumulative accumulation under sustained stress. p 21 displays a more selective pattern: high expression is confined to a narrow intermediate- S domain due to the combined stress gating and the biphasic modulator f hill ( S ). PUMA, in contrast, increases almost monotonically with both S and Σ, dominated by the Σ-dependent co-activation κ 1 (Σ); its behaviour is only weakly influenced by the Hopf contours. Oscillation amplitudes (B1–B3). p 53 amplitude rises sharply at , peaks inside the Hopf domain, and vanishes at . p 21 follows a similar but narrower ridge and decays more rapidly with Σ, mirroring the attenuation of the upstream feedback. PUMA exhibits only minor oscillatory modulation, confirming that it integrates the steady (non-pulsatile) component of the p 53 signal. The combined analysis delineates stress–response zones in the ( S , Σ) plane. Moderate S and intermediate–high Σ values place cells near the p 53 /p 21 amplitude ridge—favouring adaptive or senescent responses—whereas large S or Σ drive PUMA accumulation and suppress oscillations, biasing toward apoptosis. Thus, STAT1 acts as a bifurcation modulator that reshapes the oscillatory landscape of the p 53 network, bridging dynamical stability with cell-fate transitions. In addition to the stress level S , the model also exhibits bifurcation behaviour with respect to two other parameters: the saturation constant K 0 controlling Mdm2-mediated degradation of p53, and the degradation rate δ M of Mdm2. As shown in Fig. S4 (Supplementary Information), both parameters can delimit regions where Hopf bifurcations occur. Importantly, increasing STAT1 levels Σ modulates the extent of these bifurcation zones, but the effect is less pronounced than for the primary bifurcation parameter S . This indicates that while K 0 and δ M shape the stability boundaries of the p53 oscillator, the dominant control of oscillatory dynamics remains exerted by the stress input. A 3D version for Figure 4 is given in the supplementary Fig. S3 . Together, these results establish the dynamical consistency of the model. The p 53–Mdm2 core admits a single equilibrium for any fixed stress level S , ruling out bistability and ensuring a history–independent response. As S increases, this equilibrium successively loses and regains stability through a pair of Hopf bifurcations that delimit a well-defined oscillatory window. Within this window, the feedback loop generates self-sustained p 53 pulses whose amplitude and frequency are modulated by both stress and STAT1 levels. Beyond it, trajectories converge either to a low basal state (weak stress) or to a high steady state (extreme stress), corresponding respectively to physiological homeostasis and apoptotic commitment. Hence, the analysis of stability and bifurcation provides a rigorous foundation for interpreting the transition from stable to oscillatory regimes and, ultimately, for linking dynamical behaviour to cell-fate outcomes. 3. Cell-fate decision: predicting apoptosis or cell-cycle arrest from the model In complex signalling systems such as the p53–Mdm2–p21–PUMA network, cellular outcomes emerge from the interplay of multiple regulatory inputs. Here, both the genotoxic stress level ( S ) and the STAT1 activity (Σ) act as upstream modulators, dynamically shaping the p53 response. While qualitative patterns (e.g., oscillatory versus steady behaviour) can often be related to stress intensity, predicting the ultimate cell fate—whether adaptive, senescent, or apoptotic—remains challenging when two or more control axes interact. As illustrated below, even for fixed parameter sets, distinct combinations of ( S , Σ( t )) may yield temporally complex p 21 and PUMA profiles whose shape alone does not immediately reveal the underlying decision. This motivates the development of a formal, quantitative decision criterion based on measurable dynamical thresholds, enabling systematic classification of single-cell trajectories into distinct fate domains. We define a generic cell-fate criterion under stress , applicable to any single cell whether stress S and STAT1 activity Σ are constant or time-varying. In other words, the rule applies equally to arbitrary trajectories S ( t ) and Σ( t ) and to the particular cases S ( t ) ≡ S and/or Σ( t ) ≡ Σ. Several studies have attempted to extract fate predictors from p53 dynamics, relating pulsatile versus sustained responses to transient arrest, senescence, or apoptosis ( Lahav et al., 2004 ; Geva-Zatorsky et al., 2006 ; Batchelor et al., 2011 ; Purvis et al., 2012 ; Stewart-Ornstein & Lahav, 2017 ). Quantitative descriptors such as the number of p53 pulses n p53 , their amplitude A p53 , or the cumulative exposure have been used to classify outcomes ( Purvis et al., 2012 ). Downstream markers provided complementary criteria: p 21 thresholds θ 21 or time-integrated levels were associated with arrest and senescence ( El-Deiry et al., 1993 ; Xiong et al., 1993; Abbas & Dutta, 2009 ; Barr et al., 2017 ; Hanson et al., 2019 ), while strong and sustained PUMA induction signalled apoptotic commitment ( Jeffers et al., 2003 ; Nakano & Vousden, 2001 ; Sun et al., 2009 ). More integrative formulations compared these effectors directly, using ratios such as R 21 = ⟨ p 21⟩ late / max t p 21( t ) ( Stewart-Ornstein & Lahav, 2017 ) or ρ = ⟨ PUMA ⟩ / ⟨ p 21⟩ ( Kracikova et al., 2013 ), where large ρ values indicate apoptosis and small ones stable arrest. Although these indices capture important qualitative trends, they remain one-dimensional and context-specific. In contrast, our criterion introduces a two-dimensional and temporally resolved mapping from (⟨ p 21⟩, ⟨ PUMA ⟩) to discrete fate domains, defined by quantitative thresholds µ 21 and µ PU evaluated over early and stationary windows. This formulation unifies dynamic, molecular, and bifurcation-based perspectives into a continuous decision landscape that can be applied to heterogeneous cell populations and compared across archetypes. Early decision (before 60 h) . We identify stable activation patterns that persist for at least the minimal required duration. If PUMA remains continuously high while p21 stays below its upper senescence threshold over a given time , apoptosis will be triggered. If p 21 lies in the arrest band (between and ) with PUMA for , the cell cycle will be halted. If p 21 lies in the senescence band (between and ) with PUMA low for , or is stably very high for , the cell will enter senescence. Early index (PUMA–p21 overlap) . When both PUMA and p21 rise together, we look for a time window where PUMA stays above its early threshold and p21 exceeds at least one of its activation thresholds . Within that window we compute a balance index comparing peak PUMA to the sum of both peaks: We then classify using two global cutoffs : Cycling is excluded here, since both pathways are already engaged . If no qualifying window exists, no early decision is taken. Stationary decision (60–72 h) . At the end of the experiment we average PUMA(⟨ PUMA ⟩) and p21(⟨ p 21⟩) over the late window (60–72 h) and compare them to global thresholds. If both are low (below µ 21 and µ PU ), cell cycle will continue. If p21 is very high (above µ S ) with PUMA low , cell enters in senescence. If p21 exceeds its decision level (≥ µ 21 ) with PUMA low , cell cycle will be halted. If PUMA is high (≥ µ PU ) while p21 remains below its decision level , apoptosis will be triggered. Stationary index (both high late) . If, in the late window, both PUMA and p21 are high (≥ µ PU , µ 21 ), we reuse a balance index based on late averages: We classify using the same two global cutoffs: Cycling is excluded since both effectors are above activation thresholds . Effective decision (severity-first) . Finally, we combine the early and stationary decisions. If both are defined, we report the most severe outcome: Return the effective fate If only one decision is available, we report that one. Table 6 summarizes our fate criterion, a more formal version of the criterion is given in the supplement Box 1. To simulate fate decision for a generic cell, the thresholds are set from experimental dynamics and modelling practice. Typically between 50–80% of simulated maxima Geva-Zatorsky et al. (2006) ; Batchelor et al. (2008) ; Becker et al. (2010) ; Barr et al. (2017) . Early thresholds are placed around 50–60%, final senescence higher (70–80%), then fine-tuned on simulated profiles for robustness and biological plausibility Ciliberto et al. (2005) ; Purvis et al. (2012) . Throughout the criterion, J ( τ ) denotes a sliding time window of duration τ over which activation or persistence is evaluated; for instance, means That p 21 remains above its threshold for at least consecutive hours. When both PUMA and p21 are simultaneously active, we define J * as the maximal overlapping interval where their joint activation conditions hold, that is, The balance indices I p and I f are then computed over J * and over the late window [60, 72] h, respectively. Global cutoffs α f (and analogously α p ) are fixed once by bounding late-window maxima with 95th percentiles across stresses, deriving an admissible interval [ I min , I max ], and choosing cutoffs by a small-margin heuristic or a minimax search. 3.1. Simulations for single cells In this section, we apply the threshold-based decision criterion illustrated in Figure 6 to ( generic ) single-cell simulations that reproduce the behaviours summarized in Table 1 . Besides constant stress and STAT1 levels, we also consider time-dependent inputs for S ( t ) and Σ( t ). The profiles used (with their biological meaning) are listed in Table 3 . View this table: View inline View popup Download powerpoint Table 3. Representative input profiles used to model stress S ( t ), active STAT1 (Σ( t )), or both, with corresponding mathematical expressions and biological/experimental interpretation. Download figure Open in new tab Figure 5. Paired single-cell simulations of the p53–Mdm2–p21–PUMA network under heterogeneous inputs. Each row corresponds to a distinct input scenario combining a fixed stress level ( S ) with either constant or time-varying STAT1 activity (Σ( t )). Odd columns display the dynamics of the core regulatory loop ( p 53, green; Mdm2, red), while the adjacent even columns show the corresponding downstream effectors ( p 21, blue; PUMA , orange). The diversity of oscillatory and transient behaviours illustrates how similar upstream conditions can produce distinct dynamical signatures, underscoring the need for a quantitative, threshold-based criterion to decode cell-fate outcomes. Download figure Open in new tab Figure 6. Decision rule with explicit gating: the Both high column is evaluated only if neither of the first two column conditions applies. Top: stationary means on [60, 72] h with I f = ⟨PUMA⟩ / (⟨ p 21⟩ + ⟨PUMA⟩). Bottom: early rule on sliding windows J ( τ ) with I p = PUMA / ( p 21 + PUMA). Example combinations for simulations. S ( t ) and Σ( t ) are treated as independent inputs combined to emulate common experimental settings: Decreasing stress + offset STAT1 peak: acute irradiation (fast DNA repair ⇒ S decays) followed by a delayed IFN-I/STAT1 activation through cGAS–STING signaling Brzostek-Racine et al. (2011) ; Yu et al. (2015); Feng et al. (2020) . Pulsed stress + pulsed STAT1: repeated or fractionated S ( t ) and bolus-like Σ( t ) (e.g. IFN dosing). Use case: fractionated radiotherapy or chemotherapy combined with IFN- α/γ treatment Heitmeier et al. (1999) ; Giovannozzi et al. (2020) ; Kalliara et al. (2022) . Slowly varying stress + tonic STAT1: slowly increasing or periodic S ( t ) (chronic stress) combined with a tonic Σ( t ) baseline (chronic inflammation) subject to slow attenuation. Use case: inflammatory microenvironments providing a persistent STAT1 background on which stress fluctuates Marié et al. (2021) ; Giles et al. (2017) . Box 1. Global — Decision thresholds (stress-invariant) Early PUMA (APO): Early p 21 bands: Stationary gates (mean on [60, 72] h): Index cutoffs (for I p and I f ): Cell fate under time-varying STAT1 To isolate the effect of STAT1, we fixed the stress level at a single value ( S = 17) and varied only the STAT1 signal Σ( t ). We considered four representative temporal profiles—exponential decay, rising ramp, offset (asymmetric) peak, and sinusoid—each evaluated at two maximum amplitudes (Σ max ∈ {50, 120}). All other parameters and initial conditions were kept identical. For each condition we simulated the p 53– Mdm 2 module together with the effectors p 21 and PUMA, then applied our fate criterion. When Σ max = 50 (rows A–D): In the p 53–Mdm2 module, oscillations are sustained but of moderate amplitude for the decreasing and sinusoidal inputs, while the increasing profile drives a slow damping toward steady-state levels—consistent with irreversible senescence. The system remains in a low-STAT1 regime where PUMA activation is limited. No early apoptotic commitment occurs; instead, p 21 dominates the response. Depending on the input waveform, the cell either stabilizes in a moderate-amplitude arrested state or gradually transitions toward a senescent plateau. When Σ max = 120 (rows E–H): Across these conditions, p 53–Mdm2 oscillations are more strongly damped and settle around higher mean values, reflecting enhanced stabilization of nuclear p53. Elevated STAT1 strengthens PUMA induction, pushing trajectories toward apoptosis. In all cases, PUMA surpasses its apoptotic threshold early, although the timing and sharpness of commitment depend on waveform shape. The decreasing and asymmetric-peak profiles trigger the fastest apoptotic responses ( t ∼ 20–25 h), whereas the sinusoidal input produces alternating windows of susceptibility: in-phase peaks with stress favor apoptosis, while out-of-phase segments allow transient arrest or delayed senescence. Overall, these comparisons highlight that the temporal architecture of the STAT1 signal—its rise and decay kinetics, duration, and frequency—can substantially reshape fate outcomes even under the same stress intensity. By synchronizing (or desynchronizing) STAT1 peaks with p 53 pulses, cells may shift between arrest, senescence, and apoptosis despite identical stress inputs. Cell fate under time-varying STAT1 and stress We vary both inputs in time and present the simulations in two complementary grids, keeping the paired layout (left: p 53–Mdm2 with scaled S ( t ) and Σ( t ); right: p 21–PUMA with the early band and the stationary window [60, 72] h). Below is a detailed reading of Figure 8 Download figure Open in new tab Figure 7. Impact of STAT1 waveform shape and amplitude on the p53–Mdm2–p21–PUMA dynamics at fixed stress ( S = 17). Each row corresponds to a distinct STAT1 input profile (decreasing, increasing, shifted peak, sinusoidal), while paired columns show the p 53–Mdm2 module (left) and downstream effectors p 21 and PUMA (right). The upper half (Σ max = 50) represents weak STAT1 activation, leading to arrest or senescence dominated by p 21. The lower half (Σ max = 120) illustrates strong STAT1 signaling, where PUMA surpasses its apoptotic threshold and early commitment to apoptosis occurs. Colored shaded bands mark the decision windows (early vs. stationary). Together, these simulations demonstrate that not only the amplitude but also the temporal shape of the STAT1 signal—its rise, decay, and oscillatory phase—critically modulates fate determination under the same stress intensity. Download figure Open in new tab Figure 8. Coupled profiles S ( t ) and Σ( t ). Row A — Decaying S ( t ). (A1) Σ( t ) has an offset peak while S ( t ) decays, yielding initially large but damped p 53 pulses. (A2) Early APO appears when the first overlap is strong, but PUMA later falls and the [60, 72] h average stays sub-plateau, giving CY at stationarity (APO → CY by lack of duration). Row B — Decaying S ( t ) with in-phase sinusoidal Σ( t ). (B1) Sinusoid crests arrive as S ( t ) is already declining, fragmenting the duty-cycle of p 53. (B2) Early ARR from p 21 and ARR (or CY) at stationarity; phase cannot compensate for the loss of duration caused by decaying S ( t ). Row C — Increasing S ( t ).. (C1) Pulses of p 53 lengthen gradually as S ( t ) rises. (C2) The single Σ peak occurs too early, mainly feeding p 21: early ARR, and ARR (or CY) at stationarity because PUMA is insufficient within [60, 72] h. Row D — Increasing S ( t ) with in-phase sinusoidal Σ( t ). (D1) A well-timed Σ crest on high S ( t ) produces a long p 53 pulse near the stationary window. (D2) Early SEN (or ARR depending on thresholds) followed by APO final; constructive overlap in [60, 72] h carries PUMA across the plateau (SEN → APO). Row E — Double-pulse S ( t ).. (E1) Two stress bursts; the Σ peak aligns mostly with the first. (E2) No early band and CY final: the second pulse lacks STAT1 assistance, so the useful PUMA area remains too small. Row F — Double-pulse S ( t ) with in-phase sinusoidal Σ( t ). (F1) Several partial overlaps; some p 53 peaks are tall but short-lived. (F2) CY final despite decent PUMA maxima: the [60, 72] h mean stays below µ PU . Duration beats instantaneous amplitude. Row G — Staircase S ( t ). (G1) Each plateau lengthens p 53 pulses; the Σ peak is early and then decays. (G2) No early band and CY final: staircase growth without a late Σ assist is insufficient for PUMA in the stationary window. Row H — Staircase S ( t ) with in-phase sinusoidal Σ( t ). (H1) A Σ crest on a high plateau yields a very long p 53 pulse just before stationarity. (H2) Early ARR from p 21, then APO final: PUMA crosses the plateau within [60, 72] h (ARR → APO). Row I — Pulse train S ( t ). (I1) Many narrow p 53 pulses; the single Σ peak covers only a fraction. (I2) No early band and CY final: the useful integral is too small despite repeated spikes. Row J — Pulse train S ( t ) with in-phase sinusoidal Σ( t ). (J1) Repeated short overlaps filtered by the p 53–Mdm2 loop. (J2) No early band but APO final: a sufficiently long late overlap occurs inside [60, 72] h and lifts PUMA above threshold. Row K — Single Gaussian S ( t ). (K1) Very strong but brief stress; Σ then decays. (K2) No early band and CY final: large PUMA peak but insufficient duration for the stationary mean. Row L — Single Gaussian S ( t ) with in-phase sinusoidal Σ( t ). (L1) Alignment can extend one p 53 pulse, yet not long enough within [60, 72] h. (L2) No early band and CY final: confirms that, for an isolated pulse, the mean over the stationary window is decisive. The common motifs throughout the grid are: (i) Duration inside the stationary window dominates, not the instantaneous maxima. (ii) Phase matters: constructive overlap near [60, 72] h drives APO even when no early band is present. (iii) Asymmetry between PUMA and p 21: p 21 accumulates and often yields early ARR/SEN, whereas PUMA requires long, well-timed p 53 pulses, hence frequent early–final mismatches. (iv) The p 53–Mdm2 loop behaves as a low-pass filter; fast or sparse overlaps are attenuated. All full time courses for time varying stress or Stat1 profiles and different initial conditions for p 21 and PUMA are available in Supplementary Section S4 on the Figure S5 , Figure S6 and Figure S7 . 4. From single-cell rules to population outcomes (via archetypes) 4.1. From single cells to archetypes Our goal is to test whether the threshold–based fate criterion, validated on a generic single cell, also explains outcomes at the population level—even when different populations display lineage–specific behaviours. To do so, we apply the single–cell decision rule of Fig. 6 at population scale: each cell is classified independently with the same severity–first priority CY ≺ ARR ≺ SEN ≺ APO. Population fractions are simple averages across cells, and uncertainty is reported by nonparametric bootstrap. We represent a cell population by an archetype 𝒜 = (Θ, 𝒟 S , 𝒟 Σ , ℋ). Here, Θ denotes the ODE parameter set; 𝒟 S and 𝒟 Σ specify the stress and STAT1 inputs (constants or sampled distributions, e.g. log–normal micro–stress or STAT1 dispersion); and ℋ collects the decision thresholds used by the rule in Fig. 6 (early and stationary cut–offs). Differences between populations are therefore encoded by changes in Θ, in the input statistics (𝒟 S , 𝒟 Σ ), and—when distinct sensitivities are required—by archetype–specific thresholds ℋ; the dynamical equations are not modified. Single–cell heterogeneity within an archetype is introduced by drawing, for each cell, a parameter vector from independent log–normal laws centred on the archetype “Center” values in Table 5 with the listed coefficients of variation (CV). For any parameter p with mean µ p and CV c p , View this table: View inline View popup Download powerpoint Table 4. Representative population outcomes ∼48–72 h after single-dose IR. Percentages are approximate unless stated. View this table: View inline View popup Download powerpoint Table 5. Randomly sampled parameters: central value and CV, by archetype. truncated to physiological bounds (±3 s p ). CV tiers reflect expected single–cell dispersion while maintaining oscillator robustness: 8–10% for affinity/threshold constants ( K ’s), 15% for decay/catalytic rates (12% in U2OS), and 20–25% for activation gains and basal terms. STAT1 activity Σ is drawn log–normally around an archetype–specific mean with CV 20–25%; a micro–stress multiplier on S uses a small CV (5–6%) to mimic dose dispersion; initial conditions for p 21 0 and PUMA 0 are log–normal (mean 0.5, CV 20%), shifted and truncated at 0. Draws are independent; when comparing doses for the same cell, the sampled vector is held fixed to avoid Monte Carlo re–labelling. Archetype–specific decision thresholds ℋ are tuned once at S = 10 and then kept fixed for hold–out S . Pre–registered constraints on the p53–Mdm2 core To allow modest lineage–specific kinetics without over–parameterisation, each archetype permits only two bounded deviations in the p53–Mdm2 core: (i) the Mdm2 protein degradation rate δ M , and (ii) the transcription/translation delay ratio k tr /k tl . Both are restricted to ±10% around reference values, with the product k tr k tl held constant to preserve the ∼ 6 h oscillation period. All other core constants ( K 0 , k u , K u , δ p , …) are fixed. Calibration and population readouts We calibrated three archetypal lines (HCT116–like, MCF–7–like, U2OS–like) to match representative population outcomes ∼48–72 h after single–dose IR ( Bae et al. (2018) ; Bravatà et al. (2015) ). Calibration adjusts only lineage–specific central parameters and, when needed, the STAT1 co–activation strength and the threshold set ℋ; the zoned decision rule and all ODEs remain unchanged. For each archetype, fate fractions f X (𝒜) are computed as averages of single–cell labels, and uncertainty is quantified by nonparametric bootstrap. 4.2. Population simulations We simulate N =1500 cells per archetype and restrict S to literature-supported levels. We present here the results for the MCF-7–like archetype as introduced in Section S5. The same procedure is applied unchanged to the HCT116–like and U2OS–like archetypes; their complete single–cell panels, population scatters, fate fractions, and sensitivity checks are provided in Supplementary Figures Figure S9 – Figure S11 and Supplementary Boxes 2-3. A compact comparison across all three archetypes is shown in Figure 11 . The archetype analyses below reuse exactly the single-cell gates defined above; for each archetype, cutoffs and I f ( I p ) partitions are held stress-invariant (no per- S threshold tuning). Download figure Open in new tab Figure 9. MCF-7 single-cell trajectories at Σ = 80 for three stress levels. Shaded red bands mark the first early-commitment window; the grey band is the stationary window (60–72 h). (A) S = 3.5: both effectors remain low; no early gate is crossed and the stationary call is cycling (CY). (B) S = 10: p 21 pulses build and cross the senescence band while PUMA stays sub-threshold; the cell is classified SEN both early and late. This reflects the p 21 branch dominating at moderate stress. (C) S = 24: PUMA rises to a high plateau and p 21 settles at a mid/high level—its production attenuates at high stress in this MCF-7 parametrization—yielding an early SEN window and a SEN decision in the stationary phase. Download figure Open in new tab Figure 10. MCF-7-like simulation results for population). Population-level simulations for the MCF-7 archetype, showing the distribution of fates across cells. Download figure Open in new tab Figure 11. Summary of the 3 archetype simulations 4.3. MCF–7 archetype We expect a small amount of arrest and apoptosis at very low stress, but above all massive senescence instead of apoptosis at extreme stress. The specific parameter values used for the MCF-7 archetype are listed in Table 6 : View this table: View inline View popup Download powerpoint Table 6. MCF-7-specific parameters used for the population simulations ( N =1500 cells per stress level). Σ denotes the baseline STAT1 activation level (log-normal dispersion CV = 0.25), while S is the DNA-damage stress input scanned in [0, 25]. All other constants are inherited from Table 5 . Decision thresholds and α indices are given in the box bellow. Box 2. MCF -7-like — Decision thresholds (global) Early PUMA (APO): Early p 21 bands: Stationary gates (mean on [60, 72] h): µ PU = 1.2, µ 21 = 0.013, Index cutoffs (used for I p and I f ): To reproduce the experimental MCF–7 sequence (cycling → arrest → senescence) without altering the ODE structure or the global thresholds, we used a physiological zone calibration . In this scheme, only context variables that biologically depend on the stress regime were adjusted: mean STAT1 activation (Σ), branch weights ( φ pu , φ 21 ), and the fraction of cells pre-activated in p21 ( seeding fraction). Specifically, Σ was slightly increased at intermediate and high stress ( S =10, 24) to mimic the known enhancement of JAK/STAT1 signaling under DNA damage; the PUMA branch was mildly attenuated while the p21 branch was reinforced at high stress, reflecting the experimentally observed bias toward senescence in MCF–7 cells. These modulations affect only the statistical initialization of the cell population and not the dynamical equations or the decision thresholds ( α 1,2 , θ, µ ), thereby preventing overfitting. The resulting configuration remains fully predictive within each stress zone and captures the physiological transition from mixed arrest/senescence at S =10 to irreversible senescence at S =24. For each stress level S and fate k ∈ {CY, ARR, SEN, APO}, the literature specifies an admissible interval [ ℓ S,k , u S,k ]. From N =1500 simulated cells we estimate , the empirical fraction per fate. On the right-hand panels, colored bars represent the literature intervals, and the black ticks show the simulated ; a green or red mark indicates whether . MCF–7-like: why at high S we have SEN (threshold view) We keep the ODEs unchanged and modulate the outcome solely through global thresholds . (i) Early APO remains rare because the PUMA gates are stringent: θ PUMA = (1.5, 11.5, 12.0) with durations τ PUMA = (6, 4.5, 3) h; short or modest PUMA transients fail to reach these levels. (ii) Late SEN is easy to attain: the stationary gates are µ PUMA = 1.2, µ 21 = 0.013, with a very low senescence trigger ; at high stress, p21 steadily surpasses this limit while PUMA remains below µ PUMA . (iii) When both markers coexist, the overlap index uses two global cutoffs α = (0.35, 0.73), so the typical p21 ≫ PUMA ratio is classified as SEN, with only clearly PUMA-dominant overlaps promoted to APO. Net effect: senescence becomes the dominant fate at large S , consistently reproducing the MCF–7 high-stress phenotype. We assess, for each stress S 0, 3.5, 10, 24 and fate k ∈ {CY, ARR, SEN, APO}, the point estimate against the admissible literature interval [ ℓ S,k , u S,k ], using the coverage indicator c S,k , the normalized gap g S,k , and the binomial Wilson 95% confidence interval [lo S,k , hi S,k ] as defined in Sec. S5.1. S =3.5. Empirical fractions fall within their bands CY ∈ [0.25, 0.45], ARR ∈ [0.20, 0.35], SEN ∈ [0.30, 0.50], APO ∈ [0, 0.10]. Hence c S,k = 1 and g S,k = 0 for all k , giving C 3.5 = 1.00 and G 3.5 = 0. Wilson 95% CIs ( N =1500) intersect the bands in all cases (e.g. for ; for ; for , the Wilson upper bound ≲ 0.01 remains inside [0, 0.10]). Decision: PASS . S =10. Empirical fractions ≈ (0.00, 0.23, 0.61, 0.16) for (CY, ARR, SEN, APO) lie inside the literature bands [0, 0.10], [0.20, 0.35], [0.50, 0.70], [0.10, 0.25], respectively. Thus c S,k =1 and g S,k =0 for all k , yielding C 10 =1.00 and G 10 =0. Wilson CIs overlap the bands in every class (e.g. ). Decision: PASS . S =24. Empirical fractions ≈ (0.00, 0.01, 0.93, 0.07) are within the target bands [0, 0.05], [0, 0.10], [0.70, 0.95], [0, 0.10]. Hence c S,k =1, g S,k =0, C 24 =1.00, G 24 =0. Wilson CIs remain inside/overlapping the bands (e.g. ). Decision: PASS . S =0. At basal stress, the model predicts a fully cycling population ( CY ≈ 1.00) with negligible activation of p53, p21, or PUMA. This outcome does not arise from missing literature bands but from the intrinsic structure of the model: at S =0, the transcriptional drives for p53 and its downstream targets vanish, preventing protein production and hence any irreversible response. In agreement with Fig. S7 , when initial concentrations remain below the thresholds of irreversible fates (apoptosis or senescence), the system stays trapped in the cycling basin corresponding to its stress level. Accordingly, interval-based validation is not applicable here; we report and for the residual classes, acknowledging that the lack of activation reflects a physiological rather than statistical limitation of the model. Summary For MCF–7-like populations at S =3.5, 10, 24, we obtain full interval coverage with C S =1.00, G S =0, and Wilson 95% CIs overlapping the literature intervals in all classes, thus meeting the PASS criterion. At S =0, the available classes also pass , while CY is not band–validated due to missing literature range. These results confirm a consistent calibration across the MCF–7 stress sequence. 4.4. Population–level agreement for HCT116-like and U2OS–like HCT116-like Across the three stress points ( S =0, 6.5, 10), the agreement between simulated and literature fractions is uniformly strong ( grit: 12/12 ,✓). At S =0, the cohort remains fully cycling ( CY ≈ 1.00) with negligible other classes; all bins lie within their admissible ranges, Wilson 95% confidence intervals intersect the bands, and global distances are minimal. At S =6.5, the population split ( ARR ≈ 0.53, CY ≈ 0.45, SEN ≈ 0.01, APO ≈ 0.00 ) fits squarely inside the literature intervals, yielding a non-significant goodness-of-fit and small classwise residuals. At S =10, the response shifts toward apoptosis and senescence ( APO ≈ 0.55, SEN ≈ 0.24, ARR ≈ 0.12, CY ≈ 0.09 ); all values fall within the expected ranges, with tight CIs and low normalized gaps. Overall, HCT116-like simulations display excellent concordance from mid to high stress without systematic deviations. U2OS–like Over the three evaluated levels ( S =0, 8, 10), the overall agreement is high ( grid: 10/12 ✓). At S =0, the simulation yields a nearly pure cycling population ( CY ≈ 1.00), whereas the literature includes a minor arrest fraction; this explains one of the two misses ( ARR S =0 ) and reflects the absence of p53 activity at basal stress rather than a numerical error. At S =8, the simulated mix ( SEN ≈ 0.49, ARR ≈ 0.29, CY ≈ 0.19, APO ≈ 0.02 ) lies entirely within literature intervals with overlapping Wilson CIs and very small residuals, indicating near-perfect alignment. At S =10, the simulation slightly overestimates cycling ( CY ≈ 0.06 vs. band upper bound ≈ 0.05), while APO, SEN , and ARR remain inside their respective ranges. This deviation is marginal and does not alter the fate hierarchy ( SEN > APO > ARR > CY ). Global distances stay below 0.1 and all Wilson CIs intersect their bands. Across archetypes, mid-stress points again yield the tightest SIM–LIT alignment. The two isolated deviations ( ARR S =0 in U2OS, CY S =0 marginally above its band) are mechanistically interpretable: they reflect basal inactivity or residual proliferative potential rather than parameter mismatch. Elsewhere, every fate lies within its literature band with overlapping 95% Wilson confidence bounds and normalized gaps G S ≤ 0.1, confirming excellent global concordance under stress-invariant gates. Discussion Our stress–aware model couples the p53–Mdm2 core to p21 / PUMA and treats STAT1 (Σ) as a modulatory input. In doing so, it recovers hallmark features of the DNA-damage response while adding a second, cytokine-controlled axis. First, by sweeping S without re-tuning parameters, we reproduce the canonical sequence of dynamical regimes (self-sustained pulses, damped pulses, plateau), thereby confirming systems-level predictions and single-cell measurements on p53 oscillatory control ( Bar-Or et al., 2000 ; Lahav et al., 2004 ; Ciliberto et al., 2005 ; Proctor & Gray, 2008 ; Geva-Zatorsky et al., 2006 ; Batchelor et al., 2008 , 2011 ). Within the oscillatory domain, sustained p53 activity confirms the preferential engagement of p21 and arrest/senescence programs ( Harper et al., 1993 ; Xiong et al., 1993; Abbas & Dutta, 2009 ; Stewart-Ornstein & Lahav, 2017 ; Hanson et al., 2019 ). Outside this domain, the collapse of oscillations toward a steady state agrees with strong PUMA induction and an apoptotic bias ( Nakano & Vousden, 2001 ; Jeffers et al., 2003 ; Kracikova et al., 2013 ; Zhang et al., 2007 ). Second, our explicit ( S , Σ) bifurcation map complements prior descriptions of the p53–Mdm2 axis by adding a tunable STAT1 dimension that shifts the Hopf band and reshapes fate partitions. In particular, increasing Σ both weakens Mdm2 -mediated p53 removal and co-activates effector promoters, which aligns with reports that interferon/STAT signals feed back into p53 networks and IFN programs ( Levy & Darnell, 2002 ; Brzostek-Racine et al., 2011 ; Feng et al., 2020 ). Our timing experiments further show that sustained or promptly rising Σ( t ) accelerates early-gate crossing and increases transient PUMA load, while brief or delayed Σ( t ) yields milder, often reversible effects. This extends observations on stimulus-dependent decoding of p53 dynamics ( Batchelor et al., 2011 ; Purvis et al., 2012 ; Becker et al., 2010 ; Wang et al., 2023 ) by making the STAT1 time-profile an explicit gain/timing knob. At population scale, injecting realistic heterogeneity into S and Σ reproduces lineage-specific fate biases without re-fitting biochemical parameters: HCT116-like populations concentrate on p21 -dominated arrest/senescence, U2OS-like populations show enhanced PUMA responsiveness and apoptosis, and MCF-7-like populations remain mixed. This confirms and refines the idea that fate depends not only on p53 amplitude but also on target stability and feedback efficiency ( Hanson et al., 2019 ; Stewart-Ornstein & Lahav, 2017 ; Sun et al., 2009 ). It also complements literature on S-phase damage and quiescence decisions via p21 ( Barr et al., 2017 ), by embedding these outcomes in a two-input ( S , Σ) landscape that predicts cohort-level splits. Regarding STAT1, our results converge with evidence that prolonged STAT1 activity can potentiate cytotoxic and apoptotic programs ( Heitmeier et al., 1999 ; Giles et al., 2017 ; Kalliara et al., 2022 ; Townsend et al., 2004 ) and are consistent with STAT1–p53 cooperativity on damage-induced apoptosis. At the same time, the literature also documents context-dependent, sometimes pro-tumoral STAT1 roles; our model adds that the temporal form of Σ( t ) (sustained versus pulsed) is a key discriminator that may reconcile seemingly opposing observations ( Meissl et al., 2017 ). Mechanistically, placing the STAT1 dimension alongside S clarifies how upstream cytokine cues reshape intracellular decoding by the p53 network and its effectors. Finally, our framework suggests two natural extensions. A compartmental version will separate nuclear and cytosolic pools and translocation feedbacks, building on established transport models ( Dimitrio, 2012 ; Eliaš et al., 2014 ). Coupling this stress–response module to a mechanistic JAK–STAT pathway model would close the loop from extracellular ligands to intracellular Σ( t ) and fate decoding ( Levy & Darnell, 2002 ), enabling predictions testable by combining controlled interferon waveforms with DNA-damage protocols and single-cell readouts of p53, p21, and PUMA. Supplementary Information S1. The model S1.1. Calibration of the stress scale S We sought a single, dimensionless stress variable S ∈ [0, 25] that captures how UV and γ irradiation stabilize p 53 by weakening its Mdm2–mediated ubiquitination. Table S1 compiles empirical anchors from the literature: for representative UV or γ doses we list the measured p 53 half–life t 1 / 2 , the corresponding effective ubiquitination rate k u ≃ ln 2 /t 1 / 2 , and the S value assigned to that condition. These anchors cover the full range from basal conditions ( S =0) to severe damage ( S ∼25). To complement these kinetic anchors, Table 1 in the main article summarizes the phenotypic regimes associated with increasing S : qualitative biological effects and the characteristic p 53 dynamics observed near the corresponding γ -dose benchmarks. UV exposures were placed on the same S scale by matching p 53 stabilization features (half–life/amplitude) to the nearest γ -dose regime. View this table: View inline View popup Download powerpoint Table S1. Mapping between dose, p 53 half–life ( t 1 / 2 ), inferred ubiquitination rate ( k u ), and the assigned stress level S . Calibration procedure (summary) In the model, stress modulates the effective ubiquitination activity through the factor so that higher S reduces p 53 degradation. We calibrated ( K 0 , n ) as follows: Use Table S1 to anchor S values to measured t 1 / 2 (and thus k u ) under UV/ γ . For each anchor and for a grid of ( K 0 , n ), simulate the model and extract summary features of p 53 dynamics (amplitude, period, damping) and stabilization. Minimize the aggregate mismatch between simulated features and the regimes in Table 1 (e.g., oscillatory vs. plateau, strength of stabilization) while respecting the kinetic constraints from Table S1 . The best compromise was K 0 = 12 and n = 6, which reproduces the transition from basal to pulsatile to stabilized p 53 across the S range. This procedure operationalizes the stress scale: the first table fixes kinetic anchors (half–life/ubiquitination), the second enforces dynamical phenotypes, and the fitted ( K 0 , n ) provide a consistent mapping from dose to S within the model. S2. Stability and bifurcation S2.1. Existence, Uniqueness, and Positivity of Solutions Existence, positivity, and uniqueness of solutions are given by the following theorem: Existence, uniqueness, and positivity of solutions Consider the system of differential equations defined on : For any initial condition x 0 ∈ Λ and any pair of continuous functions ( S ( t ), Σ( t )), there exists a unique bounded solution x ( t ) defined on [0, +∞[, which remains in Λ for all t ≥ 0, i.e., all concentrations stay nonnegative at all times. Proof . Existence of solutions All building blocks of F in x are 𝒞 ∞ on Λ (polynomials, products, and quotients whose denominators are of the form p 53 + K u , S 4 + K 4 , p 53 4 + K 4 , etc., with strictly positive constants). Hence, for each t ≥ 0 and continuous S , Σ, the map x → F ( t , x ) is 𝒞 1 on Λ, thus locally Lipschitz. By the Picard–Lindelöf theorem, there is a unique local solution from any x 0 ∈ Λ. Positivity We examine the vector field on the boundary of the non-negative orthant At the boundary p 53 = 0, the p 53 equation gives At the boundary M dm 2 RNA = 0, the M dm 2 RNA equation gives At the boundary M dm 2 = 0, the M dm 2 equation gives At the boundary p 21 = 0, the p 21 equation gives At the boundary PUMA = 0, the PUMA equation gives with g pu ( S ) ≥ 0. In each case the vector field points inward (or is tangent) on the coordinate hyperplanes. Therefore, Λ is a positively invariant set for the system. Boundaries for the core module. On p 53 = 0 On M dm 2 RNA = 0, On M dm 2 = 0, Boundaries for the effectors . On p 21 = 0, where g 21 ( S ) ≥ 0 is the non–negative Hill term used in the model. On PUMA = 0, with g pu ( S ) ≥ 0 the stress–dependent Hill factor used for PUMA . Since each boundary derivative is non–negative, trajectories starting in Λ remain in Λ for all forward times; i.e., Λ is a positive invariant set for the system. Boundedness We first recall (differential) Grönwall’s inequality, which we will use repeatedly. The Grönwall–Bellman inequality (differential form) . Let u ( t ) be continuously differentiable on [0, T ] and satisfy Then, for all t ∈ [0, T ], If a = 0, then u ( t ) ≤ u (0) + b t . Below we use that every factor of the form lies in [0,1], and that . p 53( t ) is bounded . By Grönwall, M dm 2 RNA ( t ) is bounded . Using p 53( t ) ≤ P max gives Thus, M dm 2( t ) is bounded . Therefore, p 21( t ) is bounded . since 0 ≤ g 21 ( S ) ≤ 1 and all multipliers are ≤ 1. Thus, PUMA ( t ) is bounded . Hence, Combining (i)–(v), every component of the solution remains uniformly bounded for t ≥ 0. Together with positivity (shown above), this implies that trajectories remain in a compact, positively invariant subset of the nonnegative orthant. Moreover, since the right-hand side is 𝒞 1 in the state variables on Λ (hence locally Lipschitz), the Picard–Lindelöf theorem ensures existence and uniqueness of a local solution; the a priori bounds extend it to a unique global solution for any initial condition in Λ. This establishes that the problem is well posed both mathematically and biologically. S2.2. Global sensitivity analysis (LHS–SRC, final–mean outputs) Method We quantified how model parameters shape the dynamics of the p 53–Mdm2 core and of the effectors p 21 and PUMA using Standardized Regression Coefficients (SRC) on samples generated by Latin Hypercube Sampling (LHS). S2.3. Bifurcation loci in the (S, Σ) plane Fig. S2 shows the real part of the leading eigenvalue as S is varied for several fixed Σ. A positive real part indicates an unstable fixed point and, here, the onset of a Hopf oscillatory regime. For Σ = 0 (blue) the Hopf window spans S ≈ 3.6–22.7. At Σ = 100 (red) it shrinks to S ≈ 4.3–15.1, and at Σ = 250 (green) to S ≈ 5.9–12.1. For Σ = 500 (orange) the curve remains negative for all S explored, hence no Hopf crossing. Two systematic trends emerge: (i) the onset (dotted lines) moves to the right as Σ increases; (ii) the termination (dash-dotted lines) moves left. Consequently, the oscillatory band narrows with Σ. Download figure Open in new tab Figure S1. Global sensitivity by final-mean outputs (SRC). LHS (±30%) and SRC on ( left ) p53 metrics (amplitude, #peaks, period) and ( right ) final-mean p21 and PUMA (60–72 h). Positive bars indicate parameters that increase the output when raised; negative bars indicate the opposite. Results reveal partly distinct regulators for p53 oscillatory behaviour versus effector accumulation. Download figure Open in new tab Figure S2. Influence of Σ on the Hopf threshold S * in the p21 model. Real part of the dominant eigenpair at while scanning S for multiple Σ. Vertical dotted lines mark ; vertical dash-dotted lines mark . View this table: View inline View popup Download powerpoint Table S2. Most influential parameters from the SRC analysis (final–mean outputs). Arrows indicate the direction of the standardized regression coefficients (SRC). Only the dominant effects (largest |SRC| in Fig. S1 ) are listed. These trends can be understood as follow: STAT1 enters the core loop through the factor in the M dm 2 equation; increasing Σ lowers this factor, reducing Mdm 2 production and thus the effective p 53 → Mdm 2 negative-feedback gain. Sustained oscillations require enough loop gain and delay; when Σ is high the fixed point is more strongly damped, so a larger stress is needed to destabilize it (rightward shift of ). At the same time, in our model S decreases k cat (the M dm 2-mediated p53 removal rate), which weakens the feedback at high stress and brings the system back to stability—explaining the leftward shift of and the progressive collapse of the oscillatory window. Interpretation Left: p 53 amplitude (60–72 h) forms a Hopf “tongue” at intermediate S . As Σ increases, this oscillatory band shifts to higher S and narrows , ultimately extinguishing oscillations at large Σ. Middle: p 21 amplitude (60–72 h) mirrors the p 53 tongue but is narrower and decays faster with Σ; when Σ stabilizes p 53 on a high plateau, the mean of p 21 can remain moderate/high while its amplitude is low. Right: max(PUMA) (0–100 h) increases monotonically with S and Σ, reflecting stress sensitization and STAT1 co-activation. Fate reading: the oscillatory band aligns with arrest (ARR); its progressive closure with growing Σ promotes damped/intermediate responses (consistent with SEN) and, at high S and/or Σ, large PUMA peaks bias decisions toward apoptosis (APO). Interpretation Figure S4 shows that, in addition to the stress input S , both the Mdm2 degradation rate δ M and the saturation constant K 0 act as bifurcation parameters for the p53 network. For small values of δ M , oscillatory regimes can emerge, and their extent depends moderately on the level of STAT1 (Σ). Similarly, low K 0 values may also trigger a Hopf bifurcation, but again the influence of Σ is weaker than in the case of stress S . Overall, these results indicate that while K 0 and δ M shape the boundaries of oscillatory domains, STAT1 mainly fine-tunes them, and the dominant bifurcation control remains associated with the stress level S . Download figure Open in new tab Figure S3. Hopf bifurcation for K 0 and δ M . Download figure Open in new tab Figure S4. Hopf bifurcation for K 0 and δ M . S3. Fate criterion, Calibration details and evaluation metrics Calibration of stress-agnostic stationary-index cutoffs α f We fix a single pair α f = ( α f 1 , α f 2 ), shared by all stress levels S . Stationary window and low thresholds Define the stationary window 𝒲 = [60, 72] h and the time averages We use three global thresholds: µ 21 (low p21), µ PU (low PUMA), and µ S > µ 21 (senescence-level p21). The stationary index I f is evaluated only when both species exceed their low thresholds. Global upper bounds for stationary concentrations To remain stress-agnostic, we estimate 95th-percentile upper envelopes across stresses: Admissible interval for the stationary index When both species maintain high stationary levels, we define the balance index which measures the relative contribution of pro-apoptotic PUMA compared to the total (PUMA + p21) pool. Let us denote A = ⟨PUMA⟩ and B = ⟨p21⟩. Under the admissible ranges we seek the interval of all possible values of I f ( A, B ) = A/ ( A + B ). Since the function I f ( A, B ) is strictly increasing in A and strictly decreasing in B . Therefore, the extremal values of I f over the admissible rectangle [ µ PU , U PU ] × [ µ 21 , U 21 ] occur at opposite corners: Hence, By construction, corresponds to the most p21-dominant configuration (low PUMA, high p21), while corresponds to the most PUMA-dominant configuration (high PUMA, low p21). The admissible interval thus bounds all biologically plausible stationary balances between pro-apoptotic and cytostatic dominance, and provides a stress-agnostic reference for the calibration of the global cutoffs α f 1 , α f 2 . (D) Choosing the two cutoffs α f 1 , α f 2 Two options are considered: Heuristic (fast) . With a safety margin δ ∈ [0.02, 0.05], This choice partitions the admissible interval into three regions: I f < α f 1 , α f 1 ≤ I f < α f 2 , and I f ≥ α f 2 , which correspond respectively to cycle/viability, arrest, and apoptotic regimes. Minimax (robust) . Perform a grid search over with α f 1 < α f 2 , minimizing the worst per-stress classification error: where Err S = 1 − F1 or 1 − accuracy on the population at stress S . (E)Evaluation metrics (precision, recall, F1) For each class (one-vs-rest), For multi-class fate (CY/ARR/APO), we report macro- and weighted-F1 scores; micro-F1 equals overall accuracy. (F) Mapping from I f to fate (explicit rule) With fixed global cutoffs 0 < α f 1 < α f 2 < 1, Note: senescence is now encompassed within the arrest band when p 21 exceeds µ S but I f < α f 2 , reflecting cytostatic rather than apoptotic dominance. (G) Validation (per stress) . For each stress level, compute per-class F1 scores; accept the global pair α f if all strata exceed a floor (e.g. F1 ≥ 0.70) and the across-stress variation remains moderate (range ≤ 0.10). Leave-one-stress-out validation avoids leakage. Early index cutoffs α p (analogous procedure) . For the early overlap window J ⋆ , and the same two-cutoff strategy ( α p 1 , α p 2 ) is applied to delineate early-cycle vs early-apoptotic dominance. Validation mirrors the stationary case. Our cell fate criterion reads: Supplementary Box 1 Cell-fate decision criterion (two global cutoffs) Setup . Fix a single set of thresholds, independent of stress S : Let J ( τ ) be any contiguous window of length ≥ τ . (1) Early decision (before 60 h) . (2) Early index (PUMA–p21 overlap) . If there exists with and for all t ∈ J , define Then classify with global cutoffs 0 < α 1 < α 2 < 1: (3) Stationary decision (60–72 h) . With ⟨·⟩ = ⟨·⟩ [60,72] : (4) Stationary index (both high) . If (⟨p21⟩ ⩾ µ 21 ) ∨ (⟨p21⟩ ⩾ µ S ) and ⟨PUMA⟩ ⩾ µ PU , define Then classify with the same global cutoffs: (5) Effective decision (severity-first) . Let Early ∈ {None, ARR, SEN, APO} and Stationary ∈ {CY, ARR, SEN, APO}. Return the effective fate S4. Simulations Here are simulations made on generic single cells in addition to the ones in the main article. The following table gives the threshold values used in the single generic cell simulations. View this table: View inline View popup Download powerpoint Table S3. Thresholds and durations used in the cell-fate decision criterion, with both the global stationary-index cutoffs ( α f 1 , α f 2 ) = (0.24, 0.50) and the early-index cutoffs ( α f 1 , α f 2 ) = (0.22, 0.50) obtained via the percentile-bounded calibration described in the text. Throughout, gates are stress-invariant: all cutoffs and early windows are fixed per archetype and reused for every stress level S , without per- S retuning. S4.1. Cell fate under constant stress and constant STAT1 We here consider constant levels of both stress and Stat1 Download figure Open in new tab Figure S5. Time courses of p 21 and PUMA for different fixed STAT1 levels Σ. Horizontal lines: thresholds; grey band: 60–72 h; tags indicate the effective fate. Interpretation . With constant inputs S and Σ, trajectories relax to a fixed attractor—no late reversals occur. When a coloured vertical band appears on the PUMA/ p 21 panel, it marks the time at which the early rule is first satisfied; under constant inputs this early label always matches the stationary label (60–72 h). If no band appears, fate is read from the stationary window only. At fixed S , raising Σ (STAT1) accelerates PUMA/ p 21 production and weakens the p 53–Mdm2 negative feedback, so commitment (when present) tends to occur earlier and oscillations are more strongly damped; in some cases the attractor class (oscillatory vs. plateau) changes and the fate switches. S = 4 (low stress): The p 53–Mdm2 loop decays to a low fixed point. Both p 21 and PUMA remain well below threshold; no early band is detected. The cell keeps cycling for Σ = 30 and for Σ = 100 (A2,B2: CY). S = 11 (moderate stress): The core shows a stable pulse train (limit cycle). For Σ = 30, p 21 crosses the arrest gate while PUMA stays sub-apoptotic— arrest without an early band (C2: ARR). For Σ = 100, the same outcome is reached earlier and is flagged by an early green band (D2: ARR→ARR). S = 14 (higher stress): With Σ = 30 the system remains pulsatile and p 21 settles in the arrest range— arrest with no early band (E2: ARR). With Σ = 100 the stronger STAT1 drive damps the oscillations to an intermediate plateau and maintains higher p 21, producing senescence with an early red band (F2: SEN→SEN). S = 24 (extreme stress): Trajectories collapse rapidly to a high plateau with strong PUMA and relatively low p 21. Both conditions end in apoptosis ; however, only the higher STAT1 produces an early apoptotic band (G2: APO without early band; H2: APO with an early purple band). Overall, constant-input experiments confirm that fate is governed by the steady attractor reached under ( S , Σ), while Σ primarily modulates the timing of commitment and, when it reshapes the attractor (oscillations → plateau), can shift ARR↔SEN. S4.2. Cell fate under time-varying stress We here assume stress is a time-dependant function exhibiting different profiles. Download figure Open in new tab Figure S6. Time courses for four stress profiles S ( t ) (decreasing, increasing, peak, sinusoidal) at S max =15 and 24. Interpretation (time–varying stress, Σ = 100) Coloured vertical bands mark the first time the early rule is satisfied; when no band appears, fate is decided only from the stationary window (60–72 h). We report fates as continued cycling (CY), cell–cycle arrest (ARR), senescence (SEN), and apoptosis (APO). Decreasing exponential A rapidly decaying input is too brief to sustain effectors. At S max = 15 there is no early band and the stationary decision is CY (A2). At S max = 24 a short p 21 overshoot triggers an early ARR band, but the signal vanishes before 60–72 h, so the stationary decision is CY (B2). Increasing ramp A steadily rising input prolongs activation. At S max = 15, p 21 crosses the arrest gate and the trajectory is ARR both early and at stationarity (C2). At S max = 24, the ramp drives PUMA across its threshold; APO is reached early and remains APO in the stationary window (D2). Shifted peak (narrow pulse) A short stress pulse produces isolated spikes in the p 53–Mdm2 loop that relax quickly. At both S max = 15 and S max = 24 there is no early band and the stationary decision is CY (E2, F2). Sinusoidal waveform Repeated input windows maintain p 53 /p 21 pulsatility. At S max = 15 there is no early commitment; the stationary window classifies ARR (G2). At S max = 24 commitment often requires accumulation over several cycles—there may be no early band—yet PUMA integrates to a high plateau and the stationary decision is APO (H2). Overall, waveform shape is decisive: sustained or rising inputs promote commitment (ARR at moderate intensity, APO at high intensity), brief/decaying inputs fail to meet stationary thresholds, and periodic inputs gate commitment in discrete windows, which explains why early and stationary calls can differ under the same peak stress. S4.3. Sensitivity to initial p21 and PUMA We wanted to see whether initial conditions do have any impact on cell fate. Download figure Open in new tab Figure S7. Effect of initial conditions for p 21(0) and PUMA(0) on fate. Interpretation The grid varies the stress level across rows ( S = 4, 11, 14, 24; Σ = 100 fixed) and the initial loads of p 21/PUMA across columns (the last column is a deliberately extreme “preload”). Coloured vertical bands indicate the first time an early gate is crossed; when absent, fate is decided from the stationary window (60–72 h). For physiological initial conditions (first three columns), the effective fate is essentially set by S ; initial conditions mainly shift the timing of commitment. S = 4. Trajectories relax to a low attractor and no early gate fires in the physiological columns; the stationary decision is CY. Only the extreme PUMA preload triggers an early APO band that overrides the stationary CY by the severity priority, illustrating the rule rather than genuine multistability. S = 11. The core remains on a sustained pulse train. p 21 crosses the arrest gate around 30–45 h (green bands) and the stationary window confirms ARR in all physiological columns; initial loads mostly advance or delay the gate but do not change the fate. S = 14. Pulses damp toward a mid–level plateau. In the physiological columns, p 21 crosses the SEN gate (red band, ∼15–25 h) and the stationary call is also SEN. Under the extreme preload, PUMA briefly exceeds its early apoptotic gate, yielding an APO early band, but the stationary window classifies SEN; the effective fate is APO by the priority rule. S = 24. Dynamics collapse to a high plateau with strong PUMA. An APO gate fires within 10–20 h in every column and the stationary decision is APO. Initial conditions only modulate the latency of commitment. Overall, constant Σ with varied S reveals convergence to a unique late attractor at each dose (low oscillations → CY; sustained limit cycle → ARR; mid plateau → SEN; high plateau → APO). Initial conditions alter when thresholds are crossed but not which fate is ultimately reached, except under unrealistically large preloads where the early priority applies. S5. From single-cell rules to population outcomes Our procedure to build a cell-line archetype follows a common protocol: Phenotypic targets (literature) . For each representative S , we fix the expected Cycle/ARR/SEN/APO trends (window ∼48–72 h) and, when available, early signatures. Decision logic . We adopt the criterion (early windows, stationary gates, and indices I p /I f ) as defined in the main article Box ?? ; the θ/τ and µ • thresholds, as well as the α cutoffs, are deterministic and global. Archetype constants/nudges . Certain amplitudes are adjusted fixedly to reflect lineage traits (e.g., MCF-7: PUMA dampening and slight p 21 boost; U2OS: mild PUMA dampening). These nudges define the central value around which cells are sampled. Intra-population variability . The parameters listed in the table below are drawn independently per cell from log-normal laws with the stated center and coefficient of variation (CV). We include a multiplicative noise factor on S (micro-variability), a draw for Σ (STAT1), and initial conditions IC ( p 21 0 , PUMA 0 ). Calibration . We tune first the stationary gates ( µ pu , µ 21 , µ S ), then the early windows ( θ, τ ), and finally the indices ( α p • , α f • ). Nudges are only revised as a last resort. S5.1. Archetypes Computation Given N single-cell simulations or trajectories c = {1, …, N }, we assign an effective fate Eff c ∈ {CY, ARR, SEN, APO} per cell using our fate criterion, then report fractions Uncertainty can be summarized by bootstrap CIs or by repeating with different random seeds. Supplementary Methods Validation against literature intervals. Setup For each stress level S , the simulation (SIM) outputs class proportions p = ( p CY , p ARR , p SEN , p APO ) Δ 3 ∈ (Σ i p i = 1). The literature (LIT) provides interval targets ℛ = { q ∈ Δ 3 : L i ≤ q i ≤ U i for i = 1, …, 4} with 0 ≤ L i ≤ U i ≤ 1 and Σ i L i ≤ 1 ≤ Σ i U i Feasibility and margin . Define the per-class violation Feasible if and only if v ( p ) ≡ 0. A normalized margin (unitless) highlights which class is tightest: Distance to the interval set (best achievable agreement) . We measure the closest agreement to any vector within the LIT ranges: with JSD the Jensen–Shannon distance (base 2). Both minimizations are over a convex polytope (Δ 3 intersected with the box [ L, U ]); they can be computed by a small quadratic program (for JSD, using its convex proxy) or by projecting p onto ℛ in ℓ 1 (for TV ⋆ ). Minimal reallocation mass . TV ⋆ has a direct interpretation: it is the smallest fraction of cells that must be re-labeled (moved between classes) to bring SIM inside all LIT intervals. We also report the edit plan (which classes donate/receive). Sample-size aware coverage . With n simulated cells, let be the observed proportions and draw B bootstrap replicates (multinomial). Define the interval coverage and per-class coverages . Coverage distinguishes tight agreement from incidental hits at small n . Class-wise standardized gap (sample-size free) . To summarize how far the worst class sits from its admissible band, Here g max = 1 means being one full band-width away from feasibility. Decision rule (single stress level). Acceptable if Needs improvement if TV ⋆ > 0.05 or g max > 0.50 but Cov [0.50, 0.80). Not acceptable otherwise. (Thresholds chosen for interpretability; they can be tuned in sensitivity analyses.) Aggregating across stresses Report the above metrics per S and summarize with medians (and IQR) over S . For a single composite score, use the maximum distance to feasibility: max S TV ⋆ ( p ( S ), ℛ ( S )), alongside the minimum coverage min S Cov( S ). Notes on implementation We solve TV ⋆ as a linear program on q with constraints L i ≤ q i ≤ U i , Σ i ≥ q i = 1, q i 0. JSD ⋆ is computed by a convex solver (or by grid search, given the low dimension). Bootstrap uses B = 200 by default; results are stable for B ≥ 100. S6. Simulations S6.1. HCT116 archetype Apoptosis and senescence already emerge at moderate stress. Download figure Open in new tab Figure S8. HCT116-like single–cell dynamics at two constant stresses (Σ = 100). Left: S = 6.5. p21 (blue) and PUMA (red) trajectories with the stationary window shaded (60™72 h). No stable early band is detected ( Early: None ); the late averages satisfy ⟨p21⟩ > µ 21 and ⟨PUMA⟩ < µ PU , yielding Final: SEN. Right: S = 10. A transient PUMA-dominant early band is detected ( Early: APO ), while the late window again satisfies the senescence gate ( Final: SEN ). The severity-first rule therefore reports Eff: APO . Insets list peak concentrations and the stationary index I f = ⟨PUMA⟩ / (⟨PUMA⟩ + ⟨p21⟩); the shaded rectangles mark the decision windows (purple for the detected early band, grey for the stationary window). Download figure Open in new tab Figure S9. Population-level simulations showing the distribution of fates for the HCT116 archetype. The set of thresholds used for the simulations are given in Box ?? . Supplementary Box 2 HCT116-like — Decision thresholds (global) HCT116-decision-critere Early PUMA (APO): Stationary gates (mean on [60, 72] h): Index cutoffs (used for I p and I f ): Early p 21 bands: Validation of population fractions against literature intervals (HCT116-like) For each stress level S and fate k ∈ {CY, ARR, SEN, APO} the literature provides an admissible interval [ ℓ S,k , u S,k ]. From N simulated cells we estimate the fraction and evaluate: Interval coverage and normalized gap . Aggregate per stress: Uncertainty-aware checks . For each ( S, k ) we compute a binomial Wilson 95% CI [lo S,k , hi S,k ] for and record Decision rule (per stress) . (Thresholds 0.15 / 0.25 are conservative and can be tuned.) Results (HCT116-like; N =1500 per S) Across the three stress levels ( S =0, 6.5, 10), the validation remains globally consistent and all cases are classified as PASS . At S =0, the population is entirely cycling ( CY = 1.00), which lies well within the literature range and confirms the absence of premature activation under basal stress. At S =6.5, the fractions ARR ≈ 0.53, SEN ≈ 0.01 , and CY ≈ 0.45 all fall within their admissible intervals, with APO remaining negligible—consistent with a dominant arrest response at intermediate stress. At S =10, the mixture shifts toward apoptotic and senescent outcomes ( APO ≈ 0.55, SEN ≈ 0.24, ARR ≈ 0.12, CY ≈ 0.09 ), still inside the expected bands, reproducing the experimentally observed stress-induced transition from arrest to apoptosis. In all panels, simulated ticks (black bars) fall within or touch the colored literature intervals, and every class shows a green checkmark—indicating full interval coverage ( C S =1.00), intersecting Wilson 95% confidence intervals, and small worst normalized gaps ( G S ≤ 0.10). S6.2. U2OS archetype For this archetype, we aim to observe senescence and apoptosis already at low stress. We use a standard U2OS-like single cell; the biochemical model below (parameters and initial conditions) is fixed across simulations. Download figure Open in new tab Figure S10. U2OS-like—p21 & PUMA ( S = 8 and S = 10). Gray bands mark the stationary window (60–72 h) and, when present, the colored(purple) band indicates a validated apoptosis early window. Left ( S = 8 ) . No early window detected. Right ( S = 10 ) . An early PUMA window is crossed. View this table: View inline View popup Download powerpoint Table S4. U2OS single-cell model parameters used in the simulations (fixed Σ = 110, S swept). Decision thresholds are listed elsewhere. The threshold values used for U2OS archetype are listed in the box bellow: Supplementary Box 3. U2OS-like — Decision thresholds (global) U2OS-decision-critere Early PUMA (APO): Stationary gates (mean on [60, 72] h): Index cutoffs (used for I f ): Early p 21 bands: For S = 8. No early window is detected ( Early: None ). The 60–72 h averages classify the cell as senescence ( Final: SEN ); the effective decision is SEN . There is no valid early overlap , hence I p is undefined. The stationary index is I f ≃ 0.492. For S = 10. A PUMA early window is crossed ( Early: APO ). Stationary averaging yields arrest ( Final: ARR ); the effective decision is therefore APO . Again there is no valid early overlap , so I p is undefined. The stationary index is I f ≃ 0.534. These checkpoints summarize what the figure shows (no I p at either stress; SEN at S =8, early APO at S =10) and will be used to introduce the updated U2OS parameter table in Table S4 . Using the decision rule i, we ran population simulations (N=2000 U2OS-like cells) at the literature-supported stress levels (S=0, 8, 10) and compared the simulated fate fractions with published benchmarks; the results are shown in Fig. S11 . Download figure Open in new tab Figure S11. Simulations for the U2OS archetype. Validation of population fractions against literature intervals (U2OS-like) For each stress level S and fate k ∈ {CY, ARR, SEN, APO}, the literature specifies an admissible interval [ ℓ S,k , u S,k ]. From N =1500 simulated cells we estimate , the empirical fraction per fate. On the right-hand panels, colored bars represent the literature intervals, and the black ticks show the simulated ; a green or red mark indicates whether . Results summary Across the three stresses ( S =0, 8, 10) the validation is globally satisfactory. At S =0 the simulation is purely cycling (CY= 1.00), matching the high-end of the literature range for viable cycling cells and confirming the absence of stress-induced activation. No senescent or apoptotic response emerges at baseline. At S =8, the distribution becomes mixed, with SEN ≈ 0.50, ARR ≈ 0.31, and a small apoptotic fraction (APO ≈ 0.03). All simulated fractions fall within or touch their literature intervals ( PASS ), indicating a realistic coexistence of transient arrest and senescence. The scatter panel confirms the coexistence of multiple fates within the same ( S , Σ) range. At S =10, the population shifts further towards SEN and APO (SEN ≈ 0.43, APO ≈ 0.38), with residual ARR ≈ 0.19. All points again lie within the admissible colored bands, giving a full PASS classification. This progression (CY → SEN/APO) reproduces the canonical U2OS sequence: cycling cells dominate at low stress, mixed senescence/arrest at moderate stress, and dual senescence–apoptosis responses at high stress. Interpretive remarks Two aspects are worth noting. (i) The pure-CY outcome at S =0 slightly underestimates the minor non-cycling tail allowed by literature, but this likely stems from idealized initial conditions: zero basal p 53, p 21, and PUMA, no micro-stress jitter, and narrow STAT1 dispersion suppress early threshold crossings. (ii) At intermediate stress ( S =8), the simulated fractions sit near the edges of their admissible ranges (ARR at upper bound, SEN at lower), indicating that the model’s decision surfaces are nearly tangent to empirical limits—consistent with the modest normalized gap G S ≈ 0.2. Overall, the U2OS-like population exhibits quantitative agreement with the literature intervals across the entire stress range, with only minor edge-touching effects around the intermediate regime where senescence and arrest overlap. Funder Information Declared ARES Footnotes rebecca.walo{at}unikin.ac.cd (Rebecca Walo Omana), apondondo{at}gmail.com (Apollinaire Ndondo Mboma), didierkumwimba{at}gmail.com (Didier Kumwimba), didier.gonze{at}ulb.be (Didier Gonze) References ↵ Abbas , T. , Dutta , A. , 2009 . p21 in cancer: intricate networks and multiple activities . Nat. Rev. Cancer 9 , 400 – 414 . OpenUrl CrossRef PubMed Web of Science ↵ Abou-Jaoudé , W. , Ouattara , D.A. , Kaufman , M. , 2009 . From structure to dynamics: Frequency tuning in the p53–Mdm2 network I. Logical approach . J. Theor. Biol . 258 , 561 – 577 . OpenUrl CrossRef PubMed ↵ Bae , S.M. , Park , S.J. , Choi , M. , Song , M. , Cho , Y.E. , Do , E.-J. , Ryu , Y.-M. , Park , S. , Lee , B.-H. , Lee , S.-W. , Hwang , S.W. , Park , S.H. , Yang , D.-H. , Ye , B.D. , Byeon , J.-S. , Yang , S.-K. , Joo , J. , Kim , S.-Y. , Myung , S.-J. , 2018 . PSP1, a phosphatidylserine-recognizing peptide, is useful for visualizing radiation-induced apoptosis in colorectal cancer in vitro and in vivo . Transl. Oncol . 11 , 1044 – 1052 . OpenUrl PubMed ↵ Bar-Or , R.L. , Maya , R. , Segel , L.A. , Alon , U. , Levine , A.J. , Oren , M. , 2000 . Generation of oscillations by the p53–Mdm2 feedback loop . Proc. Natl. Acad. Sci. USA 97 , 11250 – 11255 . OpenUrl Abstract / FREE Full Text ↵ Barr , A.R. , Cooper , S. , Heldt , F. , Butera , F. , Stoy , H. , Mansfeld , J. , Novák , B. , Bakal , C. , 2017 . DNA damage during S-phase mediates the proliferation–quiescence decision in the subsequent G1 via p21 . Nat. Commun . 8 , 14728 . OpenUrl CrossRef PubMed ↵ Batchelor , E. , Mock , C.S. , Bhan , I. , Loewer , A. , Lahav , G. , 2008 . Recurrent initiation: a mechanism for triggering p53 pulses in response to DNA damage . Mol. Cell 30 , 277 – 289 . OpenUrl CrossRef PubMed Web of Science ↵ Batchelor , E. , Loewer , A. , Lahav , G. , 2011 . Stimulus-dependent dynamics of p53 in single cells . Mol. Syst. Biol . 7 , 488 . OpenUrl Abstract / FREE Full Text ↵ Becker , V. , Schilling , M. , Bachmann , J. , et al. , 2010 . Covering a broad dynamic range: information processing at the erythropoietin receptor . Science 328 , 1404 – 1408 . OpenUrl Abstract / FREE Full Text ↵ Bravatà , V. , Minafra , L. , Russo , G. , Forte , G.I. , Cammarata , F.P. , Ripamonti , M. , Candiano , G. , Cervello , M. , Giannitrapani , L. , Mazzoccoli , G. , Fazio , I. , Pellerito , O. , Agostara , B. , 2015 . High-dose ionizing radiation regulates gene expression patterns of the immune response, extracellular matrix and DNA repair in Intraoperative Electron Radiotherapy (IOERT) . Anticancer Res . 35 , 2577 – 2591 . OpenUrl Abstract / FREE Full Text Brooks , C.L. , Gu , W. , 2006 . p53 Ubiquitination: Mdm2 and Beyond . Mol. Cell 21 ( 3 ), 307 – 315 . OpenUrl CrossRef PubMed Web of Science ↵ Brzostek-Racine , S. , Gordon , C. , Van Scoy , S. , Reich , N.C. , 2011 . The DNA damage response induces IFN . J. Immunol . 187 , 5336 – 5345 . OpenUrl Abstract / FREE Full Text ↵ Ciliberto , A. , Novak , B. , Tyson , J.J. , 2005 . Steady states and oscillations in the p53/Mdm2 network . Cell Cycle 4 , 488 – 493 . OpenUrl CrossRef PubMed Web of Science ↵ Coleman , M.L. , Marshall , C.J. , Olson , M.F. , 2003 . Ras promotes p21Waf1/Cip1 protein stability via a cyclin D1–imposed block in proteasome-mediated degradation . EMBO J . 22 , 2036 – 2046 . OpenUrl Abstract / FREE Full Text ↵ Dimitrio , L. , 2012 . Modelling nucleocytoplasmic transport with application to the intracellular dynamics of the tumor suppressor protein p53 . Ph.D. Thesis, Université Pierre et Marie Curie (Paris VI) & Università Sapienza Rome . ↵ El-Deiry , W.S. , Tokino , T. , Velculescu , V.E. , et al. , 1993 . WAF1, a potential mediator of p53 tumor suppression . Cell 75 , 817 – 825 . OpenUrl CrossRef PubMed Web of Science ↵ Eliaš , J. , Dimitrio , L. , Clairambault , J. , Natalini , R. , 2014 . The dynamics of p53 in single cells: physiologically based ODE and reaction–diffusion PDE models . Phys. Biol . 11 , 045001 . OpenUrl PubMed ↵ Feng , X. , Tubbs , A. , Zhang , C. , Tang , M. , Sridharan , S. , Wang , C. , Jiang , D. , Su , D. , Zhang , H. , Chen , Z. , Nie , L. , Xiong , Y. , Huang , M. , Nussenzweig , A. , Chen , J. , 2020 . ATR inhibition potentiates ionizing radiation–induced interferon response via cytosolic nucleic acid–sensing pathways . EMBO J . 39 , e104036 . OpenUrl CrossRef PubMed ↵ Geva-Zatorsky , N. , Rosenfeld , N. , Itzkovitz , S. , Milo , R. , Sigal , A. , Dekel , E. , Yarnitzky , T. , Liron , Y. , Polak , P. , Lahav , G. , Alon , U. , 2006 . Oscillations and variability in the p53 system . Mol. Syst. Biol . 2 , 2006 .0033. OpenUrl CrossRef PubMed ↵ Giles , E.M. , Sanders , T.J. , McCarthy , N.E. , et al. , 2017 . Regulation of human intestinal T-cell responses by type I interferon–STAT1 signaling is disrupted in inflammatory bowel disease . Mucosal Immunol . 10 , 184 – 193 . OpenUrl CrossRef PubMed ↵ Giovannozzi , S. , Boccitto , M. , Gattazzo , C. , et al. , 2020 . Live cell imaging demonstrates multiple routes toward a STAT1 gain-of-function phenotype . Front. Immunol . 11 , 1114 . OpenUrl PubMed ↵ Hanson , R.L. , Porter , J.R. , Batchelor , E. , 2019 . Protein stability of p53 targets determines their temporal expression dynamics in response to DNA damage . Cell Rep . 26 , 3709 – 3718.e6 . OpenUrl CrossRef PubMed ↵ Harper , J.W. , Adami , G.R. , Wei , N. , Keyomarsi , K. , Elledge , S.J. , 1993 . The p21 Cdk-interacting protein Cip1 is a potent inhibitor of G1 cyclin-dependent kinases . Cell 75 , 805 – 816 . OpenUrl CrossRef PubMed Web of Science ↵ Heitmeier , M.R. , Scarim , A.L. , Corbett , J.A. , 1999 . Prolonged STAT1 activation is associated with interferon-γ-induced cytotoxicity in islets of Langerhans . J. Biol. Chem . 274 , 8338 – 8344 . OpenUrl ↵ Hunziker , A. , Jensen , M.H. , Krishna , S. , 2010 . Stress-specific response of the p53–Mdm2 feedback loop . BMC Syst. Biol . 4 , 94 . OpenUrl CrossRef PubMed ↵ Jeffers , J.R. , Parganas , E. , Lee , Y. , Yang , C. , Wang , J. , Brennan , J. , MacLean , K.H. , Han , J. , Chittenden , T. , Ihle , J.N. , McKinnon , P.J. , Cleveland , J.L. , Zambetti , G.P. , 2003 . PUMA is an essential mediator of p53-dependent and -independent apoptotic pathways . Cancer Cell 4 , 321 – 328 . OpenUrl CrossRef PubMed Web of Science ↵ Jin , Y. , Lee , H. , Zeng , S.X. , Dai , M.S. , Lu , H. , 2003 . MDM2 promotes p21Waf1/Cip1 proteasomal turnover independently of ubiquitylation . EMBO J . 22 , 6365 – 6377 . OpenUrl Abstract / FREE Full Text ↵ Kalliara , E. , Lázár , E. , Morris , R.H. , et al. , 2022 . Post-transcriptional regulatory feedback encodes JAK/STAT desensitisation after transient IFN-γ . Front. Immunol . 13 , 947213 . OpenUrl PubMed ↵ Kracikova , M. , Akiri , G. , George , A. , Sachidanandam , R. , Aaronson , S.A. , 2013 . A threshold mechanism mediates p53 cell-fate decision between growth arrest and apoptosis . Cell Death Differ . 20 ( 4 ), 576 – 588 . OpenUrl CrossRef PubMed Web of Science ↵ Lahav , G. , Rosenfeld , N. , Sigal , A. , Geva-Zatorsky , N. , Levine , A.J. , Elowitz , M.B. , Alon , U. , 2004 . Dynamics of the p53–Mdm2 feedback loop in individual cells . Nat. Genet . 36 , 147 – 150 . OpenUrl CrossRef PubMed Web of Science ↵ Levy , D.E. , Darnell , J.E. , 2002 . STATs: transcriptional control and biological impact . Nat. Rev. Mol. Cell Biol . 3 , 651 – 662 . OpenUrl CrossRef PubMed Web of Science ↵ Marié , I.J. , Liao , H.-K. , Wang , Y. , et al. , 2021 . Tonic interferon restricts pathogenic IL-17–driven inflammation . eLife 10 , e68371 . OpenUrl PubMed ↵ Meissl , K. , Macho-Maschler , S. , Müller , M. , Strobl , B. , 2017 . The good and the bad faces of STAT1 in solid tumours . Semin. Cancer Biol . 45 , 12 – 20 . OpenUrl Michael , D. , Oren , M. , 2003 . The p53–Mdm2 module and the ubiquitin system . Semin. Cancer Biol . 13 ( 1 ), 49 – 58 . OpenUrl CrossRef PubMed Web of Science Nag , S. , Qin , J. , 2013 . The MDM2-p53 pathway revisited . J. Biomed. Res . 27 ( 4 ), 254 – 271 . OpenUrl CrossRef PubMed ↵ Nakano , K. , Vousden , K.H. , 2001 . PUMA, a novel proapoptotic gene, is induced by p53 . Mol. Cell 7 , 683 – 694 . OpenUrl CrossRef PubMed Web of Science ↵ Ouattara , D.A. , Abou-Jaoudé , W. , Kaufman , M. , 2010 . From structure to dynamics: Frequency tuning in the p53–Mdm2 network II. Differential and stochastic approaches . J. Theor. Biol . 264 , 1177 – 1189 . OpenUrl CrossRef PubMed Web of Science ↵ Proctor , C.J. , Gray , D.A. , 2008 . Explaining oscillations and variability in the p53–Mdm2 system . BMC Syst. Biol . 2 , 75 . OpenUrl CrossRef PubMed ↵ Purvis , J.E. , Karhohs , K.W. , Mock , C. , Batchelor , E. , Loewer , A. , Lahav , G. , 2012 . p53 dynamics control cell fate . Science 336 , 1440 – 1444 . OpenUrl Abstract / FREE Full Text ↵ Schwanhäusser , B. , Busse , D. , Li , N. , Dittmar , G. , Schuchhardt , J. , Wolf , J. , Chen , W. , Selbach , M. , 2011 . Global quantification of mammalian gene expression control . Nature 473 , 337 – 342 . OpenUrl CrossRef PubMed Web of Science ↵ Stewart-Ornstein , J. , Lahav , G. , 2017 . p53 dynamics in response to DNA damage vary across cell lines and are shaped by DNA repair efficiency and ATM activity . Sci. Signal . 10 , eaah6671 . OpenUrl Abstract / FREE Full Text Stewart-Ornstein , J. , Cheng , H.W.J. , Lahav , G. , 2017 . Conservation and divergence of p53 oscillation dynamics across species . Cell Syst . 5 , 410 – 417 . OpenUrl PubMed ↵ Sun , T. , Mu , X. , Cui , J. , 2009 . Modeling the role of p53 pulses in DNA damage-induced cell death decision . BMC Bioinformatics 10 , 190 . OpenUrl CrossRef PubMed ↵ Townsend , P.A. , Scarabelli , T.M. , Davidson , S.M. , Knight , R.A. , Latchman , D.S. , Stephanou , A. , 2004 . STAT-1 interacts with p53 to enhance DNA damage-induced apoptosis . J. Biol. Chem . 279 ( 7 ), 5811 – 5820 . OpenUrl Abstract / FREE Full Text ↵ Vousden , K.H. , Ryan , K.M. , 2009 . p53 and metabolism . Nat. Rev. Cancer 9 , 691 – 700 . OpenUrl CrossRef PubMed Web of Science ↵ Wang , P. , Wang , H.-Y. , Gao , X.-J. , Zhu , H.-X. , Zhang , X.-P. , Liu , F. , Wang , W. , 2023 . Encoding and decoding of p53 dynamics in cellular response to stresses . Cells 12 , 490 . OpenUrl Xiong , Y. , Hannon , G.J. , Zhang , H. , Casso , D. , Kobayashi , R. , Beach , D. , Yu , Q. , Katlinskaya , Y.V. , Carbone , C.J. , et al. , 2015 . DNA-damage-induced type I interferon promotes senescence and inhibits stem cell function . Cell Rep . 11 , 785 – 797 . OpenUrl CrossRef PubMed ↵ Zhang , T. , Brazhnik , P. , Tyson , J.J. , 2007 . Exploring mechanisms of the DNA-damage response: p53 pulses and their possible relevance to apoptosis . Cell Cycle 6 , 85 – 94 . OpenUrl CrossRef PubMed Web of Science References [1]. Abuetabh , Y. , Wu , H.H. , Chai , C. et al : DNA damage response revisited: the p53 family and its regulators provide endless cancer therapy opportunities . Exp Mol Med 54 , 1658 – 1669 ( 2022 ). doi: 10.1038/s12276-022-00863-4 . OpenUrl CrossRef PubMed [2]. Michael , D. , Oren , M. : The p53–Mdm2 module and the ubiquitin system . Semin Cancer Biol . 13 ( 1 ): 49 – 58 ( 2003 ). OpenUrl CrossRef PubMed Web of Science [3]. Nag , S. , Qin , J. : The MDM2-p53 pathway revisited . J Biomed Res . 27 ( 4 ): 254 – 271 ( 2013 ). OpenUrl CrossRef PubMed [4]. Cheng , Q. , Chen , J. : Mechanism of p53 stabilization by ATM after DNA damage . Cell Cycle . 9 ( 3 ): 472 – 478 ( 2010 ). doi: 10.4161/cc.9.3.10556 . OpenUrl CrossRef PubMed Web of Science [5]. Choudhary , H.B. , Mandlik , S.K. , Mandlik , D.S. : Role of p53 suppression in the pathogenesis of hepatocellular carcinoma . World J Gastrointest Pathophysiol . 14 ( 3 ): 46 – 70 ( 2023 ). doi: 10.4291/wjgp.v14.i3.46 . OpenUrl CrossRef PubMed [6]. Hirao , A. et al : DNA damage-induced activation of p53 by the checkpoint kinase Chk2 . Science . 287 ( 5459 ): 1824 – 1827 ( 2000 ). doi: 10.1126/science.287.5459.1824 . OpenUrl Abstract / FREE Full Text [7]. Bar-Or , R.L. , Maya , R. , Segel , L.A. , Alon , U. , Levine , A.J. , Oren , M. : Generation of oscillations by the p53–Mdm2 feedback loop . Proc Natl Acad Sci USA . 97 ( 21 ): 11250 – 11255 ( 2000 ). OpenUrl Abstract / FREE Full Text [8]. Lahav , G. et al : Dynamics of the p53–Mdm2 feedback loop in individual cells . Nat Genet . 36 ( 2 ): 147 – 150 ( 2004 ). OpenUrl CrossRef PubMed Web of Science [9]. Abbas , T. , Dutta , A. : p21 in cancer: intricate networks and multiple activities . Nat Rev Cancer . 9 ( 6 ): 400 – 414 ( 2009 ). OpenUrl CrossRef PubMed Web of Science [10]. Stewart-Ornstein , J. , Lahav , G. : Dynamics of CDKN1A in single cells defined by an endogenous fluorescent tagging toolkit . Cell Rep . 19 ( 8 ): 1965 – 1976 ( 2017 ). OpenUrl [11]. Townsend , P.A. et al : STAT-1 Interacts with p53 to Enhance DNA Damage-induced Apoptosis . J Biol Chem . 279 ( 7 ): 5811 – 5820 ( 2004 ). OpenUrl Abstract / FREE Full Text [12]. Geva-Zatorsky , N. , Rosenfeld , N. , Itzkovitz , S. , Milo , R. , Sigal , A. , Dekel , E. , Yarnitzky , T. , Liron , Y. , Polak , P. , Lahav , G. , Alon , U. : Oscillations and variability in the p53 system . Mol Syst Biol . 2 : 2006 .0033 ( 2006 ). OpenUrl CrossRef PubMed [13]. Batchelor , E. , Loewer , A. , Lahav , G. : Stimulus-dependent dynamics of p53 in single cells . Mol Syst Biol . 7 : 488 ( 2011 ). OpenUrl Abstract / FREE Full Text [14]. Jeffers , J.R. et al : Puma is an essential mediator of p53-dependent and -independent apoptotic pathways . Cancer Cell . 4 ( 4 ): 321 – 328 ( 2003 ). OpenUrl CrossRef PubMed Web of Science [15]. Purvis , J.E. , Karhohs , K.W. , Mock , C. , et al : p53 dynamics control cell fate . Science . 342 ( 6160 ): 1355 – 1359 ( 2013 ). OpenUrl [16]. Kracikova , M. et al : A threshold mechanism mediates p53 cell fate decision between growth arrest and apoptosis . Cell Death Differ . 20 ( 4 ): 576 – 588 ( 2013 ). OpenUrl CrossRef PubMed Web of Science [17]. Ciliberto , A. , Novak , B. , Tyson , J.J. : Steady states and oscillations in the p53/Mdm2 network . Cell Cycle . 4 ( 3 ): 488 – 493 ( 2005 ). OpenUrl CrossRef PubMed Web of Science [18]. Spencer , S.L. , Gaudet , S. , Albeck , J.G. , Burke , J.M. , Sorger , P.K. : Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis . Nature . 459 : 428 – 432 ( 2009 ). doi: 10.1038/nature08012 . OpenUrl CrossRef PubMed Web of Science [19]. Youlyouz-Marfak , N. et al : Identification of a novel p53-dependent activation pathway of STAT1 by antitumour genotoxic agents . Cell Death Differ . 15 ( 3 ): 376 – 385 ( 2008 ). OpenUrl CrossRef PubMed Web of Science [20]. Haronikova , L. et al : The p53 mRNA: an integral part of the cellular stress response . Nucleic Acids Res . 47 ( 7 ): 3257 – 3271 ( 2019 ). doi: 10.1093/nar/gkz124 . OpenUrl CrossRef PubMed [21]. Hu , W. , Feng , Z. , Levine , A.J. : The Regulation of Multiple p53 Stress Responses is Mediated through MDM2 . Genes Cancer . 3 ( 3-4 ): 199 – 208 ( 2012 ). doi: 10.1177/1947601912454734 . OpenUrl CrossRef PubMed [22]. Takagi , M. , Absalon , M.J. , McLure , K.G. , Kastan , M.B. : Regulation of p53 Translation and Induction after DNA Damage by Ribosomal Protein L26 and Nucleolin . Cell . 123 ( 1 ): 49 – 63 ( 2005 ). OpenUrl CrossRef PubMed Web of Science [23]. Brooks , C.L. , Gu , W. : p53 Ubiquitination: Mdm2 and Beyond . Mol Cell . 21 ( 3 ): 307 – 315 ( 2006 ). OpenUrl CrossRef PubMed Web of Science [24]. Pant , V. , Lozano , G. : Limiting the power of p53 through the ubiquitin proteasome pathway . Genes Dev . 28 ( 17 ): 1739 – 1751 ( 2014 ). doi: 10.1101/gad.247452.114 . OpenUrl Abstract / FREE Full Text [25]. Hunziker , A. , Jensen , M.H. , Krishna , S. : Stress-specific response of the p53Mdm2 feedback loop . BMC Syst Biol . 4 : 94 ( 2010 ). doi: 10.1186/1752-0509-4-94 . OpenUrl CrossRef PubMed [26]. Proctor , C.J. , Gray , D.A. : Exploring the mechanisms of p53 oscillations in a model of DNA damage . J Theor Biol . 253 ( 2 ): 323 – 338 ( 2008 ). OpenUrl CrossRef PubMed [27]. Maltzman , W. , Czyzyk , L. : UV irradiation stimulates levels of p53 cellular tumor antigen in non-malignant mouse cells . J Biol Chem . 259 ( 4 ): 2451 – 2457 ( 1984 ). OpenUrl [28]. Price , B.D. , Calderwood , S.K. : G1 arrest by UV irradiation is associated with increased levels of p53 and reduced cyclin B1 in hamster and human cells . J Biol Chem . 268 ( 16 ): 11927 – 11931 ( 1993 ). OpenUrl [29]. Meek , D.W. , Simon , S. , Kikkawa , O. : Modulation of p53 and its effector proteins in response to DNA damage . Oncogene . 21 ( 52 ): 9032 – 9046 ( 2002 ). OpenUrl [30]. Loewer , A. , Batchelor , E. , Gaglia , G. , Lahav , G. : Basal dynamics of p53 reveal transcriptionally attenuated pulses in cycling cells . Mol Cell . 40 ( 3 ): 367 – 380 ( 2010 ). OpenUrl [31]. Azzam , E.I. , Jay-Gerin , J.P. , Pain , D. : Ionizing radiation-induced metabolic oxidative stress and prolonged cell injury . Cancer Lett . 327 ( 1-2 ): 48 – 60 ( 2012 ). doi: 10.1016/j.canlet.2011.12.012 . OpenUrl CrossRef PubMed Web of Science [32]. Purvis , J.E. , Lahav , G. : Encoding and decoding cellular information through signaling dynamics . Cell . 152 ( 5 ): 945 – 956 ( 2013 ). doi: 10.1016/j.cell.2013.02.005 . OpenUrl CrossRef PubMed Web of Science [33]. Stewart-Ornstein , J. , et al : Cell-cycle dynamics of p53 in response to DNA damage . Proc Natl Acad Sci USA . 114 ( 29 ): E6020 – E6029 ( 2017 ). doi: 10.1073/pnas.1705643114 . OpenUrl Abstract / FREE Full Text [34]. Schröder , K. , Hertzog , P.J. , Ravasi , T. , Hume , D.A. : Interferon-γ: an overview of signals, mechanisms and functions . J Leukoc Biol . 75 ( 2 ): 163 – 189 ( 2004 ). doi: 10.1189/jlb.0603252 . OpenUrl CrossRef PubMed Web of Science [35]. Riley , J.S. , Tait , S.W.G. : Mitochondrial regulation of cell death: a new dawn . Nat Rev Mol Cell Biol . 20 ( 12 ): 721 – 737 ( 2019 ). OpenUrl CrossRef PubMed [36]. Hernandez-Segura , A. , Nehme , J. , Demaria , M. : Hallmarks of cellular senescence . Trends Cell Biol . 28 ( 6 ): 436 – 453 ( 2018 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted October 30, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following STAT1–p53 Dynamics Program Cell Fate through p21 and PUMA Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share STAT1–p53 Dynamics Program Cell Fate through p21 and PUMA M.K. Francis Didier Tshianyi , Rebecca Walo Omana , Apollinaire Ndondo Mboma , Didier Kumwimba , Didier Gonze bioRxiv 2025.10.28.685172; doi: https://doi.org/10.1101/2025.10.28.685172 Share This Article: Copy Citation Tools STAT1–p53 Dynamics Program Cell Fate through p21 and PUMA M.K. Francis Didier Tshianyi , Rebecca Walo Omana , Apollinaire Ndondo Mboma , Didier Kumwimba , Didier Gonze bioRxiv 2025.10.28.685172; doi: https://doi.org/10.1101/2025.10.28.685172 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41970) Biophysics (21463) Cancer Biology (18605) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19911) Epidemiology (2067) Evolutionary Biology (24329) Genetics (15615) Genomics (22514) Immunology (17743) Microbiology (40424) Molecular Biology (17194) Neuroscience (88650) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4827) Physiology (7648) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)