Senescent cell turnover slows with age providing an explanation for the Gompertz law

A causal factor in mammalian aging is the accumulation of senescent cells (SnCs). SnCs cause chronic inflammation, and removing SnCs decelerates aging in mice. Despite their importance, turnover rates of SnCs are unknown, and their connection to aging dynamics is unclear. Here we use longitudinal SnC measurements and induction experiments to show that SnCs turn over rapidly in young mice, with a half-life of days, but slow their own removal rate to a half-life of weeks in old mice. This leads to a critical-slowing-down that generates persistent SnC fluctuations. We further demonstrate that a mathematical model, in which death occurs when fluctuating SnCs cross a threshold, quantitatively recapitulates the Gompertz law of mortality in mice and humans. The model can go beyond SnCs to explain the effects of lifespan-modulating interventions in Drosophila and C. elegans, including scaling of survival-curves and rapid effects of dietary shifts on mortality.

or requests for clarifications.
1. The "saturating removal" effect introduces the only nonlinearity in X that is considered in this model: does this go some way to explain the fit. One could imagine other nonlinearities, for example a logistic growth term (rather than a linear term) in \eta_2. Would such a model fit the data? The best-fit model approximates a Gompertzian curve for age-related mortality: is this also in part simply due to the form of the saturating removal term. Another way to put this question would be: how confident are the authors that they have explored the space of possible models? 2. In Fig. 1 it would be very helpful to label the arrows of the model schematic with their associated parameters.
3. In Fig. 2, at 80 weeks the variance in TBL is large, with many mice still exhibiting low levels. Is it by chance that the two trajectories shown in Fig. 2A do not lie close to the mean trajectory of the data? It would be helpful to plot more (all?) of these trajectories overlapping, to get a better sense of these. It seems that many mice do not show increases in TBL over their lifetimes. This needs further investigation, as it may suggest that while the SR model fits the mean SnC dynamics, trajectories for individual mice may be fit better by alternative mechanisms? 4. Lines 181-186: hard to follow. Please expand on this section as it is unclear. Also, here (and elsewhere) it can be difficult to follow the references to the SI -please make these more specific wherever possible, e.g. to a specific SI figure rather than to a section. 5. Supplementary Section 2 Eq. 3 needs further justification/clarification: should the term \beta \kappa log(\kappa + X) be \beta X log(\kappa + X)? This has implications for the following quasisteady state analysis.
6. Please clarify the meaning of the term "ExpIntegralE" in the SI.

Reviewers' comments:
Reviewer #1 (Remarks to the Author): This is a very interesting paper, with an original and valuable approach to the understanding of senescence in vivo and its contribution to ageing.
We thank the reviewer for this endorsement.
I have some concerns that I think the authors could consider before publication: 1. The experimental data used to generate the mathematical model come from a mouse that is heterozygous for p16 (Burd et al. 2013). If we assume that p16 is a good read out of senescence, then, the mice used are only monitoring half the amount of p16; also if we assume that p16 is important for the implementation of senescence, then, the mice used are partially impaired in the activation of senescence. I think that the authors should acknowledge these caveats and should address them.
2. There is another mouse model that retains normal endogenous p16, but carries a transgenic allele of human p16 fused to luciferase (the resulting human p16-luc protein is inactive regarding hp16 but active regarding luc) (Yamakoshi et al., 2009: https://www.ncbi.nlm.nih.gov/pubmed/19667129 We thank the reviewer for these comments and suggestions. We addressed both comments by providing new analysis and additional data in the revised results and SI. We obtained primary data from Yamakoshi et al, as suggested, and find that it shows similar SnC dynamics to the Burd et al dataset. We now address these points in the main text with a new paragraph: "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  (2), which included luciferase-based measurements of SnC abundance in 33 individual mice, from early age (8 weeks) to middle-late adulthood (80 weeks). The measurements were based on a knock-in allele, in which the firefly luciferase gene was targeted into one of the endogenous p16 loci.
While the luciferase output retains the cis-regulatory elements of p16, the resulting mouse is heterozygous for p16. We therefore tested longitudinal measurements of p16 based on another method.
For this, we obtained data from Yamakoshi et al. (12), which created a transgenic mouse model with a human p16 gene tagged with luciferase, maintaining the native p16 loci. The dataset contains 7 mice (3 male and 4 female), whose total body luminescence was monitored every month ( Figure S2A).
Luciferase output increased with age for these mice. Unlike in the Burd et al. dataset, there were marked differences in p16 LUC between male and female mice in this dataset. While young females in this dataset had very low luciferase output (as in the Burd dataset), young males had luciferase output that was already similar to that of old females. The dataset for the males also had many missing values, with a third of the time-points containing only a single mouse. For this reason, we chose not to aggregate the measurements and to focus our analysis on the female mice.
The Yamakoshi dataset is much smaller than the dataset of Burd et al. (4 mice compared with 33 mice), and we therefore cannot fit it with a dynamical model. However, it resembles the Burd et al.
dataset in several important quantitative aspects. The mean luciferase measurements increase at a similar rate with age to the Burd et al. dataset (10% ± 2% per month compared with 10% ± 1%), as does the standard deviation (12% ± 3% per month compared with 9% ± 1%) ( Figure S2B). In addition, as predicted by the SR model, the Yamakoshi et al. mice show more persistent variation at old ages than at young ages. This can be quantified by the mean rank, which mixes rapidly in the young (<=10mo) showing a mean rank around 2-3, but stays more persistent in the old (>10mo): the individuals with high SnC remain high and those with low SnC remain low, Figure S2C. These similarities indicate that the Burd et al and Yamkoshi et al constructs may report similar underlying dynamics.

I don´t fully grasp the difference between models (iii) and (iv). In (iii), removal slows down because the
immune system slows down, but this reduced activity of the immune system can be consequence of the paracrine actions of senescent cells, which then sounds very similar to model (iv).

Is model (iv) the same as the "saturation removal" (SR) model? Please, clarify the difference between
(iii), (iv) and SR. In lines 123 and 141, authors refer to their model as: "rapid SnC turnover and removal slowdown". I like this way of referring to their model because I understand it immediately. However, I am not sure if this is (iii), (iv) or SR, or none? I suggest that if authors have no data on "saturation" of the immune system, that they simply say "impaired immuneclearance" (by saturation, or because senescent cells make themselves invisible, or because the extracellular matrix is impenetrable, or because the bone marrow does not produce healthy active immune cells, etc).

These two comments helped us to clarify the difference between (iii) and (iv). We now added the following explanation of this point in the introduction:
"(iii) SnC removal decreases with age due to age-related decline in immune surveillance functions 30 , and (iv) SnCs reduce their own removal rate, which can be due to SnC-related signaling such as SASP, downregulation of immune surveillance by SnCs, SnCs saturating immune surveillance mechanisms (similar to saturation of an enzyme by its substrate), or to disruption of tissue and extracellular matrix architecture that interferes with removal. " "Mechanism (iv) is distinct from mechanism (iii) because the decline in removal rate in (iv) depends on SnC abundance, rather than on age directly. Although (iv) can arise from various biological processes, we denote it for simplicity 'saturation of removal'." We further clarified the meaning of the SR model by the following revised text in the results: "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 ( Figure 2B): 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 Figure 2B. "

I find unlikely that senescent cells may reach such high levels in old mice as to saturate the immune system (innate or adaptive). Is there any evidence for this?
We thank the reviewer for this question. As far as we know it is unclear what the absolute abundance of senescent cells is in old tissues. One can estimate how high it needs to be in order to saturate the immune system: the known cell types that remove SnCs include NK cells and macrophages, which together make up on the order of 0.1% of the body's cells. Thus, if senescent cells turn out to approach these levels at old ages, saturation is possible. Otherwise, saturation is unlikely and slowdown of removal may be due to other mechanisms such as SASP. We now added a paragraph in the discussion which addresses it: "The present analysis of longitudinal p16 trajectories suggests that SnC slow down their own removal rate. This effect may be due to several mechanisms, including SASP, disruption of tissue architecture, or SnC abundance exceeding immune capacity. For the latter effect, SnC abundance at old age needs to be comparable to the abundance of the immune cells that remove them, which make up on the order of 0.1% of the body's cells 59,60 . Further research is needed to characterize these effects." 6. In Fig. 2C, which are the models without saturation removal?
We now clarify this in the figure legend of Figure 2, and by adding the interactions i-iv present in each of the models: best-fit of all models without saturation mechanism iv, that have an age-related increase in SnCs, best-fit parameters are in Supplementary Section 1. Mean and standard error (shaded red, orange regions) are from bootstrapping.

How good or bad is model (iii) compared to (iv)?
To address this question, we compared the models using Bayesian information criterion (BIC) and added the results of this comparison to the main text and to the supplementary section 1. Very strong evidence is considered to be ΔBIC above 10. The ΔBIC between the SR model and the models with mechanism iii (and without iv) is at least 44.3, suggesting very strong support for mechanism iv compared to iii. We also show in Fig 2C the

best-fit curves to the Burd dataset with all models that have (iii) instead of (iv), and detail what interactions the models have next to their
curves, for easy comparison.
8. In Fig. 3D: authors show how bad it is the best fit model without SR. Again, model (iii), how bad is it?
In Fig 3D we compare the SR model (mechanism i+iv) to the equivalent model with mechanism iii (mechanism i+iii). The latter mechanism shows best-fit parameters that are not compatible with the observed turnover times (p<0.01). We added the following sentence to the legend of figure 3: The best-fit model without mechanism (iv), the USR model (mechanisms i+iii), shows a poor prediction (orange). For both ages, the USR prediction is different from the observed half-life with p<0.01). We now add this to the discussion:

I miss discussion of other papers that have looked at the rate of senescent cell elimination in
The rapid removal of senescent cells that we observe following bleomycin-induced DNA damage is in line with studies that showed efficient removal of senescent cells in-vivo following liver fibrosis or induction by of senescence by mutant Ras 55-57 . On the other hand, when senescence was induced in the skin by directly activating the cell-cycle inhibitor p14ARF, which was not associated with an increase in tissue cytokine expression or inflammation, the induced senescent cells persisted in the tissue for several weeks 58 . Clearance may thus depend on the tissue, on the method of senescence induction, and on the presence of SASP.

My particular interpretation of the Gompertz effect is that those individuals with high fluctuations (in this case senescent cells) are those that die first, so that there is a selection for those that fluctuate less and this reflects as a slow down in death risk. Would the authors agree with this? I am trying to suggest an interpretation that is easier to communicate.
We thank the reviewer for this suggestion. We agree that deceleration of aging occurs because of selection of those with low SnC. The intuitive reason involves the extended persistence of SnC due to the slowdown of removal. We now added this to the results: The deceleration of mortality rates at very old ages occurs in the model due to the increased persistence of To address this, we scanned alternative models with nonlinearities other than saturation, including the suggested logistic growth term in the autocatalysis term. These results are summarized in a new supplementary section: Additional models with nonlinear mechanisms.
In the circuit scan described so far, we showed that saturation of the SnC removal mechanism is required to explain the longitudinal trajectories of Burd et al. The saturation effect is modeled using a Michaeles-Menten term for the SnC-dependent reduction in removal rate, . This term introduces a non-linear dependence on SnC abundance, . The autocatalysis effect, on the other hand, is modeled using a linear dependence on SnC abundance. To test whether other non-linear effects may explain the longitudinal trajectories, we tested several other plausible models. We kept the removal term linear in X, and added non-linearity to the production term.
First, we tested models where there is no saturation, and the non-linearity is introduced from logistic effects. In the first model, the logistic effect is introduced in the total production of SnC, and in the second model, it is introduced specifically in the autocatalysis factor: Both models do not improve on the USR model (same equations with eta3=0). They have a log-likelihood of -500, and the best-fit values for , are 0.
We also tested a model where autocatalysis has a quadratic term : This model improves on the other USR models, and has a maximal log-likelihood of = −486 (BIC=1012). This model is much worse than the best-fit SR model (ΔBIC=34). We therefore conclude that these models are insufficient to explain the longitudinal trajectories. Fig. 1 it would be very helpful to label the arrows of the model schematic with their associated parameters.

In
We have now labeled the arrows in Fig. 1 as suggested. 3. In Fig. 2, at 80 weeks the variance in TBL is large, with many mice still exhibiting low levels. Is it by chance that the two trajectories shown in Fig. 2A  We thank the reviewer for this comment. We compared with slope distribution predicted from the parametrized SR model.
So far, we modeled individual variation in SnC accumulation as resulting from stochasticity in which we described using a white noise term √2 . Variation in senescent cell accumulation may also result from inter-individual variation in model parameters. This inter-individual variation in parameters is relevant for the case of humans, who vary in genotype and environment. There is also likely to be some variation in model parameters between the inbred mice, whose SnC levels were measured longitudinally by Burd et al ( Figure S1A). To test whether it is justified to model variation in model parameters between individuals, we estimated the persistent variation between individuals. For this, we calculated the rank of every individual at each time point, and estimated the mean rank of every individual throughout their lifetime ( Figure S1B). In the most extreme case of inter-individual variation, where individuals have separate trajectories of SnC accumulation throughout their lifetime, we expect this distribution to be uniform (every individual maintains its rank throughout its lifetime, a 'deterministic individuality' model). On the other hand, the stochastic variation described by the parametrized SR model suggests a much narrower variation in mean rank, which is similar to what is observed in the dataset.
Finally, we estimated the variation in SnC accumulation by calculating the slope of SnC accumulation for every individual ( Figure S1C). This corresponds well with the variation in SnC accumulation rates described by the parametrized SR model ( Figure S1C). We conclude that stochastic models are adequate for describing variation in SnC dynamics for the Burd et al. dataset.
We note that the comparison of models to the longitudinal data in the main text was done with the full data trajectories, and not to their average statistics.
4. Lines 181-186: hard to follow. Please expand on this section as it is unclear. Also, here (and elsewhere) it can be difficult to follow the references to the SI -please make these more specific wherever possible, e.g. to a specific SI figure rather than to a section.
We improved and expanded the text according to this suggestion. We defined temporal scaling, and added references to the appropriate supplementary figure panels.
We further tested whether the SR model can explain the survival curves of C. elegans under different life-extending genetic, environmental and diet perturbations. 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 (Figure 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 η ( Figure   4F, Supplementary Section 9 and Supplementary Figure S10A).
Temporal scaling cannot be explained by models without rapid turnover (Supplementary Figure   S10B), or by varying any other parameter except η in the SR model. Thus, we predict loss of temporalscaling of survival curves when a perturbation affects other SR-model parameters such as removal rate β or noise ϵ (Supplementary Figure S10CD). This prediction may apply to exceptional perturbations in which temporal scaling is not found, such as the eat-2 and nuo-6 mutations (Supplementary Figure   S10EF). 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.

Supplementary
Section 2 Eq. 3 needs further justification/clarification: should the term \beta \kappa log(\kappa + X) be \beta X log(\kappa + X)? This has implications for the following quasi-steady state analysis.
We have rechecked this calculation. The potential function is: