Dynamic Modelling of Pathways to Cellular Senescence Reveals Strategies for Targeted Interventions

Cellular senescence, a state of irreversible cell cycle arrest, is thought to help protect an organism from cancer, yet also contributes to ageing. The changes which occur in senescence are controlled by networks of multiple signalling and feedback pathways at the cellular level, and the interplay between these is difficult to predict and understand. To unravel the intrinsic challenges of understanding such a highly networked system, we have taken a systems biology approach to cellular senescence. We report a detailed analysis of senescence signalling via DNA damage, insulin-TOR, FoxO3a transcription factors, oxidative stress response, mitochondrial regulation and mitophagy. We show in silico and in vitro that inhibition of reactive oxygen species can prevent loss of mitochondrial membrane potential, whilst inhibition of mTOR shows a partial rescue of mitochondrial mass changes during establishment of senescence. Dual inhibition of ROS and mTOR in vitro confirmed computational model predictions that it was possible to further reduce senescence-induced mitochondrial dysfunction and DNA double-strand breaks. However, these interventions were unable to abrogate the senescence-induced mitochondrial dysfunction completely, and we identified decreased mitochondrial fission as the potential driving force for increased mitochondrial mass via prevention of mitophagy. Dynamic sensitivity analysis of the model showed the network stabilised at a new late state of cellular senescence. This was characterised by poor network sensitivity, high signalling noise, low cellular energy, high inflammation and permanent cell cycle arrest suggesting an unsatisfactory outcome for treatments aiming to delay or reverse cellular senescence at late time points. Combinatorial targeted interventions are therefore possible for intervening in the cellular pathway to senescence, but in the cases identified here, are only capable of delaying senescence onset.


Table of Contents
. In vitro data set used for estimating the model parameters. Table S2. Legend of the model variables.                           Up to 7 rounds of alternated parameter estimation and MOTA identifiability analysis were computed in order to progressively determine all the kinetic rate constant parameters. The internal round columns include the following labels: Assumed, Fixed and Locked. A parameter was termed assumed when it was assumed a priori, and therefore was not estimated. The only assumed parameters were related to the IKK-β dynamics. A parameter was termed fixed when it could be estimated and identified based on MOTA analysis within a confidence of variance lower than 25% or a correlation coefficient lower than 0.9. A parameter p was termed locked when it was fixed without being completely identifiable according to MOTA analysis. This was done only when 1) p belonged to a tuple of related parameters and each parameter in this tuple only related to the same parameter tuple (e.g. case for k8, k9), or 2) the other parameters in this tuple were not found to significantly relate with p, creating a one-way correlation (e.g. case for k12, see Figure S8). Since these correlations were local and completely confined to the tuple parameters, they did not affect the other unrelated parameters. The parameters in the tuple were only linearly affected. MOTA correlation matrices and plots are provided in Figures S4-S17, and computed from the best 30% of 20,000 fits (see Materials and Methods for more details). For each parameter, the final value, mean, standard deviation and coefficient of variance are reported. Interestingly, parameter estimation reported an extremely low value for the parameter k5, suggesting a poor ATP production by dysfunctional mitochondria. Among the three modalities of mTORC1 activation (amino acids only, amino acids + insulin, amino acids + IKK-β; parameters k6, k7 and k29, respectively), it was detected that mTORC1 activation was significantly stronger in the presence of insulin. Finally, parameter estimation also suggested that the increase in SA-β-gal (see k30, k31) was largely dependent on ROS rather than mitophagy. Table S5. List of antibodies. The antibodies used in this study are shown, with their respective catalogue numbers and the dilutions used for the respective staining (western or immunofluorescence).    In this round the kinetic rate constants k3, k4, k17, k18, k19, k20, and k37 were fixed since they were identifiable using MOTA identifiability analysis. * : Correlation Coefficient (CC) > 0.9 and Coefficient of Variation (CV) > 0.25; ** : Correlation Coefficient (CC) > 0.9, Coefficient of Variation (CV) > 0.25 and number of tuples showing this correlation (#) > 1. Format: ParameterCode < CC CV # ("Tuple of related parameters")>.

Figure S5. Correlation plots for the round 1 of parameter estimation, as detected by MOTA identifiability analysis.
Plots for the tuples of related parameters reported in the MOTA identifiability matrix in Figure S4. In this round the kinetic rate constants k1, k2, k15, k16, k21, k30, k35, and k36 were fixed since they were identifiable using MOTA identifiability analysis.

Figure S7. Correlation plots for the round 2 of parameter estimation, as detected by MOTA identifiability analysis.
Plots for the tuples of related parameters reported in the MOTA identifiability matrix in Figure S6. In this round the kinetic rate constants k25, k26, and k32 were fixed since they were identifiable using MOTA identifiability analysis. The parameters k8, k9, and k12 were locked as reported in Table S4. The parameters k8 and k9 were part of two couples of related parameters. The former k8 with k7, the latter k9 with k10. All these two couples only related internally with themselves and therefore formed two locally defined correlations. The parameter k12 only related with k11, although k11 was not dependent on k12. Hence, this also formed a locally defined correlation.  To further investigate the source of variability, principal component analysis (PCA) of the model was computed at each round of parameter estimation. Interestingly, this analysis indicated that only about half of the estimated parameters did not significantly contribute to the overall variance for each round. This suggested that a group of parameters could have been potentially identified at each round. Figure S19. Simulated tools for model inhibition or over-activation over time.
(A) To inhibit ROS over time, a simulated ROS inhibitor species was added to the model. This new species reduced ROS levels and acted as an in vitro ROS scavenger. The abundance for this species was estimated in order to achieve ROS inhibition to 10% (blue) as compared to the control (white) throughout the time course. (B) In analogy, two new species were created for over-activating mitophagy and FoxO3a, respectively. The abundance for these two species were estimated to achieve a mitophagy or FoxO3a over-activation of 150% (magenta) as compared to the corresponding control (white) throughout the time course. Each boxplot represents the median and two quartiles, whereas the bars indicate the minimum and the maximum values, as estimated from day 1 to day 21.

Figure S20. Inhibition of ROS and mTOR in vitro.
Torin and ROS inhibition efficacy. (A) Cells were irradiated with 20Gy X irradiation and then treated with Torin1 or DMSO as described in Methods. Lysates were probed with total mTORC1 and mTORC1-S2448 antibodies. Band intensities were quantified relative to Tubulin loading control and plotted as mean +/-SD ratios of mTORC1-S2448 to total mTORC1 (3 repeats). (B) Cells were irradiated with 20Gy X irradiation and then treated with SOD and catalase as described in Methods. Cells were stained with DHE to measure intracellular superoxide levels by flow cytometry. Time course data over 21 days are plotted (n=3).

Figure S21. Analysis of the two mitochondrial sub-populations upon ROS-mTOR combined intervention.
The internal states for the new and old mitochondrial sub-populations were also investigated upon combined ROS-mTOR simulated intervention. Regarding the mitochondrial membrane potential, the global effect of these interventions mainly resulted from changes in the sub-population of new mitochondria, which showed a strong synergistic response at particular levels of mTOR-ROS inhibition. The effect upon old mitochondrial ym was almost entirely dependent upon mTOR inhibition, but still changed their potential by an insignificant amount compared to the new mitochondrial population. Concerning the mitochondrial mass, the perturbation of combined ROS-mTOR acted predominantly on the young population, although these changes were largely hidden in the overall population ( Figure 4A) due to the larger proportion being comprised of the old population. The point (0, 0) indicates the control (no inhibition).

Figure S22. Analysis of senescence-associated b-galactosidase staining upon upon ROS-mTOR combined intervention.
(A) Cells were fixed at the indicated timepoints and assayed for senescence-associated b-galactosidase followed by counterstaining with nuclear fast red prior to imaging. An average of ~300 cells were counted per coverslip, with any blue cytoplasmic staining being considered positive. Data are n=2 ± SD. Combined intervention produced significantly lower positive cells at 3 and 6 days post irradiation relative to DMSO control (P < 0.05). (B) Example images of stained coverslips showing the decreased intensity of staining observed in treated cells. Figure S23. Additional readouts upon AMPK, FoxO3a or mitophagy simulated over-activation. The model predicted a decrease in the DNA damage/oxidative stress response pathways upon over-activation of AMPK, FoxO3a or mitophagy. AMPK was predicted to achieve the strongest inhibition throughout the time course as compared to FoxO3a or mitophagy. Interestingly, an over-activation of FoxO3a was predicted to initially increase and then reduce CDKN1A levels as compared to the control (black line), suggesting differential regulation of cell cycle arrest over time. Model stochastic simulations up to 20 days graphically showed increased stochasticity for the oxidative stress/DNA damage signalling species. Moderate increased variance was also detected for the species Mitophagy and the species representing mitochondrial mass and membrane potential. Number of stochastic runs: 500; black line indicates the means, dark grey area indicates 95% confidence interval of the mean and grey area indicates a standard deviation. (A) Species associated to in vitro data (observable variables). (B) Mitochondrial internal states (derived variables).