Excitability in the p53 network mediates robust signaling with tunable activation thresholds in single cells

Cellular signaling systems precisely transmit information in the presence of molecular noise while retaining flexibility to accommodate the needs of individual cells. To understand design principles underlying such versatile signaling, we analyzed the response of the tumor suppressor p53 to varying levels of DNA damage in hundreds of individual cells and observed a switch between distinct signaling modes characterized by isolated pulses and sustained oscillations of p53 accumulation. Guided by dynamic systems theory we show that this requires an excitable network structure comprising positive feedback and provide experimental evidence for its molecular identity. The resulting data-driven model reproduced all features of measured signaling responses and is sufficient to explain their heterogeneity in individual cells. We present evidence that heterogeneity in the levels of the feedback regulator Wip1 sets cell-specific thresholds for p53 activation, providing means to modulate its response through interacting signaling pathways. Our results demonstrate how excitable signaling networks can provide high specificity, sensitivity and robustness while retaining unique possibilities to adjust their function to the physiology of individual cells.

-Wavelet based peak detection method A An exemplary single cell p53 trajectory, the detected pulses are color coded. B The wavelet transform corresponding to the trajectory shown in A. C The ridge lines extracted from the wavelet transform (B). The color-coding corresponds to the pulses in (A). See also material and methods section in the main text.

Figure S2 -Phase portrait of the excitable toy model
Three subthreshold trajectories, which rapidly settle back to the stable fixed point, and three excitation loops (pulses) making the full phase space excursion are shown for the NPF system. The direction dependent threshold visibly originates from the saddle point and separates the phase space into subthreshold and excitable dynamics. It is called a separatrix and coincides with the stable manifold of the saddle point. Such a phase space structure with a saddle as organizing center is typical for type I excitable systems.

Figure S3 -qPCR quantification of potential positive feedbacks ATM inhibitor studies
Time resolved mean expression of p53 (A), the known negative feedback regulator Wip1 (B), the effector p21 (C) and potential transcriptional positive feedbacks (D) as indicated on the y-axis. Timing and fold change of the three potential feedbacks PIDD, 14-3-3 sigma and PTEN do not support an effect on p53 pulse formation. Error bars indicate standard deviation.

Figure S4 -ATM inhibitor studies
A Differential effect of signal vs. species X inhibition for the generic negativepositive feedback (NPF) system presented in Figure 2 of the main text. Signal inhibition leads to an "all-or-none" systems response. X inhibition leads to a time-of-inhibition dependent collapse of the excitation loop, reflected by the smaller response amplitudes. B Single cell p53 first pulse amplitudes after strong stimulation (10Gy) and kinase inhibitor Wortmannin addition at time points as indicated. The earlier ATM is inhibited, the smaller are the detected p53 pulses. C Fraction of responsive cells for the experimental conditions described in F. Early ATM inhibition abrogates the p53 pulse formation for a substantial amount of cells below the detection limit. Figure S5 -Bifurcation analysis of the p53 model A One parameter bifurcation analysis. The control parameter is the stimulus strengths (S). For low stimulation the system is in the excitable regime. Upon the simultaneous homoclinic (SNIC) and saddle-node bifurcation (LP1) a limit cycle suddenly appears with full amplitude. The arrow marks the point in parameter space for the phase portrait shown in B See also Supplementary Text. B Phase portrait for the p53 model in the excitable regime as marked in A. Shown are five subthreshold trajectories which rapidly settle back to the stable fixed point and three excitation loops (pulses) making the full phase space excursion. Neither high or low Wip1 concentrations alone can trigger a pulse, only strong ATM activation crosses the threshold. C Two parameter bifurcation analysis, the first control parameter is the stimulus strengths (S), the second parameter Wip1 protein maturation rate (TW). The Cusp Point (CP) spawns two saddle-node curves which enclose the excitable regime. A Bogdanov-Takens Point (BT) is very close to the Cusp Point and spawns a Hopf curve enclosing the oscillatory regime. For the one fixed-point regime a response was classified as excitable if it shows at least a five-fold ATM response upon a small initial perturbation. See also Supplementary Text.

Figure S6 -Local sensitivity analysis of the deterministic core model
Model parameter variations are given as log2 fold change from the standard values given in supplementary section 4. The red bars indicate the extent of the excitable regime in parameter space. It was determined by stimulating the system with a perturbation in ATM * and checking for a pulsatile response. Within the excitable regime, the maximal and minimal effect of the respective parameter variations on the p53 pulse amplitude (width) is indicated in blue (yellow). Parameters associated with the Wip1 phosphatase show the highest sensitivity, which underlines its proposed function as a major modulator of the p53 signaling pathway.  C Western blot analysis of ATM and Chk2 phosphorylation upon DNA damage in A549 cells overexpressing Wip1 compared to control cells. Quantification was performed by densitometry; signal intensities were corrected for unequal loading using GAPDH and normalized to undamaged control cells (Control 0h). D Wip1 induction for highly (10Gy) stimulated cells. The number of pulses is clearly lower for cells with increased Wip1 expression. E Mean Wip1 mRNA levels upon transfection of scrambled control and Wip1targeting siRNA were measured by quantitative RT-PCR (3d after transfection). Wip1 expression was downresgulated below 20% of control levels. Error bars indicate standard deviation.  The negative-feedback (NF) system presented in Figure 2 of the main text is given by: This network was inspired by simple models used to describe the p53-Mdm2 core negative feedback loop, for example as in Ref. [1]. The variable y describes an mRNA species whose production is activated by a protein X acting as transcription factor. The negative feedback closes with the degradation of X by the maturated protein Y . The negative-positive feedback (NPF) system is given by the equations: Here species X exhibits a self-activating term with a Hill coeffecient of two acting as the positive feedback, e.g. describing processes like auto-phosphorylation. This term alone would introduce bistability into the system. The additional negative feedback is topologically the same as for the NF system presented above, albeit the degradation of X by Y follows mass action kinetics. This NPF system resembles a so-called bistable frustrated unit as discussed by the authors of Ref. [2]. The parameters for both models are given as tables in section 2.
2 Tables of parameters for the toy models  Table 2: Positive-negative feedback system model parameters

p5model equations
The core deterministic p53 model was formulated as a system of coupled ordinary differential equations (ODEs). The associated regulatory network is shown in Figure 3A of the main text. We modeled only the amount of activated ATM (AT M ⇤ ), where the positive feedback described in the main text is modeled phenomenologically by a self-activation term with a Hill coeffecient of two. The phosphatase Wip1 not only dephosphorylates AT M ⇤ directly, but also other species (i.e. H2AX, MRN complex) which are involved in the positive feedback [3]. This is modeled by an inhibition of the AT M ⇤ self-activation term by the Wip1 protein. The constitutive expression of p53 ( Figure S2A) was described by the constant production rate C for the protein P 53. The p53 protein itself is a target of AT M ⇤ and Mdm2 is reported to have a lower binding affinity to phosphorylated P 53 . As with the present data phosphorylated p53 was not quantified, no additional species was introduced into the model. Therefore, this mechanism was also indirectly captured phenomenologically by inhibiting the Mdm2 dependent P 53 degradation via AT M ⇤ . The relative strength of this second AT M ⇤ interaction is given by the parameter R. Transcriptional activation of the two mRNA species mdm2 and wip1 by P 53 were modeled by saturating kinetics. All other terms of the ODE system follow mass-action kinetics or first-order kinetics. The dependence of the effective input signal S on the number of DSBs is modeled by a saturating kinetic to control both the maximal stimulation and the sensitivity towards DSBs. For the stochastic forcing ( Figure 3D,E and F, main text) the number of DSBs is time-dependent and given by the stochastic DSB process defined below in section 5. The naming of the parameters follows general conventions: parameters containing the letter "T" denote production rates concerning translation or transcription, parameters starting with a "d" denote degradation rates and parameters containing "k" are the respective Michaelis constants. The subscripts encode the affiliation to the modeled species, e.g. a "m" in the subscript is associated with the Mdm2 mRNA and a "M" with the respective protein species. Thus the parameter k P m , for example, describes the Michaelis constant of the Mdm2 mRNA production induced by the p53 protein. Rate parameters which do not follow this nomenclature are: A which gives the maximal activation rate of ATM, P which describes the dephosphorylation of AT M ⇤ by W ip1, g which gives the maximal Mdm2 dependent degradation of P 53 and finally C which is the maturation rate of new P 53 entering the system. An overview of all model parameters and their values is given in the table 4 below. d dt  Michaelis constant for the production of mdm2 by P 53 1.

Detailed description of the Markovian DSB process
The dynamics of the DNA double strand breaks (DSBs) were considered as the primary input of the p53 model. Fortunately, as outlined in the main text, single cell trajectories of damage foci are available. While the direct relation between foci number and number of DSBs remains unclear, they can still serve as a quantitative marker [4]. Therefore, the number of foci was treated as the number of DSBs. Typically the DSB dynamics are described by an exponential model [5], which sufficiently describes the repair process. However, the half-lifes of the foci show large variability across a cell population [4,5]. The sources of the observed stochasticity of the DSB dynamics lie both within the processes which cause the DNA damage as well as the cellular repair mechanisms. The former include e.g. the errors arising from mitotic division, the occurrence of radical metabolic byproducts and cosmic radiation, which are intrinsically noisy. The latter appear to be irregular also due to various reasons, e.g. different severity of the DNA lesions or different copy numbers and spatial availability of the proteins involved in the repair process.

The time-homogeneous repair process
A simple approach to model the stochastic repair process is to define two rates b and r by the following assignments of probabilities, given that there are n DSBs present at time t: b dt = probability that a new DSB occurs in [t, t + dt] This is a Markovian birth-death process, named the payroll process by Gillespie [6]. The asymptotic analytical solutions for the mean and variance read: These relations only provide the ratio between the two rates b and r. To effectively reverse calculate the rates from the foci data, dynamical properties of the process have to be taken into account. The aim here is not to claim or provide an exact parameter estimation, but to give a reasonable estimation given the simplistic stochastic model used and the data available. In the following, the background damage level N b is used for the estimation, this already gives b = N b r. The first moment time-evolution function for the payroll process given in ref. [6] reads: Substituting with the asymptotic mean yields The notion of N b and hN b i might be a little confusing. For the stochastic process, these two are identical, N b being the asymptotic time average of one realization and hN b i being the asymptotic ensemble average. For the time series foci data, these two hardly coincide. Reasons for that are the finite and rather short period of sampling, some additional cell-to-cell variability and probably some systematic error in the measurements, as the number of foci is only a proxy for the number of DSBs. Nevertheless, to advance with the estimation, the l.h.s. of equation 7 is treated as an ensemble average at time t. Solving for the repair rate one obtains: The time dependence of the rate is obviously artificial, as the rate should be time independent for a time homogeneous process. And indeed, when using this formula to reverse calculate the repair rate out of an ensemble of synthetic data, r is time independent for all 0 < t < t c . An example is shown in figure  S4A and B. However, when using equation 8, one carefully has to observe the denominator. From a critical time on, the ensemble average at time t c gets very close to the asymptotic average: hDSBi tc ⇡ N b . This quickly leads to numerical precision issues, as the denominator converges exponentially fast to zero and might even turn negative due to fluctuations. Practically this means, that the estimation of the repair rate r should be carried out in a time domain, where the trajectories on average are still decaying. To improve the estimate, an averaging within the respective time domain [0, t c ] can be carried out and was done for the inference from the real data. It is noteworthy, that the relaxation dynamics are needed for the parameter estimation. An unstimulated ensemble of trajectories corresponding to stationary dynamics alone would not allow for the described method as this would correspond to t c = 0. Having demonstrated that the proposed method for rate parameter estimation works with synthetic data, a dataset containing cells irradiated with 5Gy of gamma radiation was analyzed accordingly. As the DSB induction for irradiated cells is very fast, only the repair dynamics are observable. Because every foci trajectory has generally a different N 0 , the data is binned to form subensembles with comparable initial amount of DSBs. As can be seen in figure S4C, the initial amount of DSBs present in the cells has a broad distribution despite the population received a fixed damage dose. This makes the parameter estimation even more difficult, due to the low sample numbers in the subensembles. The results for two different subensembles for the 5Gy data set are shown in figure  S4D. Given the broad initial distribution and that there are only 63 cells in total, there is no satisfactory pooling available. But nevertheless, the repair rate was estimated to be in the order of r ⇡ 0.315h 1 per hour, which yields an average lifetime of a damage locus to be around 3 hours ( Figure S4E). This is in good agreement with the results reported experimentally [4,5]. Given the mean damage background level to be hN b i ⇡ 2.3, the corresponding basal break rate for the DSB process is b ⇡ 0.8h 1 .

Modeling the DSB induction by NCS
The dataset used for the main results in the main text consists of cells which were stimulated with the radiomimetic drug neocarcinostatin (NCS). As opposed to irradiated cells, the DSB induction is directly observable and lasts around an hour in these experiments (Figure 3 B and C). To address this phase of rapid DSB induction within our Markovian modeling framework, we augmented the stochastic process with a time dependent break rate: Here b b is the basal break rate as estimated in section 5.1 above, and b s > b b is the NCS dose dependent stimulation break rate. The switching between both rates occurs at T s . To devise a simulation algorithm we first recall the next-jump density function from the classical Gillespie algorithm (SSA) [6]: where the a 0 k s are the propensities. Using integration by parts and employing inversion sampling gives the expression: here r is a random number drawn from the uniform distribution U (0, 1). Solving this equation for ⌧ allows for the SSA algorithm. For the time-homogeneous case a k (n, t) ⌘ a k (n) this is straightforward: where for the last expression we substituted the propensities for our homogeneous DSB process. For arbitrary time-inhomogeneous processes equation 11 often has to be solved numerically, which greatly increases the cost of the algorithm. However, as our time-inhomogeneous rate b(t) is a simple step-function, equation 11 can still be solved analytically. Integration yields a family of piecewise linear functions and depending on the systems time t the next jump occurs at: here s = T s t is the distance in systems time to the break rate switch and z = ln(1/r). The additional conditions on z for the first two cases (t  T s ) ensure, that the correct section of the piece wise linear function is inverted.
In summary, augmenting the cost effective SSA algorithm with a test for two conditions per step makes it readily applicable to time dependent rates formulated as step functions. An example realization of the DSB process upon NCS stimulation is shown in Figure S4F. The good agreement of our DSB process with an ensemble of NCS stimulated cells is shown in Figure 3C of the main text.

Bifurcation analysis of the p53 model
The one parameter (or codimension one) bifurcation analysis of the deterministic p53 model is shown in Figure S3A. The control parameter is the signal strengths S which appears on the r.h.s. of the equation for AT M ⇤ (see section 3), describing DSB induced ATM activation. For a low input signal, equivalent to no or very low damage levels, the system is in the excitable regime. This is characterized by the presence of one stable and two unstable fixed points or steady states. The lowest fixed point is a stable spiral, giving rise to small subthreshold oscillations for small perturbations. The middle fixed point is a saddle, whose stable manifold separates the phase space into regions of different qualitative behavior. This separatrix acts as an direction dependent threshold for the system. The bifurcation leading to the oscillations in the p53 model is the saddle-node homoclinic bifurcation [7]. This is a global bifurcation involving a local saddle-node bifurcation. What happens at the bifurcation point is, that the heteroclinic connection of the saddle to the stable fixed point becomes a homoclinic orbit of the merged saddle-node. This saddle-node then disappears via the local saddle-node bifurcation (denoted as limit point in figure S3A) and a limit cycle appears near the former homoclinic orbit (called a saddle node on invariant circle (SNIC) bifurcation). From within the excitable regime, which is considered as operating point of the model for no DSBs present, the oscillatory regime can only be reached by increasing the signal S. The limit cycles of this positive feedback oscillator are born with huge amplitudes, and show only very little dependence on the signal strength up to S ⇡ 0.5. In effect, for no or low signal the system resides in an excitable state capable of showing isolated pulses. In addition, for a stronger signal after the saddle-node homoclinic bifurcation the system undergoes stable sustained oscillations with pulse shapes very similarto excitory pulses ( Figure 3D-F main text). These are exactly the characteristics found in the p53 single cell data for low vs. high damage input, which were discussed in the main text.
To further extend the bifurcation analysis, a search for codimension-2 bifurcations was performed as shown in figure S3B. Here a Cusp bifurcation point (CP) was found. This marks the appearance of the phasespace structure required for the excitability class I regime: the co-existence of a saddle as organizing center, one stable and one unstable fixed point. There is additionally a Bogdanov-Takens (BT) point very close by, spawning a Hopf bifurcation curve which effectively destabilizes the upper fixed point and makes the system excitable. The specific bifurcation parameters used are T W , the Wip1 protein production rate, and again S. However, there are actually only 2 symmetric types of these codimension-2 diagrams, the other type and further examples are presented in Ref. [8]. The actual region of excitability in parameter space extends beyond the three fixed point regime (see figure S3B). Although the two upper fixed points vanished via a saddle-node bifurcation, the phase space still has memory about their presence. What is meant by that is, that the global flow (the excitation loop) is not seriously affected by this local bifurcation. The exact definition of the threshold as separatrix is, however, lost in this regime. The oscillatory regime here is reached via a subcritical Hopf bifurcation fol-lowed by a saddle-node bifurcation of limit cycles. The stable limit cycle again suddenly appears with huge amplitudes.
Notably, besides the excitable and oscillatory regimes there is a bistable regime for very low maximal transcriptional activation of Wip1 and a not too strong signal. This is due to the now weak negative feedback which is facilitated by the Wip1-ATM interaction. Therefore the positive feedback predominates and induces bistability. Interestingly, for medium to high signal strengths the bistability is lost via a saddle-node bifurcation and the system becomes monostable, with the highest steady state now being the only attractor of the system. This resembles the behavior found for Wip1 knockout cells [9].