Synergistic effects of repair, resilience and retention of damage determine the conditions for replicative ageing

Accumulation of damaged proteins is a hallmark of ageing, occurring in organisms ranging from bacteria and yeast to mammalian cells. During cell division in Saccharomyces cerevisiae, damaged proteins are retained within the mother cell, resulting in an ageing mother while a new daughter cell exhibits full replicative potential. The cell-specific features determining the ageing remain elusive. It has been suggested that the replicative ageing is dependent on the ability of the cell to repair and retain pre-existing damage. To deepen the understanding of how these factors influence the life of individual cells, we developed and experimentally validated a dynamic model of damage accumulation accounting for replicative ageing on the single cell level. The model includes five essential properties: cell growth, damage formation, damage repair, cell division and cell death, represented in a theoretical framework describing the conditions allowing for replicative ageing, starvation, immortality or clonal senescence. We introduce the resilience to damage, which can be interpreted as the difference in volume between an old and a young cell. We show that the capacity to retain damage deteriorates with high age, that asymmetric division allows for retention of damage, and that there is a trade-off between retention and the resilience property. Finally, we derive the maximal degree of asymmetry as a function of resilience, proposing that asymmetric cell division is beneficial with respect to replicative ageing as it increases the lifespan of a given organism. The proposed model contributes to a deeper understanding of the ageing process in eukaryotic organisms.


S1.1 Data of cell growth for individual cells of S.cerevisiae
For the purpose of validation, the model is fitted to experimental data of cell growth (Fig S1).
Using microscropy, the increases in cell area A (µm) 2 over time t [min] have been measured for numerous both young and old wildtype cells of S.cerevisiae. As the model is qualitative, it describes the dynamics of intact and damaged proteins regarded as bulk quantities. Thereby, the model describes overall features of the cell as oppose to physical quantities, and in order to conduct the validation both the data and the model must be rendered dimensionless. For both the modelling output and the experimental data, the "area" and the time must be rendered dimensionless.
To exemplify the scaling procedure, the left hand figure displays the raw data of a growing daughter cell (Fig S1). The studied daughter is young in the sense that it is born without damage, and in the raw data it grows until cell division occurs which corresponds to the "jumps" in the data where the mother cell looses mass due to the formation of a daugther cell before continuing growing. Hence, using the cell mass before and after the jump, the size proportion parameter s can be estimated from the data by taking the quotient of the upper and lower data point, i.e. s 1 = y 1 y 2 where the subscript 1 corresponds to the first cell division to use the notation introduced in Fig S1. By taking the average of multiple "jumps", i.e. multiple cell divisions, for multiple cells it is possible to estimate the mean size proportion s which is the one that has been implemented in the model validation procedure. Remarkably enough, the jump size is rather constant which thus validates the assumption that the mother cell obtains a proportion s of the cell mass before division while the daughter obtains the corresponding proportion of (1 − s). Dimensionless time, τ = µ · t Dimensionless "Area"= P

Processed data
Cell area over time for a WT cell of S. cerevisiae Figure S1. Non-dimensionalisation of data of the increase in cell area over time for a single daughter cell of S.cerevisiae. The increase in cell area over time for a wild type cell of S.cerevisiae. The left hand figure corresponds to the raw data and the right hand figure corresponds to the processed dimensionless data. The processed data is obtained by scaling the cell area on the y-axis by the value (P div ) −1 which is estimated from the data and the time on the x-axis is rendered dimensionless by scaling it by µ = 0.5 h −1 . The value of P div is calculated by assuming that the first data point equals the value (1 − s) · P div where s is the size proportion. The size proportion s is calculated by the two "jumps" associated with each cell division.
Using the size proportion s, it is possible to estimate the value P div used to scale the area. The estimation of the parameter P div is estimated from newly budded daughter cells originating from young (in terms of RLS) mother cells. If we assume that such a mother cell retains damage, it follows that the corresponding daughter cell will be close to "damage-free" at the point of birth. Therefore the first data point, that is the initial area of these daughter cells, should be A(t = 0) = (1 − s) P div where s is the average size proportion estimated from the "jumps" associated with cell division described previously. Thus, the "dimensionless area" is obtained by scaling the data by the value Morever, the time variable is rendered dimensionless by scaling it by the growth rate µ. The choice of this scaling factor is motivated by the fact that the time variable for the models are rendered dimensionless by using this factor (see section S1.2). The resulting dimensionless time variable is defined as τ = µ · t, and the implemented value of the growth rate for the model validation is µ = 0.5 h −1 for S.cerevisiae based on the experiments by Snoep et al. [13]. The above procedure explains how the data is rendered dimensionless, however from the point of view of the experimental design it is important that the ageing phenomena is captured in the data.
As the competing damage accumulation models describe replicative ageing, it is important that the data takes ageing into account. Therefore, we have followed the cell growth of young daughter cells until the first cell division occurs and old mother cells until the last division before cell death occurs. A successful model of replicative ageing should be able to reproduce the growth curves of both a damage-free daughter cell and an old mother cell close to cell death. It is for this reason that we have chosen the following experimental design. By using the "averaging approach" (Fig S2) that will be described next, we generate two separate sets of timeseries data corresponding to an average young daughter cell and an average old mother cell. Thereafter, we say that the model that best can replicate both these two data sets is the one deemed the most realistic.
The four steps of the "averaging"-processing of the data are the following. First, we collect the data of, in this case, three damage-free daughter cells and five old mother cells close to cell death ( Fig S2A). Secondly, we extract the growth curves of a single division corresponding to the growth until the first division in the case of the daughter cells and the growth until the last division before cell death in the case of the mother cells ( Fig S2B). Thirdly, the area is non-dimensionalised by scaling it by the value P div (Fig S1) as described previously and the time is scaled so that each cell division occurs after 1 unit of time ( Fig S2C). Then, by using interpolation we average the cell area for the different cells and we multiply the the time with the average generation time for all cells in minutes before scaling the time with the growth rate µ in order to obtain the dimensionless time τ (Fig S2D). Thus, in the last step the time is "stretched back" by scaling it by both the average generation time and the growth rate µ. The results are the two time series corresponding to an average "damage-free" daughter cell and an average old mother cell close to cell death. The non-dimensionalisation of both the data and the model results in the possibility of quantifying the three modelling parameters s, P div and Q (Tab S1).  . Figure S2. Procedure for averaging the data. The left hand column shows the "damage-free" daughter cells and the right hand column shows the old mother cells. The four rows shows the "data-averaging" procedure. Subsequently, the non-dimensionalisation of the competing models, which is compatible the methodology implemented for processing of the data, is presented.

S1.2 Non-dimensionalisation of the competing models
The aim is to render the modelling output dimensionless in a similar manner to the non-dimensionalisation of the data. As the modelling output is defined asŷ(θ, t) = P (θ, t) + D(θ, t), this achieved by scaling both the time variable t and the two states P, D appropriately.
The dimensionless time τ is defined as follows.
Furthermore, the dimensionless states are the intact and damaged proteins as functions of their respective thresholds P div and D death respectively.
It is important to emphasise that since µ is used to introduce the dimensionless time variable τ , it is removed from the presented model during the non-dimensionalisation procedure. Therefore, the growth rate µ cannot be estimated using the dimensionless model and in order to compare the dimensionless model to the data it is essential that the data is rendered dimensionless in the same way. For this purpose, we have implemented the value µ = 0.5 h −1 when the time series data was rendered dimensionless (see section S1.1). Now, among the previous models there are two that explicitly focus on the bakers yeast S.cerevisiae. These models were constructed by Erjavec et al. [8] and Clegg et al. [6]. Using the above non-dimensionalisation procedure, the dimensionless versions of the Erjavec-model (Eq (S3)) and the Clegg-model (Eq (S4)) are obtained. The differences between the various models are the assumptions underlying the construction of the ODEs which involves the forces cell growth, formation and repair of damage.
In the case of the Erjavec-model (Eq (S3)), the cell growth term is not exponential and the growth rate declines with increasing amounts of damage according to the fraction P K M +P +QD . Moreover, a fraction of the intact proteins contributes to the formation of damage while the remaining part of the intact proteins is degraded. Lastly, the damage that is repaired does not contribute to the amount of intact proteins but it is degraded instead.
In the case of the Clegg-model (Eq (S4)), the cell growth term is exponential and the growth rate declines as previously with increasing amounts of damage according to the fraction P P +QD . The formation of damage term is the same as in the current model (Eq (S7)). Lastly, a fraction of the damage that is repaired contributes to the amount of intact proteins while the remaining part is degraded and the rate of repair declines with increasing amounts of damage according to the fraction P βP +QD .
Using the scaled forms of the variable (Eq (S1)) and the states (Eq (S2)), the three (i.e. the presented model, the model by Erjavec et al. [8] and the model by Clegg et al. [6]) non-dimensionalisation procedures are presented subsequently before the structural identifiability analysis is described.

S1.2.1 Non-dimensionalisation of the presented model
The aim is to render the ODE-model (Eq (S5)), with the following initial conditions (Eq (S6)) dimensionless. Start by introducing the dimensionless form of the time (Eq (S1)) and the states (Eq (S2)). This results in the fact that the rate parameters becomes scaled by the growth rate µ and the introduction of the damage resilience parameter defined as follows.
Now, substituting the time τ into Eq (S5) yields the following equations and simplifying the above equations in addition to introducing all scaled components results in the following dimensionless system of ODEs which is the desired result.
In a similar fashion, the initial conditions are scaled in the following way and substituting these expressions into Eq (S6) yields the following equations.
Now, simplifying the above results in addition to using the previous scaling arguments yields which is the desired result.

S1.2.2 Non-dimensionalisation of the Erjavec model
The continuous part of the original Erjavec model [8] is the following.
In order to compare the three competing models, the non-dimensionalisation procedure must be conducted in the same way. The reason for this is that the data must also be non-dimensionalised in the same way as the models in order to estimate the rate parameters in the ODEs from the measured data.
To this end, we start off by changing notation of the above equation to obtain the following (Eq (S9)).
Above the names of the states are relabelled: P ← P int and D ← P dam . Also the coefficients and rate parameters are relabelled as follows.
Now, the same non-dimensionalisation procedure is conducted on the dynamical system (Eq (S9)) where the states are redefined as follows and the time variable is rendered dimensionless as follows.
Introducing the above states, the model is non-dimensionalised as follows.
Now, redefining the states and the time variable as previously described in addition to redefining the parameters as follows results in the following system of ODEs (Eq (S10)) which is the desired result.

S1.2.3 Non-dimensionalisation of the Clegg model
In a similar fashion as for the Erjavec model, the aim is to render the Clegg model dimensionless [6].
The original model is formulated as follows.
Now, the parameters and states are relabelled as follows.
which results in the following model (Eq (S11)).
As previously, the states are rendered dimensionless. Substituting these dimensionless components into Cleggs model (Eq (S11)) yields the following.
Introducing the following parameters in addition to the previously defined damage resilience quotient Q yields the following system of ODEs (Eq (S12)) which is the desired result.

S1.3 Structural identifiability (SI) of the models reveals that Q cannot be estimated in the model fitting
Although Q can be estimated from the actual data points, it is of interest to see which parameters in the various models that can be estimated from a least square fitting of the models to the times series data of cell area over time. From a theoretical point of view, the formulation of the various ODEs determines which involved parameters that can be estimated from a given output and this can be systematically investigated using what is called a structural identifiability (SI) analysis.
For the purpose of conducting such an analysis, the three candidate damage accumulation models have the following structure where the functions f 1 and f 2 govern the dynamics of the system and these functions are determined by the forces cell growth, formation and repair of damage. The initial conditions are determined by the functions f P and f D (Eq (S8)) and the vector θ contains all rate parameters. The structure of the system of ODEs (i.e. the functions f 1 and f 2 ) determines which parameters in θ that can be identified from the outputŷ (θ, τ ). The parameters that can be identified from the outputŷ (θ, τ ) for a given model are called structurally identifiable and we denote these by ∼ θ. The results of the structural identifiability analysis (Tab S2) of the three competing models have been obtained using the Mathematica software developed by Karlsson et al. [10]. Table S2. The structural identifiable parameters of the three competing models. The first column lists the three models (the current, Erjavec et al. [8] and Clegg et al. [6]), the second column lists all parameters θ in each model and the third column lists the structurally identifiable parameters ∼ θ from the the outputŷ (θ, τ ).

S1.4 Model validation shows that the presented model outperforms other similar models but that the parameters are not numerically identifiable
The models are validated by fitting the SI-parameters to the experimental data. The model that best replicates the data of both the young and old average yeast cell ( Fig S2) is deemed most accurate. However, even though the optimal parameters in the selected model are structurally identifiable, it is not necessarily the case that these parameters can be identified from the data and then the actual parameters cannot be used for predictions. For the parameter estimation, start by denoting the non-dimensionalised measured data y(τ ). Further, assume that the data has been measured for m ∈ Z + distinct time points τ i , i = 1, 2, . . . , m. Then, the underlying assumption of the classical least square parameter estimation procedure is that the error e between the measured data y and the simulated outputŷ is normally distributed with zero mean and variances σ i .
Thus, the idea is to select the structurally identifiable parameters ∼ θ that minimise the magnitude of the errors e i for all time points (i.e. i = 1, 2, . . . , m) or equivalently selecting the parameters that minimise the variance σ. This is equivalent of selecting the parameters that minimise the least square cost function denoted LS ∼ θ and accordingly the parameter estimation problem can be formulated as follows.
When selecting models, it is often necessary to take the complexity of the models into account in addition to the fit. An example of such a criteria for selecting models is the Akaike Information Criteria (AIC) and this penalises a large number of parameters [9]. It is defined as follows where L is the likelihood function and p is the number of parameters. Under the normality assumption (Eq (S14)), the loglikelihood can be written as follows.
In our case, we have conducted a classic least square estimation using the built in function "lsqnonlin" in Matlab and in this case only the least square (Eq (S15)) is obtained. Since the error variances σ i (Eq (S14)) often are correlated it can be hard to reconstruct the loglikelihood function (Eq (S17)) from the least square LS. To this end we have assumed a constant variance s 2 estimated for the obtained errors e i as follows.
In other words, we assume no correlation between the measurment errors e i (Eq (S14)) which entails that we assume that s = σ i ∀i ∈ {1, . . . , m}. Thus, we have approximated the loglikelihood (Eq (S17)) as follows and our AIC-values (Eq (S16)) for the three models are calculated using the above approximation (Eq (S19)). Using the "least square"-methodology (Eq (S15)), the three rivalling models have been validated (Tab S3, Fig S3). Two major conclusions can be drawn. Firstly, the presented model (Eq (S7)) has a substantially better fit compared to the rivalling models. Therefore, the structure of the model is therefore the most biologically realistic. Secondly, the estimated parameters are not numerically identifiable (NI). The numerical identifiability of the parameters are determined by two factors, namely the relative size of the parameter variance compared to the parameter value itself and the correlation between the various parameters [7]. An unbiased estimator of the parameter variance is called the coefficient of variation defined as the quotient of the parameter variance and the parameter value itself, i.e. c = σ k k . In order for the parameter k to be numerically identifiable, the value of c should be close to zero percent. For all three models the coefficients of variation are very large (Tab S3) and therefore the parameter values are not numerically identifiable. In addition, the values of the correlations between the parameters should be close to 0 as well in order for the parameters to be numerically identifiable, which is not the case for any of the three models either. For the non-dimensionalisation of the time τ , the growth rate µ = 0.5 h −1 was implemented [13]. Table S3. Model validation and numerically identifiability. The five columns are from left to right: The model for which the validation has been conducted, the best fit in terms of the least square (LS), the Akaike Information Criteria (AIC), the optimal parameters obtained from the estimation and the coefficient of variation for each optimal parameter in percent. In addition to the high coefficient of variation, it is also important that the correlation between the parameters are close to zero [7]. The correlations between the parameters for the presented model are the following.

Model Best
The correlations between the parameters for the model by Erjavec et al. [8] are the following. 1 1.00 1. The correlations between the parameters for the model by Clegg et al. [8] are the following.
We next examined the ability of the models to account for replicative ageing by examining if they satisfy the criteria that RLS ∝ Q for a fixed set of rate parameters. Thus, by picking parameters giving rise to a finite RLS for all three models, an increase in Q while keeping the remaining parameters fixed should increase the life span, while a decrease in Q should decrease the life span. The simulations generated by the models reveal that the presented model together with the model by Erjavec et al. [8] satisfy this criteria while the model by Clegg et al. [6] does not (Fig S4).  In summary, the following conclusions can be drawn. The current model outperforms the rivalling counterparts in terms of the fit to the experimental data (Tab S3) and it can moreover describe ageing in yeast ( Fig S4). However, all of the estimated parameters in all three models are NI indicating that none of the values of the optimal parameters are precise. This implies that there are infinitely many sets of parameters that can give rise to the same fit for all three models. Despite the issues with the identifiability of the parameters, the good fit of the presented model entails that at least parts of the underlying assumptions of the model is true in some sense. Therefore, it is of interest to analyse the mathematical properties of the model at hand.

S2 Mathematical analysis
Mathematical analysis of the damage accumulation model is conducted to verify its validity. Previously, similar models [1,8,12,3,4,6,11] have answered research questions by means of simulation. However, there are mainly two main problems associated with the implemented methodology in these articles. Firstly, in the construction of these models there has been a lack of motivation as to why a particular model is constructed as it is and whether or not the choice of model parameters is biologically realistic. Secondly, attempts to estimate values of model components and parameters that are very hard or even impossible to measure have been done by searching through the literature. Both these problems can be addressed by investigating the properties of the model at hand by using mathematical analysis. Consequently, the presented damage accumulation model is much more complete as it addresses these two problems by using mathematical analysis as a method in addition to simulations. As the model is divided into two parts, the discrete (Eq (S8)) and the continuous (Eq (S7)) so is the analysis of the model. The first section concerns the analysis of the discrete part of the models where mathematical relations between the parameters s, re and Q are the focus. The second section concerns mainly the parameters g, k 1 and k 2 controlling the dynamics of the continuous part of the model. Here, various constraints on the damage formation rate k 1 and the damage repair rate k 2 are derived in order to elucidate the allowed sets of parameter that give rise to reasonable dynamics in terms of RLS.

S2.1.1 The relation between retention, resilience and asymmetry
We wish to derive the following two analytical results under the following two assumptions.
1. P 0 given by f P (Eq (S8)) is greater or equal to P 0,min = (1 − s) For the first result, we know that any mother cell obtains the following proportion of intact proteins after cell division.
Since re ≥ 0 it is the upper bound that must be proven (Eq (S20)). Now, starting from the first assumption listed above yields the following calculations which is the desired result (Eq (S20)).
For the second result, we use the second assumption entailing that the maximum amount of damage that the mother cell can retain corresponds to 100%, i.e. re max = 1. As re ∝ s this implies that the following holds.
By combining these two findings, the following calculations result in which yields the desired result (Eq (S21)).

S2.1.2 The increase in generation time as a function of damage resilience
In the article, the effect of the damage resilience parameter on the RLS of individual cells is presented ( Fig 1B in the article). The number of divisions is plotted against the time for low (Q = 2.6), medium (Q = 2.8) and high (Q = 3.0) resilence to damage. The length of the "stairs" in these plots corresponds to the generation time for that division and the generation time increases with the number of divisions. The effect of an increse in resilience to damage on the generation time is shown in the table below (Tab S4).

S2.2.1 Theoretical framework for the mathematical analysis
In order to pick the parameters (k 1 , k 2 ), two biologically motivated requirements are imposed.
1. Given enough food in the system, the cells should grow and as a consequence of cell growth damage should be accumulated within the cell 2. Every cell within the population has a finite replicative life span and should therefore die after a finite number of cell divisions These biological requirements are now converted into a mathematical description in the following way. The first condition is ensured by using stability analysis and the second by calculating the so called lethal initial damage threshold, D T which determines whether or not a cell will undergo cell death. The idea behind these two methodologies are subsequently presented in chronological order.
The first condition implies that both states P (τ ) and D(τ ) should be increasing functions of time which entails calculating the steady states of the system. The system, Eq (S7), has two steady states, namely one trivial steady state x 1 = (P 1 , D 1 ) = (0, 0) at the origin, and one non-trivial steady state at x 2 = (P 2 , D 2 ) = k 2 Q k 1 g, g . Thus, the first requirement can be mathematically translated into two parts. Firstly, it is essential to choose the parameters so that the non-trivial steady lies outside the square with side length 1 in the phase portrait, that is the point x 2 should be located outside the square S (Eq (S22)).
Secondly, the parameters (k 1 , k 2 ) should be chosen based on the stability of the fixed points. More precisely, the parameters should be chosen such that the trivial steady state x 1 is unstable and the non-trivial steady state x 2 is stable. As a consequence, this choice of parameters results in solutions to the system in Eq (S7) that moves away from the trivial steady state x 1 and towards the non-trivial steady state x 2 . In order to illustrate the mentioned dynamics, the proportion of intact proteins P can be plotted against the damage D in a so called phase portrait (Fig S5).
Cell death

Cell Division
Intact proteins, P Damage, D For the second requirement it is necessary to derive the so called lethal initial damage threshold D T . In fact, the mentioned threshold corresponds to a particular trajectory in the phase portrait which determines whether a cell with certain rate parameters will undergo cell death or not. More specifically, if a cell has less damage then D T at any given time it will undergo at least one more division. Conversely, if the cell in question contains more damage than D T it will undergo cell death.
In order to select parameters that ensure a finite RLS, the values of D T for both the daughter and mother cells are of interest. For parameter values such that the non-trivial fixed point x 2 (Fig S5) is located far away from the cell death side D = 1 of the square S (Eq (S22)) this trajectory will be approximately linear and can in this case be approximated using the methodology introduced by Rashidi et al. [12]. More precisely, for these values the trajectory corresponding to the threshold D T can be obtained by approximating it with numerous lines in the phase portrait. However, for small values of g where g ≈ 1, the trajectories will be very curved as can be seen in Fig S5. We still do not know how to know how to calculate this trajectory analytically and we ensure the second requirement on the RLS by means of simulations. In other words, given the rate parameters of the model at hand we simulate the RLS by solving the ODEs numerically.
Using the described approach, the aim of the mathematical analysis is to find biologically relevant parameters. Subsquently, the stability analysis is described.

S2.2.2 Stability analysis of the damage accumulation model
By using stability analysis it is possible to predict the long term dynamical behaviour of a given system. For non-linear systems of ODEs such as the damage accumulation model in this article (Eq (S7)) this goal is achieved in two steps. Firstly, the system is linearised around each steady state and secondly the stability of these linearised systems is determined. For the purpose of predicting the long term dynamics in order to ensure that the first previously posed biological requirement is satisfied, the stability of the formulated damage accumulation model is analysed.
The stability analysis is divided into four steps.
1. Find the two (i.e. the trivial and non-trivial) steady states of the damage accumulation model 2. Derive algebraic expressions for the Jacobian matrix, its determinant and trace 3. Find constraints on the model parameters such that the non-trivial steady state is stable

Find constraints on the model parameters such that the trivial steady state is unstable
The stability analysis begins with finding the steady states to the damage accumulation model (Eq (S7)). Using the notation (P , D ), the steady states of interest satisfy the following equations.
Now, Eq (S26) entails that there are two steady states. Combining Eq (S26) and (S25) yields the following two equations Thus, the damage accumulation model, Eq (S7), has two steady states: a trivial steady state x 1 (Eq (S27)) at the origin and a non-trivial steady state x 2 (Eq (S28)).
Next, a general expression for the Jacobian matrix is derived. Using the definition of the Jacobian matrix on the damage accumulation model (Eq (S7)) yields the following result.
Calculating the trace and determinant of the Jacobian matrix (Eq (S29)) results in the following two equations.
Using the formulas for the trace and determinant (Eq (S30) and (S31)) it is possible to determine the stability of the trivial (Eq (S27)) and non-trivial (Eq (S28)) steady state respectively.
Next, we would like to find parameters such that the non-trivial steady state is stable. By plugging in the non-trivial steady state, Eq (S28), into Eq (S30) and (S31) in order to obtain expressions for the trace and determinant yields the following.
According to the definition of stability, the non-trivial steady state (Eq (S28)) is stable for positive rate parameters, i.e. k 1 , k 2 > 0.
Lastly, we would like to find parameters such that the trivial steady state is unstable. Again plugging in the trivial steady state, Eq (S27), into Eq (S30) and (S31) the following expressions are obtained.
Given these expression, the aim is to find parameters such that trajectories move away from the trivial steady state x 1 . Thus, in order to ensure instability we could choose parameters such that the trace of the Jacobian at the trivial steady state is positive. This is exactly what we call the starvation constraint in the triangle of ageing which is formulated as follows.
The reason why it is called the ''starvation constraint" is due to the fact that it connects the uptake of food, i.e. g, to the rate of formation of damage and repair, i.e. k 1 and k 2 . For these parameters the trivial steady state will definitely be unstable.
However, potential problems might occur since the determinant Det [J (x 1 )] is always negative. The negativity of the term is caused by the fact that the uptake of food is positive, i.e. g > 0, as a consequence of Eq (S32) as well as the repair rate, i.e. k 2 > 0. Consequently, for parameters satisfying the starvation constraint the trace will be positive, Tr [J (x 1 )] > 0, and the determinant will be negative, Det [J (x 1 )] < 0. These parameters will therefore give rise to one positive and one negative eigenvalue and consequently the trivial steady state x 1 is characterised as a saddle point. This implies that solutions will be drawn towards the trivial steady state along the eigenvector corresponding to the negative eigenvalue and move away from the trivial steady state along the other eigenvector corresponding to the positive eigenvalue. Since the damage accumulation model is restricted to the square S (Eq (S22)) it is important that the eigenvector corresponding to the negative eigenvalue does not influence the dynamics within S. Therefore, it is necessary to ensure that the ''stable" eigenvector is not located in either the first or third quadrant.
More precisely, assume that the negative eigenvalue of J ( Then the task at hand can be formulated as follows. In order to ensure that trajectories move away from x 1 (Eq (S27)) and towards x 2 (Eq (S28)) within the square S (Eq (S22)) the sign of the first, i.e. v 1 , and second, i.e. v 2 , element of the eigenvector v must have opposite signs. Since it is always possible to scale the eigenvector v, we can set v 1 = 1 without loss of generality. Provided this notation, we need to prove that the other element of v is always negative, i.e. v 2 < 0. By definition, the eigenvector v satisfies the matrix equation which can be written as follows.
The above matrix equation is equivalent to the following two equations.
Substituting, the value v 1 = 1 into the upper of the two equations and solving for v 2 yields the following result.
The eigenvector of interest can thus be written as a scalar multiple of the following vector.
As a result of the starvation constraint (Eq (S32)) the term (g − k 1 ) must be positive. Since the eigenvalue is negative it follows that the numerator in the element v 2 is negative and the denominator is positive. Accordingly the second element of the eigenvector associated with the negative eigenvalue λ is negative, i.e. v 2 < 0, and thus the eigenvector v (Eq (S33)) lies in the fourth quadrant. Thus, for positive parameters which satisfies the starvation constraint the trajectories will move away from the trivial steady state x 1 and towards the non-trivial steady state x 2 within the square S (Eq (S22)).

S3 Numerical implementations
The numerical implementation underlying the results from the simulations yields four key outcomes. Firstly, the RLS of an individual cell is simulated. Secondly, the ageing region of the parameter space is generated. Thirdly, the comparison between the life prolonging strategies is conducted. Fourthly, the selection of parameters by fitting the model to data (see subsection S1.1) is conducted. The scripts which generated the results have been written in Matlab and are available at request. However, in this section the code is not attached but instead the ideas behind the numerical implementations are presented in terms of pseudocode algorithms.

S3.1 Simulate the RLS of a single cell
Imagine that we have a differential equation of the following type.
For stiff problems for which we have a large eigenvalue in terms of its absolute value, the dynamics initially changes quickly and later it stabilises. In these cases it is appropriate to use an adaptive method such as Runge-Kutta Fehlberg 45, RK45 method [5]. These types of methods are based on a Taylor expansion of the function f where we step forward in time given a discretisation of the time line. For example a first order approximation with a step size ∆t proceeds as follows where t k+1 = t k + ∆t and ∆t is rather small in order to get a stable algorithm for stiff problems.
In order to get a faster solution algorithm it is possible to derive a fourth x 4 and a fifth x 5 order approximation where the error is given by e = x 5 − x 4 . Subsequently the step size is adapted as long as the error e is small enough, and this is the idea behind the RK45 method. In Matlab the solver "ode45" uses this algorithm which given a time interval [t 0 , t max ] and initial conditions x 0 = x(t 0 ) as inputs the solver returns a vector with time points t and corresponding solutions x at these time points as outputs.
For our model, we would provide the maximum time τ max and the initial conditions (P 0 , D 0 ) as inputs and then the solver would return the three vectors τ , P & D. Assume that the number of time steps to reach τ max which is the length of these vectors is m. Then provided that the τ max is large enough there should be an index i ∈ {1, . . . , m} for which the i th element of the vector P , denote this P i , satisfies the following meaning that cell division occurs between these entries in time. Then using linear interpolation we can calculate the step length h = 1 − P i P i+1 − P i and then we would get the generation timeτ , the intact proteins at division P and the damage at division D according to the following assignments.
By plugging in D, the size proportion s and the retention value re into the functions f P and f D (Eq (S8)) the initial conditions for the next generation is obtained. This process can be continued until D reaches the value 1 before P does implying that cell death occurs. We present a flowchart of the above scheme below (Alg 1). , re, K S , Q, g k 1 , k 2 , µ(S) ← S S + K S 3: Initialise:

Algorithm 1 Simulation of the RLS for a single cell
Initiate the RLS, i.e. RLS ← 0 4: Solve ODEs: Solve the system of ODEs until τ max is reached with the initial conditions (P 0 , D 0 ) 5: if P ≥ 1 for some time τ i then 6: Interpolate P to 1, τ i to the generation timeτ and D to its corresponding value 7: end if 8: if D ≥ 1 for some time τ i and P < 1∀τ i then 9: Clonal senescence has occurred=⇒Terminate 10: end if 11: Define a help index, i ← 1, and the maximum number of iterations, maxIter ← 600 12: while D < 1 & i < maxIter do 13: Increment the number of generations: RLS ← RLS + 1

14:
Save the generation timeτ in the generation time vector τ , that is τ (i) ←τ

S3.2 Generating the replicative ageing region of the parameter space
The ageing are of the parameter space is generated by simulating the RLS of individual cells (Alg 1) for different parameter values. Given the various model parameters it is possible to create two vectors with equidistant values of the rates k 1 and k 2 where the number of elements is predetermined. Then by looping through these vectors it is possible to save the parameter pairs (k 1 , k 2 ) that satisfies all constraints with their corresponding RLS. The algorithm for generating the ageing region of the parameter space is presented below (Alg 2).
(P 0 , D 0 ) = (y 1 (τ = 0), 0) = (0.3680, 0.00) were y 1 (τ = 0) is the first data point of the time series corresponding to the dimensionless area of the average daughter. Above, the initial conditions build on the assumption that the daughter cell is born with no damage.
For the time series correponding to the old mother cell, we assume that the mother cell has accumulated damage. If we denote the first data point of the second time series by y 2 (τ = 0), according to the initial conditions, the damage before division is the following. D = 1 Q y 2 (τ = 0) s − 1 = 0.3012 Using the above notation, the initial conditions (Eq (S8)) are the set to following. With the above initial conditions for the two respective time series, the results from the model validation part have been obtained (Tab S3). From the solver lsqnonlin, it is also possible to conduct the identifiability analysis. This procedure relies on the Jacobian matrix which contains the time series of the derivatives of the model output with respect to the model parameters. In our case, we have the following output as an approximation of the "dimensionless area" y(θ, τ ) = P (θ, τ ) + Q D(θ, τ ) and in the case of the presented model (Eq (S7)) the following parameters are estimated from the data. where the columns consist of the so called sensitivities [2]. The sensitivities are crucial for the local algorithms that most numerical optimisers use in their iterative search for a local optima, and by using these it is possible to conduct the numerical identifiability analysis. Assume that the estimated variance s 2 of the measured data is given by where m is the number of data points and y is the average dimensionless area. In this case, we use the Jacobian matrix (Eq (S34)) to calculate the Fisher Information matrix as follows [7].
The inverse of the Fisher information matrix, F −1 , constitutes a lower bound of the covariance matrix COV(θ) and hence we approximate the covariance matrix accordingly.

COV(θ) ≈ F −1
In Matlab, it is possible to obtain the correlation matrix by using the function corrcov on the above matrix. The coefficient of variation for each parameter is obtained by dividing the corresponding diagonal element of the covariance matrix with the actual parameter values. For further reading see Chapter 12 in the book by DiStefano [7].