Single-cell analysis of p53 transitional dynamics unravels stimulus- and cell type-dependent signaling output motifs

Background To understand functional changes of complex biological networks, mathematical modeling of network topologies provides a quantitative measure of the way biological systems adapt to external stimuli. However, systemic network topology-based analysis often generates conflicting evidence depending on specific experimental conditions, leading to a limited mechanistic understanding of signaling networks and their differential dynamic outputs, an example of which is the regulation of p53 pathway responses to different stress stimuli and in variable mammalian cell types. Here, we employ a network motif approach to dissect key regulatory units of the p53 pathway and elucidate how network activities at the motif level generate context-specific dynamic responses. Results By combining single-cell imaging and mathematical modeling of dose-dependent p53 dynamics induced by three chemotherapeutics of distinct mechanism-of-actions, including Etoposide, Nutlin-3a and 5-fluorouracil, and in five cancer cell types, we uncovered novel and highly variable p53 dynamic responses, in particular p53 transitional dynamics induced at intermediate drug concentrations, and identified the functional roles of distinct positive and negative feedback motifs of the p53 pathway in modulating the central p53-Mdm2 negative feedback to generate stimulus- and cell type-specific signaling responses. The mechanistic understanding of p53 network dynamics also revealed previously unknown mediators of anticancer drug actions and phenotypic variations in cancer cells that impact drug sensitivity. Conclusions Our results demonstrate that transitional dynamics of signaling proteins such as p53, activated at intermediate stimulus levels, vary the most between the dynamic outputs of different generic network motifs and can be employed as novel quantitative readouts to uncover and elucidate the key building blocks of large signaling networks. Our findings also provide new insight on drug mediators and phenotypic heterogeneity that underlie differential drug responses. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-022-01290-7.


I. Mathematical model and analysis of the ATM/p53/Mdm2 regulatory module in response to Etoposide
Dynamic output of the p53-Mdm2 motif in response to Etoposide is formulated by the following simplified delay differential equations (DDEs). [ where [ ] denotes dimensionless concentrations of the total proteins (p53, Mdm2, ATMt) or the active, phosphorylated form of ATM (ATMp). Briefly, ATM is phosphorylated and activated by different concentrations of Etoposide with rate constant ea and ATMp is dephosphorylated with rate constant a0 (Equation (1)). ATMp subsequently phosphorylates p53 and Mdm2, leading to decrease in p53-Mdm2 binding and Mdm2-mediated p53 degradation, as described by the ATMpdependent rate parameter mp . p0 and p describe the rate of basal production and Mdm2independent degradation of p53, respectively (Equation (2)). The transcriptional activation of Mdm2 by tetrameric p53 is characterized by a Hill function of 4 th order in Equation (3)  Equations (1)-(3) were derived with the assumption of rapid equilibrium of phosphorylation and dephosphorylation, as they are much faster than transcription and degradation. 3 This assumption allowed us to simplify the kinetic equations for total p53 and Mdm2 without specifying their respective phosphorylated and unphosphorylated forms, as well as to obtain the ATMp dependence of the kinetic rates as follows (refer to the supplementary information of ref. [15] for the detailed derivation).
where p0 ( 0 ) and p0 ( 0 ) are the basal phosphorylation and dephosphorylation rate constants of p53 (Mdm2), ap ( ) are the ATM-mediated phosphorylation rate of p53 (Mdm2), mp0 is the Mdm2-mediated p53 degradation rate constant, pm0 is the Michaelis constant for p53-induced upregulation of Mdm2, and 0 and 1 are the degradation rate constants of the unphosphorylated and phosohorylated Mdm2, respectively.
Similar to the modeling results for the p53-Mdm2 negative feedback motif under Nutlin-3a, p53 levels simulated for the ATM/p53/Mdm2 module produced a unique p53 transitional dynamics at intermediate Etoposide concentration, i.e., an initial large p53 pulse followed by an elevated plateau (supplementary Fig. S1A). This again did not agree with our experimental data for the Etoposide response and suggested additional regulatory components/interactions are involved beyond the ATM/p53/Mdm2 interactions. Figure S1B showed the distribution of parameter values that can generate periodic pulsing of p53 at 1 M Etoposide as well as monotonic p53 induction at high drug dose, as we observed experimentally for the etoposide-sensitive cell 4 lines. The values for rate constants of Mdm2-mediated p53 degradation ( mp0 ), p53-induced Mdm2 production ( pm ) and ATM-mediated Mdm2 degradation ( 1 ) were the most broadly distributed, while the Michaelis constant pm0 spanned a smaller, but still 10-fold range.
As expected, the dynamic output of ATM/p53/Mdm2 is also regulated by the time delay in p53-mediated Mdm2 upregulation, m . By varying m , we found the oscillatory feature of p53 level at low drug dose is evident only when m  0.7 hour (hr) (Fig. S1C, left panel). More specifically, when 0.7 hr  m < 1.7 hr, p53 dynamics exhibit damped oscillation, while 1.7 hr  m < 4.3 hr gives rise to steady sinusoidal p53 oscillation. For m  4.3 hr, p53 oscillation becomes non-sinusoidal, but maintains the periodicity. Across the whole oscillatory regime, period of the p53 oscillation is proportional to the time delay m , and the oscillation amplitude also showed largely linear dependence on m , when m < 4.3 hr (Fig. S1C, right panel). However, varying m again did not alter the transitional dynamics of p53 induction as discussed above.

II. Numerical simulation and parameter space search
The above model and the models discussed in the main text were numerically simulated using the Matlab built-in function dde23. We set the initial concentration of the different protein Tables S1-S4. These specific parameter sets were chosen based on the quantitative features of our single-cell imaging data as well as the western blot data published in a previous study [15]. That is, these parameter sets produced simulation results most similar to our experimental observations. 6