SnC dynamics during ageing in mice

To find which of the model mechanisms best describes SnC dynamics, and with which rate constants, we compared the models to longitudinal data on SnC abundance in mice collected by Burd et al. 1. SnC abundance was measured using a luciferase reporter for the expression of p16INK4a, a biomarker for SnCs. Total body luminescence (TBL) was monitored every 8 weeks for 33 mice, from early age (8 weeks) to middle–late adulthood (80 weeks) (Fig. 2a).

Fig. 2 Saturated-removal (SR) model captures longitudinal SnC trajectories in mice. a Total body luminescence (TBL) of p16-luciferase in mice (n = 33). Gray lines connect data from the same individual mice (green and purple lines are examples of individual trajectories). b SR model equations and their approximate analytical solutions. The SR model (red line) captures c the mean SnC abundance, d standard deviation of SnC abundance, e skewness, and f shape of the distributions among equal-aged individuals, and g correlation between subsequent measurements on the same individuals. TBL was normalized to give a mean abundance of 1 at young ages. Maximum-likelihood parameters for the SR model are: η = 0.15 day−1 year−1, β = 0.27 day−1, κ = 1.1, ε = 0.14 day−1. Pink lines in c: best-fit of all models without saturation mechanism iv, that have an age-related increase in SnCs, best-fit parameters are in Supplementary Note 1. Mean and standard error (shaded red, pink regions) are from bootstrapping. Source data are provided as a Source Data file. Full size image

The luciferase in these mice was introduced into one of the p16 loci, causing the mice to be heterozygous for p16, which may impair proper activation of the senescence program. We therefore also tested longitudinal measurements of SnCs based on another method. For this we obtained longitudinal data from Yamakoshi et al. 30, who measured SnC abundance by creating a transgenic mouse model with a human p16 gene tagged with luciferase, retaining the native p16 loci. Although this dataset has much fewer mice, it shows similar dynamics to the dataset of Burd et al.1 (Supplementary Note 1, Supplementary Fig. 2), suggesting a similar underlying dynamical process.

We tested how well each model describes the longitudinal SnC trajectories of Burd et al.1 by finding the maximum-likelihood parameters for each of the 16 models, adjusting for number of parameters (Supplementary Notes 1 and 2, Supplementary Tables 1–4). A principle emerges from this comparison: in order to capture the longitudinal dynamics, the mechanism must have rapid turnover of SnCs on the timescale of a few days in young mice, and it also must include mechanism (iv), which represents a decline in removal that depends on SnC abundance rather than directly on age. The simplest model that describes the data thus has only two interactions (Fig. 1c): SnC production rate increases linearly with age (mechanism i), and SnCs slow down their own removal rate (mechanism iv). We call this model the saturating removal model (SR model), whose equation is given in Fig. 2b.

The SR model captures the accelerating rise of mean SnC abundance with age in the longitudinal data (Fig. 2c and Supplementary Figs. 3, 4): as SnCs accumulate, they slow their own removal, leading to even higher SnC levels. The SR model also explains the SnC variability between individuals which accelerates with age (Fig. 2d), and the SnC distributions among equal-aged individuals (Fig. 2e), which are skewed to the right (Fig. 2f).

Importantly, the SR model captures the fact that SnC fluctuations become more persistent with age, as evidenced by an increasing correlation between subsequent measurements (Fig. 2g, F-test for linear regression, p-value 0.0047; F-statistic 16.5): individuals with higher (or lower) than average SnC levels stay higher (or lower) for longer periods with age. This increased persistence is due to the effect of SnCs on their own removal rate. Models without mechanism iv (saturation of removal) show a poor overall fit (pink lines in Fig. 2c, ΔBIC > 44.3).

SnC lifetime is days in young mice and weeks in old mice

The maximum-likelihood parameters of the SR model (listed in the caption of Fig. 2) provide quantitative predictions for SnC half-lives: SnC turnover is rapid in young mice, with a half-life of about 5 ± 1 days at 3 months of age; Turnover slows with age, so that SnC half-life is about 25 ± 6 days at 22 months.

We tested these predictions using experiments in mice by inducing SnCs and analyzing their dynamics. To induce senescence in mice lungs we used intra-tracheal bleomycin administration (Fig. 3a), a DNA-damaging agent that induces cellular senescence in the lung epithelium a few days after treatment5,31.

Fig. 3 SnC half-life measurements in mice support SR model predictions. a Bleomycin or PBS was introduced by intratracheal installation to mice on day 0. Lungs were analyzed on the indicated days thereafter. Representative images of lung cells analyzed by imaging flow cytometry show how senescent epithelial cells were identified, using SA-β-Gal, Pan-Cytokeratin (pCK), and DAPI staining. SnC removal rate was estimated by log-linear fit. b The SR model predicts that SnCs rapidly return to baseline in young mice and that removal is slower in old mice. c Fraction of SnCs in mouse lungs after treatment with bleomycin (1.5 U/kg). In young mice, SnC levels return to baseline with a half-life of about 5 days. In old mice, baseline SnC levels are about five-fold higher, and SnC removal rate is slower than in young mice . d SnC removal rates (half-life−1) for young and old mice (mean and standard error from bootstrapping, black) agree with the SR model predictions (red, mean and SE were calculated by bootstrapping, see the “Methods” section). The best-fit model without mechanism (iv), the USR model (mechanisms i + iii), shows a poor prediction (pink). For both ages, the USR prediction is different from the observed half-life with p < 0.01 from bootstrapping. Source data are provided as a Source Data file. Full size image

We quantified the fraction of senescent lung epithelial cells at different time points following bleomycin administration (Fig. 3a) using imaging flow cytometry. Epithelial SnCs were defined as cells positive for a senescent cell marker (SA-β-Gal) and an epithelial marker (pan-Cytokeratin, pCK). This cell population was also HMGB1 nuclear negative, as expected in SnCs5,32, and previously shown5 to correspond to non-proliferative cells (negative Ki67 assay, see Supplementary Note 3, Supplementary Fig. 6).

In 3-month-old mice, SnC levels decayed with a half-life of τ = 4.7 days (τ−1 = 0.21 +/− 0.07 days−1) and reached their baseline level within less than a month (Fig. 3b, c), as predicted. SnC levels in young mice lungs are thus in a rapid dynamic balance of production and removal.

To test the prediction that removal slows with age (Fig. 3b), we performed the bleomycin treatment in old mice (22-month old). In these mice, the baseline level of SnCs was about five-fold higher than in young mice (Fig. 3d). SnCs decayed with a half-life of τ = 18 days, τ−1 = 0.055 +/− 0.035 days−1), slower than that of young mice as predicted (p = 0.038 from bootstrapping, Fig. 3b).

These turnover measurements quantitatively agreed with the predictions of the SR model (Fig. 3d, Supplementary Note 4, Supplementary Fig. 7) with no additional fit. This agreement occurred despite the use of distinct SnC markers in the two data sets (SA-β-Gal in the bleomycin experiment vs. p16INK4A-luciferase in the longitudinal experiment), suggesting consistency between the measurement methods.

Our results suggest a core mechanism in which SnC production rate rises linearly with age, and SnCs slow their own removal (Supplementary Note 5, Supplementary Fig. 8). This slowdown of removal accelerates SnC accumulation with age. Slowdown of removal also amplifies fluctuations in SnC levels at old ages. This amplification, known as critical slowing down33,34, results in long-lasting differences among individuals at old ages. In other words, young mice have large spare removal capacity of SnC; old mice have much smaller spare removal capacity. This smaller removal capacity means that any addition of SnCs takes longer to remove, causing larger and more persistent variation in SnC levels among individuals (Fig. 2g).

The SR model quantitatively recapitulates the Gompertz law

In the remainder of the paper, we use mathematical analysis to explore the implications of these findings for the question of variability in mortality. Mortality times vary even in inbred organisms raised in the same conditions, demonstrating a non-genetic component to mortality35,36. In many species, including mice and humans, risk of death rises exponentially with age, a relation known as the Gompertz law37,38,39, and decelerates at very old ages. The Gompertz law has no known explanation at the cellular level.

To connect SnC dynamics and mortality, we need to know the relationship between SnC abundance and risk of death1. The precise relationship is currently unknown. Clearly, SnC abundance is not the only cause for morbidity and mortality. It seems to be an important causal factor because removing SnCs from mice increases mean lifespan25, and adding SnCs to mice increases risk of death and causes age-related phenotypes23. We therefore explored the simple possibility that death can be modeled to occur when SnC abundance exceeds a threshold level X C , representing a collapse of an organ system or a tipping point such as sepsis (Fig. 4a). Thus, death is modeled as a first-passage time process, when SnC cross X C . We use this assumption to illustrate our approach, because it provides analytically solvable results. We also show that other dependencies between risk of death and SnC abundance, such as sigmoidal functions with various degrees of steepness, provide similar conclusions.

Fig. 4 SR model can explain the variability in mortality between individuals. a To model the relation between risk of death and SnC levels, we assumed a simple threshold model where death occurs when SnC abundance exceeds a critical threshold X C . b Mouse mortality (C57BL/6J mice obtained from the Mouse Phenome Database60, black line) is well fit by the SR model (red line) with parameters consistent with the data of Figs. 1, 2, with death defined when SnC exceed a threshold (η = 0.084 day−1 year−1, β = 0.15 day−1, κ = 0.5, ε = 0.16 day−1, X C = 17). c Similar results are obtained by assuming a more general sigmoidal dependency between SnC abundance X and risk of death: \(h = \left( {1 + {\mathrm{{e}}}^{ - \theta \left( {X - X_{\mathrm{C}}} \right)}} \right)^{ - 1}\). Parameters are the same as b, except that X C is adjusted according to the steepness parameter θ (inset). d The SR model with added age-independent extrinsic mortality of 0.4 × 10−3 year−1 (red) matches human mortality statistics61 (black). Inset: approximate analytical solution for the first passage time in the SR model shows the Gompertz law and deceleration at old ages. The parameters are similar to b, except a ~60-fold decrease in η: η = 0.00135 day−1 year−1, β = 0.15 day−1, κ = 0.5, ε = 0.142 day−1, X C = 17. e SR model describes rapid shifts in mortality when fully fed Drosophila transition to a lifespan-extending dietary intervention (LE), (inset: experimental data from Mair et al.45), with β = 1 h−1, κ = 1, ε = 1 h−1, η = 0.03 day -1 h−1 and X C = 15. LE was modeled by a decrease in η: η = 0.02 h−1 day−1 (changes in other parameters lead to similar conclusions, see Supplementary Note 8). f Lifespan of C. elegans raised at different temperatures varies by an order of magnitude, but survival curves collapse on a single curve when time is scaled by mean lifespan (inset: data from Stroustrup et al. 35). The SR model provides scaling for perturbations that affect η, but not other parameters (β = 1 h−1, κ = 1, and ε = 1 h−1, η = 0.07 h−1 day−1, X C = 20, Supplementary Note 9). Source data are provided as a Source Data file. Full size image

The SR model analytically reproduces the Gompertz law, including the observed deceleration of mortality rates at old ages (Fig. 4b–d, Supplementary Note 2, Supplementary Fig. 5, Supplementary Table 5). Notably, most models without both rapid turnover and slowdown of removal do not provide the Gompertz law (Supplementary Note 2). The deceleration of mortality rates at very old ages occurs in the model due to the increased persistence of SnC at old age. Those with high SnC have already died, whereas those with low SnC retain low SnC levels for long periods of time and avoid death. The SR model gives a good fit to the observed mouse mortality curve (Fig. 4b, c, Supplementary Note 1) using parameters that agree with the present experimental half-life measurements and longitudinal SnC data (Supplementary Note 1). Thus, turnover of days in the young and weeks in the old provides SnC variation such that individuals cross the death threshold at different times, providing the observed mortality curves.

The SR model can describe the observed increase in mean lifespan of mice in experiments in which a fraction of SnCs are continually removed (Supplementary Note 6). More generally, the SR model can address the use of drugs that eliminate SnCs, known as senolytics40. To reduce toxicity concerns, it is important to establish regimes of low dose and large inter-dose spacing41. The model provides a rational basis for scheduling senolytic drug administrations. Specifically, treatment should start at old age, and can be as infrequent as the SnC turnover time (~month in old mice) and still be effective (Supplementary Note 6).

We also adapted our results from the mouse data to study human mortality curves. In humans, mortality has a large non-heritable component42,43. A good description of human mortality data, corrected for extrinsic mortality, is provided by the same parameters as in mice, except for a 60-fold slower increase in SnC production rate with age in the human parameter set (Fig. 4d, Supplementary Note 7, Supplementary Table 6). This slower increase in SnC production rate can be due to improved DNA maintenance in humans compared to mice44. We conclude that the critical slowing-down described by the SR model provides a possible cellular mechanism for the variation in mortality between individuals.

SR-type dynamics and ageing of Drosophila and C. elegans

The generality of the SR model suggests that it might also apply to organisms where ageing may be driven by factors other than SnCs, such as Drosophila melanogaster and C. elegans, in which lifespan variation is well-studied35,45. In these organisms, the present approach can be extended by considering X as a causal factor for aging, that accumulates with age and has SR-type dynamics46, namely turnover that is much more rapid than the lifetime, increasing production and self-slowing removal. One clue for the identity of such factors may be gene-expression variations in young organisms that correlate with individual lifespan47,48,49, and the actions of genes that modulate lifespan39,50,51,52,53.

Work in C. elegans and Drosophila provides constraints to test the SR model. For example, Drosophila shows rapid switches between hazard curves when transitioned between normal and lifespan-extending diets (Fig. 4e, inset). These switches are well-described by the SR model, due to its rapid turnover property (Fig. 4e and Supplementary Note 8, Supplementary Fig. 9). The rapid turnover property entails that the level of X can adjust after a change in any of the parameters of the model. A model without rapid turnover could not explain these results.

We further tested whether the SR model can explain the survival curves of C. elegans under different life-extending genetic, environmental, and dietary perturbations35. These perturbations change mean lifespan by up to an order of magnitude. The survival curves show a remarkable feature called temporal scaling: the survival curves collapse onto approximately the same curve when age is scaled by mean lifespan (Fig. 4f insets). That is, the entire distribution of death times, including its mean and standard deviation, is determined by a single parameter, which depends on the perturbation. We find that the SR model provides the shape of the survival curves, as well as their temporal scaling feature. Temporal scaling is found in the SR model by assuming that the perturbations affect the accumulation rate η (Fig. 4f, Supplementary Note 9 and Supplementary Fig. 10A).

Temporal scaling cannot be explained by models without rapid turnover (Supplementary Fig. 10B), or by varying any other parameter except η in the SR model. Thus, we predict loss of temporal-scaling of survival curves when a perturbation affects other SR-model parameters, such as removal rate β or noise ∈ (Supplementary Fig. 10C, D). This prediction may apply to exceptional perturbations in which temporal scaling is not found, such as the eat-2 and nuo-6 mutations (Supplementary Fig. 10E, F). We conclude that the SR model of rapid turnover with critical-slowing down is a candidate explanation for the temporal scaling of survival curves in C. elegans.