A systematic approach to decipher crosstalk in the p53 signaling pathway using single cell dynamics

Abstract

The transcription factors NF-κB and p53 are key regulators in the genotoxic stress response and are critical for tumor development. Although there is ample evidence for interactions between both networks, a comprehensive understanding of the crosstalk is lacking. Here, we developed a systematic approach to identify potential interactions between the pathways. We perturbed NF-κB signaling by inhibiting IKK2, a critical regulator of NF-κB activity, and monitored the altered response of p53 to genotoxic stress using single cell time lapse microscopy. Fitting subpopulation-specific computational p53 models to this time-resolved single cell data allowed to reproduce in a quantitative manner signaling dynamics and cellular heterogeneity for the unperturbed and perturbed conditions. The approach enabled us to untangle the integrated effects of IKK/ NF-κB perturbation on p53 dynamics and thereby derive potential interactions between both networks. Intriguingly, we find that a simultaneous perturbation of multiple processes is necessary to explain the observed changes in the p53 response. Specifically, we show interference with the activation and degradation of p53 as well as the degradation of Mdm2. Our results highlight the importance of the crosstalk and its potential implications in p53-dependent cellular functions.

Fig 1. Activation status of NF-κB affects the p53 response to DSBs.

A549 reporter cells were tracked via live-cell time-lapse microscopy and the p53 median nuclear fluorescence intensity was measured in cells treated with 10 Gy irradiation (IR) in combination with DMSO (a), TNFα (b) or IKK2i (c). Light colored lines denote single cell trajectories whereas bold lines indicate the median of trajectories. d) Quantification of selected features (timing maxima, timing minima, inter-peak-intervals and dampening factor) of the first four p53 pulses. To test the significance of changes, we used the Wilcoxon rank sum test in combination with the Bonferroni-Holm method to correct for multiple testing, *p<0.05, **p<0.01, ***p<0.001. e) mRNA expression of p53 target genes was measured upon 10 Gy IR in A549 cells treated with DMSO or IKK2i by qRT-PCR. β-actin was used as an internal control. Data were normalized to values of the sample harvested 0 h post IR. Error bars indicate the minimum and the maximum value for the relative quantity of three technical replicates.

Fig 2. Scheme of the subpopulation modeling framework and the model describing the p53 network.

a) Single cell trajectories are clustered into subpopulations with similar dynamics (upper left panel). For each cluster, a peak-based mean is calculated, representing subpopulation-specific dynamics. To account for cellular heterogeneity and to reproduce the different dynamics of subpopulations, a pool of models is generated (upper right panel). In this example, the model pool comprises two models, one is specific for subpopulation a, the other one for subpopulation b. The production rates of Z and X are assumed to be susceptible to noise. Hence, parameters of these two processes are considered subpopulation-specific and therefore specific for an individual model. While k1_a and k2_a are specific for subpopulation a and thus assigned to one model, k1_b and k2_b are specific for subpopulation b and assigned to the second model. Importantly, both models share parameters such as k3, assuming that for example phosphorylation rates or degradation rates are not affected by noise. In a final step, the model pool is simultaneously fitted to the peak-based mean of each subpopulation. b) The ODE model consists of seven variables representing phosphorylated ATM (ATM-P), p53 in its inactive and transcriptionally active (act) form, as well as mRNAs and proteins of Mdm2 and Wip1. Moreover, it contains 22 parameters; those marked in red are defined as subpopulation-specific since the corresponding processes are assumed to be susceptible to noise.

Fig 3. Model pool reproduces heterogeneous p53 dynamics.

a) The ten models of the model pool were fitted to the peak-based mean of the ten subpopulations (blue lines). The red line indicates the simulation of the best fit. The number of cells assigned to a subpopulation is stated in the upper right corner of each subpopulation plot. To visualize the simulation and experimental data on the population level (lower right corner), the weighted mean over all subpopulations was calculated. The weight is determined by the number of assigned cells. b) L1 regularization-based evaluation of subpopulation-specific parameters. The parameters correspond to the six biochemical processes that were assumed to be subpopulation-specific. The processes include basal transcription of the Mdm2 gene (β_mt), translation of Mdm2 mRNA (β_mtm), synthesis of p53 (β_p), activation of ATM (β_s), basal transcription of the Wip1 gene (β_wt) and translation of Wip1 mRNA (β_wtw). Subpopulation-specific parameters are represented by colored dots. Each color denotes an individual subpopulation. Parameters with a value close to zero are identified by L1 regularization as unspecific (circles).

Fig 4. Sensitivity analysis predicts existence of multiple interactions of IKK2 and the p53 network.

a) Expected change in features upon perturbation of a parameter, which indicates an interaction of IKK2 and the p53 network. Each box represents the effect of a parameter perturbation on a specified feature for an individual subpopulation. A red box indicates a consistent increase in the specified feature, in all considered peaks. Blue boxes represent a decrease. b) Sensitivity analysis of single parameters. If the absolute value of sensitivity coefficients of all considered peaks are below 1·10⁻⁴, the specified parameter change has no considerable effect on the respective feature (grey boxes). Subpopulations with sensitivity coefficients of peaks which are not consistently changed are depicted in white. c) For the sensitivity analysis of parameter pairs, 13 pairs were selected from 462 possible combinations based on their ability to reflect the expected change in features in at least seven out of ten subpopulations. Arrows indicate a perturbation of parameters by increasing or decreasing the respective parameter value by 1%.

Fig 5. Validation of the 30 best ranked parameter combination fits using time-variant IKK2 inhibition.

a) A549 reporter cells were tracked and the p53 median nuclear fluorescence intensity was measured in cells treated with DMSO or IKK2i at different time points. Sample trajectories are shown in which IKK2 was inhibited at the indicated time points. In an additional experiment, IKK2 was inhibited at 1.5h, 2.5h and 5h after IR. b) Each box depicts the color coded log₁₀ χ² value for model simulations with individual parameter combinations and time points of IKK2 inhibition. For better visualization, the plotted log₁₀ χ² value summarizes the values of two data sets with similar time points of IKK2 inhibition (1.5h/1.5h, 2.5h/3h, 5h/5h). The inhibitor was applied at accumulating levels of p53, resulting in the first peak (t_Inh = 1.5h), decreasing p53 levels (t_Inh = 2.5h and t_Inh = 3h) as well as at the accumulation resulting in the second peak (t_Inh = 5h). The parameter combinations are sorted based on the corresponding summarized log₁₀ χ² value. Indices of the parameter combinations are given by the numbers on the right hand side of the plot.

Fig 6. IKK2 inhibition leads to a delayed Wip1 and Mdm2 accumulation and increased pChk2 levels on the population level.

Upon 10 Gy IR in A549 cells treated with DMSO or IKK2i, the levels of p53 (a), Mdm2 mRNA (b), Mdm2 protein (c), Wip1 mRNA (d), Wip1 protein (e) and pChk2 (f) were measured and quantified. Protein levels were determined by Western blot analysis. Data from independent experiments were normalized to values of the sample harvested 0 h post IR. Error bars indicate standard deviation of two biological replicates. mRNA expression was measured using qRT-PCR. β-actin was used as an internal control. Error bars indicate the minimum and the maximum value for the relative quantity of two biological replicates, each consisting of three technical replicates. Simulated dynamics of total p53 (g), Mdm2 mRNA (h), Mdm2 (i), Wip1 mRNA (j), Wip1 (k) and pATM (l) are shown for eight parameter combinations and IR alone.

Fig 7. Scheme of the two identified mechanisms allowing to reflect the effects of IKK2 inhibition on p53 dynamics.

a) The first mechanism comprises the two parameter combinations α_mpi & α_m & β_sp and α_mpi & α_m & T_s. The dotted lines indicate that the parameter combinations differ in the corresponding parameters (β_sp and T_s), but the affected process is the same. In both parameter combinations the process is predicted to be reduced, which is denoted by the blue color. In contrast, the red solid lines represent the two processes shared among the two combinations, which are predicted to be increased. b) The parameter combinations α_mpa & α_sm & β_sp and α_mpa & α_sm & T_s constitute the second mechanism. In accordance to the first mechanism, the two parameter combinations differ in parameter β_sp and T_s. The two parameters affecting ATM-induced degradation of Mdm2 (α_sm) and Mdm2-mediated degradation of activated p53 (α_mpa) are predicted to be reduced upon IKK2 inhibition.