Mathematical modelling of telomere length dynamics

Telomeres are repetitive DNA sequences located at the ends of chromosomes. During cell division, an incomplete copy of each chromosome’s DNA is made, causing telomeres to shorten on successive generations. When a threshold length is reached replication ceases and the cell becomes ‘senescent’. In this paper, we consider populations of telomeres and, from discrete models, we derive partial differential equations which describe how the distribution of telomere lengths evolves over many generations. We initially consider a population of cells each containing just a single telomere. We use continuum models to compare the effects of various mechanisms of telomere shortening and rates of cell division during normal ageing. For example, the rate (or probability) of cell replication may be fixed or it may decrease as the telomeres shorten. Furthermore, the length of telomere lost on each replication may be constant, or may decrease as the telomeres shorten. Where possible, explicit solutions for the evolution of the distribution of telomere lengths are presented. In other cases, expressions for the mean of the distribution are derived. We extend the models to describe cell populations in which each cell contains a distinct subpopulation of chromosomes. As for the simpler models, constant telomere shortening leads to a linear reduction in telomere length over time, whereas length-dependent shortening results in initially rapid telomere length reduction, slowing at later times. Our analysis also reveals that constant telomere loss leads to a Gaussian (normal) distribution of telomere lengths, whereas length-dependent loss leads to a log-normal distribution. We show that stochastic models, which include a replication probability, also lead to telomere length distributions which are skewed.


Introduction
Repetitive DNA sequences at the end of chromosomes-known as telomeres-are shortened when cells divide, leading to one aspect of cellular aging. In this paper we derive and analyse mathematical models which describe the evolution of the distribution of telomere lengths over many generations of the cell-cycle.
In the 1960s Hayflick and Moorehead (1961) performed a series of experiments which overturned the prevailing view that normal cells were immortal. In experiments involving normal human fibroblasts cells, they found that after a finite number of divisions cell numbers reached a finite size, which is now called the Hayflick limit. Once this limit is reached, the cells become senescent: that is, they stop replicating but remain functional (Cristofalo et al. 2004). In later work, Hayflick and Moorehead observed that, the number of senescent cells in the mouse lens epithelium increases with age. Processes that contribute to senescence include oxidative stress, as shown by Muller et al. (2007) and Wei and Lee (2002); mitochondrial dysfunction (Passos and von Zglinicki 2005); somatic mutation (Kirkwood and Proctor 2003); and, of particular relevance here, telomere shortening, which was demonstrated by Allsopp and Harley (1995).
The main roles of telomeres are to protect the chromosomes against the loss of genetic material and to prevent fragments of chromosomes from rejoining (Cooper and Hausman 2009). Kirkwood (2011) has given compelling evidence that telomeres play an important role in ageing, with telomere length being a key factor in determining a cell's replicative potential. When a cell divides, its chromosomes are duplicated via a process called DNA replication. During replication, one of the daughter chromosomes is shortened at the 5' end due to the unidirectional synthesis of the new DNA chain (Olovnikov 1973). Figure 1 illustrates the effect of this process, where m, n describe the lengths of telomeres at each end of the chromosome. This process continues until the telomere length falls below a critical level and then the cell becomes senescent. We use the terms 'telomere loss' and 'telomere-shortening' interchangeably, to refer to the reduced length of telomeres in daughter cells when compared to the parent cell.
In healthy human cells, telomeres are typically 3 to 15 kilobasepairs (kbp) in length. They shorten at rates of 50-200 basepairs per replication, and undergo 30-60 popula-

5'
Daughter Chromosome 2 Fig. 1 Illustration of the chromosome replication process: on the left, the two strands of parent DNA are shown, with telomeres of lengths m, m, n and n − y; the two offspring chromosomes are shown to the right of the arrow. In each case the strand inherited from the parent is shown with a thick line, and the thin lines show the newly synthesised strand, with the arrow indicating the direction of synthesis. Daughter chromosome 2 is seen to be identical to the parent; having inherited the longer strand from the parent, and synthesised the shorter strand. However, daughter chromosome 1 inherits the shorter strand of parent DNA, and synthesises an even shorter complementary strand, resulting in telomeres of lengths m, m − y, n − y and n − y tion doublings before becoming senescent (Harley et al. 1990). When the chromosomes become too short, that is, their length falls below a threshold value, telomeres lose their protective function, triggering DNA damage which can lead to end-to-end fusions, chromosome breakage, rejoins and senescence. Evidence that telomere shortening is directly related to cell senescence can be traced back to experiments by Allsopp and Harley (1995) in which the telomeres of senescent human fibroblasts were found to be shorter than those of replicating cells and the proportion of replicating cells was proportional to the mean telomere length. Under certain conditions, the amount of telomere loss per replication is large, and cells become senescent more rapidly. For example, Wyllie et al. (2000) have shown that excessive telomere shortening associated with Werner's syndrome causes patients to experience accelerated aging. In other cases, where the enzyme telomerase is active, as in cancer cells, telomere length can be maintained or even extended, enabling the cells to become immortal (Greider and Blackburn 1985). von Zglinicki et al. (2000) suggest that there are multiple mechanisms which contribute to telomere shortening, including oxidative stress. By dosing cells with hydrogen peroxide, they showed that oxidative stress accelerated telomere loss and led to shorter replicative lifespans. The effect of perceived life stress and stress due to being a long-term care-giver on telomere length has been analysed by Epel et al. (2004), who showed that increased stress correlates with reduced telomere length. We conclude that normal ageing can be characterised by telomere shortening and that certain environmental and genetic factors may accelerate telomere shortening. Werner et al. (2015) have fitted tissue-level models to data from subjects from birth to age 85 years. In their models of cellular differentiation and telomere shortening, cells are assumed to lose a fixed length of telomere per replication; and stem cells, divide either asymmetrically (giving rise to another stem cell and a cell with shorter telomeres) or symmetrically (where both daughter cells have shorter telomeres. Fitting their model to experimental data, they obtain good agreement with the observed decreasing rate of telomere loss over time.

Simulations and mathematical models of telomere loss
Telomeres are thought to adopt the G-quadruplex structure at the 3'-end, which affects the end replication of DNA. This alters the susceptibility to the alternative lengthening of telomere (ALT) mechanism (Tan et al. 2015). Recent theoretical work of Bodova et al. (2012), Kollar et al. (2014) and Hirt et al. (2014) have focused on understanding the detailed dynamics associated with these structures and processes.
Various models of telomere shortening in a population of independently replicating chromosomes have been proposed previously: Levy et al. (1992) developed a deterministic model for telomere shortening of individual chromosomes. Their model shows good agreement with experimental data, predicting that the average telomere length decreases linearly over time. In separate work, Arino et al. (1995) assumed constant telomere loss per replication and viewed cell proliferation as a branching process. A convincing fit of their model to independent experimental data is provided. Olofsson and Kimmel (1999) extended Arino's model to account for cell death, assuming a fixed probability of cell death for senescent cells.
While many mechanisms are known to cause telomere shortening, the rate of shortening remains unclear. Many mathematical models assume a constant rate of telomere loss, for example, Levy et al. (1992). In their model, however, Buijs et al. (2004) assumed that telomere loss depends linearly on telomere length: a constant loss is attributed to the end-replication problem and an additional term is attributed to a shortening factor which they estimate from experimental data. Portugal et al. (2008) took a different approach, developing a stochastic model in which telomere shortening occurs at a constant rate but the probability of cell division depends linearly on the telomere length. The resulting stochastic model is well-described at the population level by a Gompertzian growth law.

Outline
This paper is structured as follows: in Sect. 1.1 we summarise previous theoretical work on telomere shortening. In Sect. 1.3 we outline our modelling approach. This is based on deriving equations for the distribution of telomere lengths in a population of chromosomes or cells which are repeatedly undergoing replication, causing their telomeres to shorten over time. In Sects. 2, 3, 4 and 5 we investigate four such models: in the first two, we consider populations of chromosomes that replicate independently of each other. In Sects. 4 and 5, the earlier models are scaled up to the cell level: we consider populations of cells, where each cell comprises N = 46 chromosomes. This model extension enables us to investigate how the lengths of subpopulations of telomeres evolve in cell populations. In all the models we develop, cells are assumed to exist in either a replicative or a senescent state. Once senescent, a cell remains senescent forever; we do not model cell death. Finally, in Sect. 6 we summarise the results and draw conclusions.
In each of Sects. 2, 3, 4 and 5, we start with a discrete model which describes how the numbers of cells with a particular telomere length changes from one generation, g, to the next. Since we are interested in the evolution over many generations, where all numbers will be large and, in general, the amount of telomere lost per replication is significantly smaller than the telomere length, it is appropriate to replace all these quantities with continuous variables. By considering the evolution of slowly-varying solutions of the discrete equation, we derive continuum models, namely partial differential equations (pdes) which have the same dynamics. The advantage of this approach is that, in general, the resulting pdes are more amenable to theoretical analysis than discrete systems. We use asymptotic techniques to construct and analyse solutions of the resulting pdes.
Our approach unifies and extends existing models: in particular, our approach includes, deterministic approaches where telomere loss is constant (Levy et al. 1992), dependent on telomere length (Buijs et al. 2004), as well as stochastic models, in which there is a probability of cell-replication which depends on telomere length, as considered by Portugal et al. (2008). Our aim is to explain the way in which telomere shortening occurs leads to senescence, for example, to predict how the fraction of a population which is senescent changes over time, and to determine whether it is possible to deduce from such data the mechanisms that regulate telomere shortening.
We present results which show the approximate solutions of the pdes are in good agreement with numerical simulations from the earlier paper by Qi et al. (2014).

Preliminary concepts of model development
Our main aim in this paper is to analyse mathematical models which describe the distribution of telomere lengths and their evolution over many generations. Our models are based on the earlier work of Levy et al. (1992), Buijs et al. (2004), andPortugal et al. (2008).
Initially we view the generation number as a discrete variable which we denote by g. We denote the telomere length of a chromosome by a single variable, n measured as a number of base pairs. Typically, telomeres are initially long (∼ 6000 bps) and lose only a small number of base pairs during each replication event (∼ 50 bps). Denoting by K (g) n a chromosome at generation g with telomere length n, and the number of basepairs lost upon replication by y(n), the process of replication can be described by (1.1) In Eq. (1.1), we assume that replication of a chromosome with telomere length n at generation g, produces two 'daughter' chromosomes at generation g + 1, one with telomere length n and one with length n − y(n). If the length of the daughter chromosome n − y(n) falls below a prescribed, critical value, then replication will not happen. The value of n c will be specified later.
In our analysis, we follow Buijs et al. (2004) and assume that y(n) = y 0 + ny 1 , (1.2) with y 0 , y 1 nonnegative constants, so that the number of base pairs lost is positive and may depend on telomere length. The case of Levy et al. (1992), where telomere loss is constant, is recovered by setting y 1 = 0. In general we take y 1 > 0 so that longer telomeres shorten at a faster rate than shorter ones; however, the analysis presented below is valid for the case where longer telomeres have a lower rate of loss (y 1 < 0) provided no lengthening of telomeres occurs (y(n) ≥ 0, that is, n < −y 0 /y 1 = y 0 /|y 1 |).
We term the deterministic model given by (1.1) as Case A and analyse it in Sect. 2 below. In Sect. 3 we generalise it to a stochastic model, termed Case B, by assuming that replication is a stochastic event which occurs with probability P div (n) so that the telomere replication rate depends on its length, n; that is, (1.1) occurs with probability P div (n), and otherwise, we simply have K We define the division probability as P div (n) = an + b, where a, b are constants, chosen such that 0 ≤ P div ≤ 1 for telomere lengths, n, in the range of interest. Hence for Case B, the replication rule (1.1) is replaced by Table 1 Summary of the functional forms used to model cell division (P div (n)) and telomere shortening (y(n)) where n ≥ 0 represents the telomere length of a particular chromosome and y 0 , y 1 , a, b and α are non-negative constants Model Replication rule Shortening rule References Case A P div = 1 y(n) = y 0 + y 1 n Buijs et al. (2004) Levy et al. (1992 with probability P div (n), The deterministic Case A can be viewed as a special case of Case B, for which P div (n) = 1. In the following sections, we describe and analyse the models summarised in Table 1. In Sect. 2 we start with Case A, which describes a population of chromosomes replicating independently.

Formulation of discrete model
In this section we use the replication rule (1.1) to derive a discrete dynamical system which predicts how the distribution of telomere lengths evolves over the generations. More precisely, we propose a discrete model which relates the number of telomeres of length n (n ≥ 0) at generation (g + 1) to the number of telomeres lengths n and n (where n − y( n) = n at generation g. In general, the amount of telomere lost depends only on the current telomere length. We generate a continuum limit from the discrete system by assuming that the distribution varies slowly with respect to both generation number g and telomere length n. Finally, we construct solutions of the continuum model and discuss the properties of the solutions. We start by assuming that replication occurs with probability P div = 1 and that the number of basepairs lost during replication depends on telomere length via (1.2). The replication rule (1.1) can then be written as n−(y 0 +y 1 n) , n > n c := y 0 1 − y 1 . (2.1) Thus chromosomes of length n at generation (g + 1) arise from parent telomeres with the same length, n, or a slightly longer telomere, of length (n + y 0 )/(1 − y 1 ) > n. We remark that if n < n c , then (2.1) would produce a daughter chromosome with a negative telomere length. Since negative telomere lengths are not physically realistic, Telomere length Rescaled measure of time we do not allow such events to occur. In particular, if 0 ≤ n ≤ n c then we assume that replication does not occur, and in place of (2.1) we have K so that the existing chromosome remains in the population, but does not undergo replication. Such chromosomes are termed 'senescent'. Identical results could be obtained by allowing (2.1) to occur for all n, but only consider n ≥ 0, when analysing the solution for K (g) n . The parameters and their ranges are summarised in Table 2. We now introduce K (g) n to denote the number of chromosomes with telomere length n at generation g, that is, K (g) n represents the concentration of telomere K (g) n . Hence we model the process (2.1) via the equation Typically, we assume that initially (g = 0), there is just one chromosome with a telomere length of Q basepairs, so that where δ n,Q is the Kronecker delta function (δ n,Q = 1 if n = Q, and δ n,Q = 0 otherwise). The early generations (g = O(1)) constitute a transient timescale during which the distribution evolves from the initial condition (2.3) to a smooth distribution. The effect of replication rule (2.1) on the distribution of telomere lengths is localised, causing the distribution to spread out (1-sided) and the mean telomere length to fall. In the next section, from (2.2) we develop a second order pde approximation, which will need two boundary conditions to be applied at small and large n. Thus it is helpful to note that the initial conditions (2.3) and evolution Eq. (2.2) discussed above together imply that K (g) n = 0 for n < 0 and K (g) n = 0 for n > Q.

Derivation of continuum model
In normal human cells, telomeres range in length from 3 to 15k basepairs and the average length of telomere lost during chromosome replication is 50-200 basepairs (Harley et al. 1990), which is much less than the initial telomere length. Thus we define Q to be the maximum telomere length, and introduce the small parameter L 1 to represent the typical fraction of telomere lost per generation (per chromosome).
Since the relative length of telomere lost per generation is O(10 −2 ), we treat telomere length (n) as a continuous variable. For the special case y 1 = 0, telomere loss y(n) = y 0 occurs at a constant rate (2.1), and we define x = n/Q, and L = y 0 /Q 1, (2.4) so that x = O(1) and the governing evolution Eq. (2.2) can be rewritten as We introduce the cumulative distribution function G (g) To account for the expected exponential growth in the population, we write where F(x, τ ) is a cumulative distribution function satisfying the initial and boundary conditions F(x, 0) = H (x − 1), F(0, τ ) = 0 and F(x, τ ) → 1 as x → +∞. The boundary conditions correspond to zero-flux conditions which ensure that telomeres of length greater than Q and less than zero cannot be formed. Using (2.7) together with (2.4), Eq. (2.6) can be rewritten as Since the term in (2.8) at τ + L has a coefficient of two, and the terms at τ have a combined coefficient of two, it is natural to expand about τ + 1 2 L. And since the term at x + L only has a weighting of one, whereas the terms at x together have coefficients summing to three, the 'centre of mass' of the template (2.8) occurs at x + 1 4 L. Thus the natural point about which to perform a Taylor series expansion of (2.8) is (x + 1 4 L, τ + 1 2 L), and performing this expansion yields the pde By expanding about (x + 1 4 L, τ + 1 2 L), the coefficients of the terms involving F τ τ and F xτ vanish and we obtain the pde. We aim to solve (2.9) subject to the boundary conditions (2.18).
An alternative derivation of (2.9) is available via asymptotic analysis: we note that the leading order approximation of (2.8) for L 1 is F τ = 1 2 F x , which has a travelling wave solution of the form F(x, τ ) = F(z, T ) where z = x + 1 2 τ and T is a new, longer, timescale. To include the next order correction terms, we introduce so that (2.8) can be rewritten as A Taylor expansion of this equation about (z, T ) gives the pde F T = 1 8 F zz , which, when converted back from (z, T ) to (x, τ ), leads to (2.9).
Another method for deriving the continuum limit of (2.8) relies on considering the behaviour of solutions which are slowly-varying in both x and τ . Since (2.8) and equations of the form F τ = cF x +φ L F x x are both linear, each possess solutions of the form F(x, τ ) = e −γ τ−βx , where the relationships between β and γ are respectively and γ = cβ − φ Lβ 2 . (2.12) Rearranging the former, and expanding for L 1, we obtain Equating this with the second expression in (2.12), we obtain c = 1 2 , φ = 1 8 and thus F τ = 1 2 F x + 1 8 L F x x as the appropriate approximation of (2.8).

The general continuum model
To describe the evolution of the distribution in the general case y 1 = 0, we introduce a new, continuous variable for telomere length, x(n), and a new timescale, τ , defined, respectively, by (2.14) The motivation for the introduction of x(n) as a measure of telomere length is that under (2.1) a telomere of length x(n) gives rise to offspring of lengths x(n) and x(n), where ( 2.15) Thus we obtain a model in which the telomere loss is no longer dependent on the telomere length, n , in place of (2.2), we have in which we have assumed that both G (g) (x) and K (g) (x) change slowly with telomere length, x. As with (2.6)-(2.7), we aim to solve (2.16) by writing The appropriate initial and boundary conditions for F are (2.18) The advantage of reformulating the problem in this way is evident in the theory following (2.8) now being applicable.

Solution of the continuum model
The solution of the pde (2.9) with boundary conditions (2.18) is a moving Gaussian of the form This cumulative distribution corresponds to the Gaussian density f = F x given by (2.20) Recalling from (2.14) and (2.7) that, to leading order, we obtain the telomere length distribution in terms of the original variables (n, g) as (2.22) In Fig. 2 we use (2.22) to plot 2 −g K (g) n against n for fixed values of g in the case y 1 = 1/60, y 0 = 50, Q = 5950. As expected, the distribution is bell-shaped and widens over time. A more detailed analysis of this figure reveals that the leftward movement slows at later times, and that the distribution is skewed to the right, that is, it is not symmetric when reflected in the mode, and the decay as n → ∞ is slower than the decay as n → 0.
Having obtained solutions to the differential Eq. (2.9) which approximate solutions of the discrete model (2.2), we use (2.21) to transform back to the physically meaningful (n, g)-coordinate systems to interpret and comment on the results. We note that K (g) n is a probability distribution function, with (y 0 + y 1 n) having a log-normal distribution, that is, log(y 0 + y 1 n) ∼ N ( μ, σ 2 ) with In terms of the telomere length, n, the mode, median, mean and standard deviation of the distribution are Under the assumption that g = O(L −1 ), we have n med = n mean = n mode to leading order. Including first order correction terms in L 1 we have (2.28) The final expression (2.27), which is also obtained using this assumption, shows that the standard deviation increases with generation number for 0 < g < 1/L and decreases for g > 1/L. We denote by K (g) tot the total number of chromosomes at generation g, by μ n (g), the mean telomere length of a chromosome at generation g, and by φ div (g), the fraction of dividing chromosomes. These quantities are defined by (2.29) Guided by (2.1), in order for cell division to occur, we require n > n c (2.1), so that all daughter chromosomes have telomeres above the threshold length. Chromosomes with telomeres of length 0 < n < n c exist but cannot replicate; we term these cells 'senescent' and note that the fraction of such senescent chromosomes is φ sen (g) = 1 − φ div (g). These cells contribute to K (g) tot and μ n (g) but not to φ div (g). To illustrate features of the distribution of telomere lengths, in Fig. 3 we plot the mean telomere length μ n (g) and the fraction of dividing chromosomes φ div (g) over many generations. From the left panel of Fig. 3 we note that, in general, the loss of telomere is in general not linear over time.
Before considering Case B for which both telomere loss and the probability of replication are telomere-length dependent, we pause to consider the case for which Left: plot of the mean telomere length over 120 generations. The thicker solid line corresponds to the case Q = 5950, y 0 = 50, y 1 = 1/60 illustrated in Fig. 2, the narrower solid line to the case y 0 = 100, y 1 = 0. In both cases, dotted lines show two standard deviations above and below the mean. Right: the proportion of dividing chromosomes φ div (solid lines) decrease over time as the telomeres shorten, and φ sen (dashed). For the case y 1 = 0, y 0 = 100, we observe the formation of senescent chromosomes around generation 100 the number of base pairs lost per replication is constant, that is, y 1 = 0 in (2.22). It is straightforward to show that in place of (2.22) we have (2.30) Thus the distribution of telomere lengths is Gaussian, with a mean μ n (t) = (Q − 1 2 y 0 t) which decreases linearly with time, and a standard deviation σ n (t) = 1 2 y 0 √ t; these results support those presented by Levy et al. (1992), who assumed that the distribution was binomial, and obtained a linear dependence of the mean with time. The proportion of senescent cells is given by (2.31) In the more general case, where y 1 > 0, the fraction of senescent cells is given by n where n c = y 0 /(1− y 1 ), by restricting replication to chromosomes with n > n c we prevent the formation of physically unrealistic daughter chromosomes with nonpositive telomere lengths.

Problem formulation
We now generalise the analysis of the previous section to situations in which the probability of replication is telomere length-dependent. In more detail, we return to the formulation of models of cell division in which generation number g is discrete, telomeres are measured in terms of the number of base pairs, n and the number of chromosomes at generation g with telomere length n is denotes by K (g) n .
We suppose that a chromosome whose telomeres have length n, divides with probability P div (n) = an + b where the constants a and b, are chosen so that With P div (n) = an + b and y(n) = y 0 + y 1 n, the reaction equation for a single chromosome can be written as ( 3.1) where 0 ≤ n ≤ Q since the parent chromosome remains in the population regardless of whether replication occurs. Taken together, these replication rules lead to the following equation for the distribution of telomere lengths on generation (g + 1) which we again solve in g > 0, for lengths 0 ≤ n ≤ Q, subject to the initial data K 0 n = δ n,Q . At the start of replication, n ≈ Q and the fraction of telomere lost per division event is y 0 /Q + y 1 . By contrast, when n is small, replication ends and the fraction of the original telomere length lost per replication is y 0 /Q. As before, we introduce a small parameter L 1 and assume y 0 /Q = O(L) and y 1 = O(L), so that telomere loss is always small. Under these assumptions, we can perform a Taylor series expansion of the governing Eq. (3.2) in terms of the small parameter, L. As in Case A (see (2.4)), we introduce the scalings where we assume that the distribution is normalised, via 1 0 f (x, τ ) dx = 1. We define the mean of the distribution f (x, τ ) by (3.6) From (3.2), the governing equation is thus (3.7) We expand the delay in the argument of f in (3.7), to obtain Performing a Taylor expansion of the last term, We now Taylor expand the small differences in the τ -argument of f (·, ·), also to O(L 2 ), to obtain the pde (3.13) Due to the x-dependence in the coefficients on the rhs of (3.7), it is not possible to simultaneously remove the terms involving both f τ τ and f xτ by expanding around (x + δL, τ + σ L) for some δ, σ . However, our choice of δ = 0 = σ means that no f xτ terms are generated.
To find an expression for ξ(τ + L)/ξ(τ ), we integrate (3.10) with respect to x, using (3.6) and f dx = 1. We assume the boundary conditions This equation can be solved for ξ(τ ), once μ(τ ) has been determined; this in turn requires some knowledge of f (x, τ ). In practice, we use (3.14) to eliminate ξ(τ ) from (3.10), which yields where D, v are given by (3.11)-(3.12) and (3.17) Although our main focus is on the range 0 < x < 1, Eq. (3.15) is well-defined on the whole of R. For simplicity we work with the boundary conditions will eventually cause f (x, τ ) to become significant in x < 0, we neglect these since they correspond to the offspring of senescent telomeres and are only formed after the population has become senescent. Hence, we also have f x , f x x → 0 as x → ±∞.
In order to specify the initial condition for (3.15), we recall that K (0) n = δ Q,n , which implies that the cumulative distribution for K Since (3.15) is second order in τ , we formally require two initial conditions. However, this equation is overdamped, and we are mainly concerned with the solutions at later times, which we obtain by solving the first-order problems obtained by truncating the asymptotic expansions in L 1. Therefore we specify only (3.19). In the next subsection, we consider the leading order equation, and derive an exact expression for the mean telomere loss over many generations. By performing a more detailed asymptotic calculation in Sect. 3.3, we derive analytic expressions for the shape of the distribution at early times.

Analysis of the leading order PDE
To O(L), Eq. (3.15) reduces to the first-order wave equation which we solve using the method of characteristics.
and μ(τ ) is determined by the characteristic that passes through x = 1 at τ = 0: which can be solved implicitly, for the general case, by Expression (3.23) is valid if y 1 = 0 and a Q y 0 = b y 1 . In the special case y 1 = 0, the solution of (3.22) is and when a Q y 0 = b y 1 , we have neither of which can be recast into explicit expressions for μ(τ ) with elementary functions. An approximation of the time to senescence is given by τ sen = τ (0) which, from (3.23), yields the corresponding expressions for the time of onset of senescence in the special cases y 1 = 0 and a Q y 0 = b y 1 are respectively. The amplitude, A(τ ), in (3.21) solves (3.28) Dividing the corresponding sides of (3.28) by (3.22) and integrating with respect to μ supplies . (3.29) Using (3.29) to substitute with A(τ ) in (3.21), we deduce where μ(τ ) is given implicitly by (3.27) if y 1 = 0 or a Q y 0 = b y 1 and Eq. (3.23) otherwise. Returning to (3.14), now with expressions determining the mean μ(τ ), we evaluate the population size, ξ(τ ), by summing the leading order terms over 0 In the special case y 1 = 0, we have where Li 2 (·) is a dilogarithm ( In Fig. 4, we show how, for Case B, the mean μ(τ ) and the total number of chromosomes ξ(τ ), as defined by Eqs. (3.23) and (3.32), change over time when L 1. We note that when b = 1, a = 0 we have P div = 1 and so recover Case A; when the telomere loss is constant, with y 0 = 100 bp per generation (narrow dashed line) we observe exponential growth of the population, and an approximately linear reduction in average telomere length, confirming the earlier analysis. The thick dashed line, where telomere loss is constant, and P div = 0.2 + 0.8(n/Q) also gives approximately linear loss in average telomere length, but now with a rate of population increase which slows at later generations. The two cases for which both telomere loss and replication probability are telomere-length dependent yield population growth curves which slow at later times. Additionally, the rate at which the average telomere length shortens decreases with generation number, and also show average telomere lengths decreasing with time in a nonlinear fashion.
From the first-order pde (3.20), we have found approximations to the average telomere length, μ(τ ) and the number of chromosomes, ξ(τ ), for the general case where both the probability of replication and amount of telomere lost depend on telomere length (with L 1). In order more accurately to describe the distribution in the next subsection we retain second-order derivatives in our pde approximation to the underlying discrete model.

Early time asymptotic analysis of second-order pde
Our analysis of the (leading order) first-order pde (3.20) supplies expressions for the number of chromosomes, ξ(τ ), and the mean telomere length, μ(τ ) given by (3.22)-(3.23). Since the distribution is given as a Dirac δ-function, these expressions do not provide information about the variance of the distribution. In order to determine how the distribution evolves over time, we introduce new variables (z, T ) which are chosen to retain the second-order spatial derivative from (3.15). Thus we define In terms of the new variables, the pde (3.15)-(3.17) becomes where primes denote derivatives with respect to the argument, and the leading order (in L 1) expressions for the parameters are (3.37) To satisfy Eq. (3.36) at leading order, we note that v + θ μ = O(L 1/3 ), from (3.22) and the definitions of θ, v above. However, there will be some correction term, which, for the simplicity of later calculations, we write as v + θ μ = L 1/3 λ (T ). The function λ(T ) will be determined via the normalisation condition f dz = 1. Under these assumptions, Eq. (3.36) reduces to where D = D − θμ 2 . Since T = O(1) corresponds to an initial short timescale, in which τ = O(L 1/3 ), the mean telomere length is given by μ(τ ) = 1. In this case, the coefficients in (3.38) take the constant values θ = 1+b+a Q, D = (b+a Q)( y 0 + y 1 ) 2 and D = (b + a Q)( y 0 + y 1 ) 2 /(1 + b + a Q).
Taking the Fourier transform of (3.38) with f (k, T ) = e ikz f (z, T ) dz, we obtain where the initial data and global constraint (3.35) imply f (k, 0) = 1, and f (0, T ) = 1. The general solution of (3.39) is where C(·) is an arbitrary function and, without loss of generality, we assume λ(0) = 0 = λ (0). The initial condition and global constraint (3.35) imply C(q) ≡ 1 and λ(T ) = Da QT 3 /6θ so that the solution of the Fourier Transform Eq. This distribution is a Gaussian pulse which converges to a Dirac δ-function as T → 0 + , and which spreads with the standard z ∼ √ T scaling as T → ∞. Recalling (3.37), we obtain . Over the longer timescale (τ = O(1)), the distribution becomes skewed as shown in the numerical solution presented in Fig. 5, which is due to the effect of the x-dependence of the coefficients (3.11)-(3.17), in (3.15).

Summary
In Sects. 2 and 3 we derived two models for the evolution of telomere lengths over time in a population of replicating chromosomes. We used asymptotic methods to determine the shape of the distribution in various cases. For Case A, all chromosomes replicate, with probability one on each generation, provided their telomeres are sufficiently long. and the distribution has a Gaussian or a log-normal shape. For Case B, the probability of chromosome replication depends on telomere length. We obtain an expression for the evolution of the mean for all times, and show that at early times, the distribution is Gaussian; however, the scaling for the dependent variables (3.34) is nonstandard. The analyses for Cases A and B describe a population of chromosomes with each replicating independently. In practice, this description is overly simplistic: cells containing multiple (N = 46) chromosomes replicate and the checkpoints on replication occur at the cell level. Hence in the next section we incorporate this level of complexity and consider a population of cells, each of which contains N chromosomes.
contains N -chromosomes. Our model describes the distribution of telomere lengths of the chromosomes within each cell. In this section we consider the fully deterministic Case A, and in Sect. 5 we consider the more general probabilistic model, Case B, from Sect. 3.
We will assume that a cell can only replicate if all of its N = 46 chromosomes have telomeres which exceed the threshold length; if the length of any telomere falls below the threshold, the cell becomes senescent and it will not replicate further, although it will remain a viable member of the cell population. During cell division, each chromosome replicates, producing two daughter chromosomes, which are allocated to the two daughter cells randomly and independently of the other replicating chromosomes. In contrast to the models presented in Sects. 2 and 3, it is now highly likely that the total telomere length of each daughter cell will be less than that of their parent. We consider the same rules for chromosome replication as in (1.1); we use m, to denote the total telomere length in the cell and Y (m), to denote the total amount of telomere lost during each replication event.

Development of discrete cell-scale model
As in (2.2), when a cell divides, each chromosome produces one daughter chromosome with the same telomere length as the parent and a second with shorter telomere length. When specifying replication rules, we number chromosomes in a cell using r , with 1 ≤ r ≤ N . Following (1.1), we write where n r denotes the telomere length of the r th chromosome (1 ≤ r ≤ N ) in a given cell, g represents the generation number, and y, the amount of telomere lost from one daughter chromosome upon replication. By analogy with (1.2) we define y(n r ) = y 0 + N y 1 n r , where y 0 , y 1 are constants; this form is chosen to simplify later calculations. Combining N chromosomes into one cell, and describing the state of each cell solely by its total telomere length, m, and generation number, g, we define with the aim of summing (4.1) over N chromosomes, to form a replication rule for cell C (g) m . We now simplify the telomere loss rule. A model that includes information on the telomere length of each chromosome in a cell would not be analytically tractable. Since we propose to only retain information on the total telomere length of the cell, we assume telomere loss Y (m) depends only on the total telomere length in the cell, m. Thus, for each chromosome (r ), we write which replaces the term y(n r ) in (4.1). To confirm this expression is consistent with the individual loss term y(n r ) = y 0 + N y 1 n r , we calculate the maximum total loss of telomere, which is given by N Y (m (4.4) If we assume that shortening of the r th chromosome can be inherited by either daughter cell, then there are 2 N different ways of allocating the longer and shorter daughter chromosomes to the two daughter cells. If we assume that each arrangement is equally likely, then we have (4.5) by combining (4.1) with (4.2). The number of shortened chromosomes a daughter cell inherits determines the amount by which its telomere length is reduced, this can vary between zero and N Y (m) bps. Cells with telomere length m can be the offspring of cells with telomere length (m + j y 0 )/(1 − j y 1 ) for any 0 ≤ j ≤ N . Summing (4.5) over all possibilities, weighted according to the corresponding probabilities, gives a discrete evolution equation for C (g) m , the number of cells at generation g with total telomere length m, The simplification is due to the summand being invariant under the transformation j → N − j. Typically, we solve this model subject to the initial conditions C (0) m = δ m,N Q , which corresponds to a single cell with total telomere length m = N Q at generation zero.

Deterministic continuum cell-level model
We start by reformulating (4.6) in terms of the cumulative distribution function G (4.7) If we permit telomeres of negative lengths to form (i.e. m < 0), then the number of telomeres doubles every generation and the boundary conditions for the cumulative distribution function are (4.8) since, at the start of the simulations we have one cell with total telomere given by m = m 0 = N Q. To determine the shape of the distribution, we rescale to remove the exponential growth by writing G The dimensionless small parameter, L, is defined by the small initial relative rate of loss of telomere, that is We introduce rescaled parameters Following (2.14), we replace the discrete variable m (with 0 ≤ m ≤ N Q) by the continuous variable x, where It has already been noted that cells with telomere length m = (m + j y 0 )/(1 − j y 1 ) give rise to daughters of length m; correspondingly, daughter cells with parameter x arise from parents with parameter x given by where we have used log(1 − h) ∼ −h − 1 2 h 2 for h 1. Combining the asymptotic scalings (4.10) with the change of variable (4.11), the expansion (4.12), τ = Lg and G (g) The discrete differences in this model can, with good accuracy, be replaced by derivatives; using Taylor's series, based on L 1 and j y 1 L/N 1. Here we have introduced the shifts x → x − 1 4 y 1 L (δ = O(1)) and τ → τ − 1 2 L so that when we perform the Taylor series expansion in L 1, the coefficients of F τ τ or F xτ vanish. Thus we obtain We solve (4.14) subject to the initial conditions F(x, 0) = H (x −x max ), where H (·) is the Heaviside function, and the boundary conditions F(x, τ ) → 0 as x → −∞ and F(x, τ ) → 1 as x → +∞ which together imply (4.15) Equation (4.15) corresponds to the cumulative distribution function for a Gaussian distribution which has density (4.16) Rewriting (4.16) in terms of telomere length, m, and generation number, g, we deduce that the distribution of telomere lengths is given by a log-normal distribution. In more detail, we have C (4.17) In the special case y 1 = 0, we define x = m/N Q so that 0 ≤ x ≤ 1, and in place of (4.12), we have x = x + j y 0 /N Q = x + j y 0 L/N . When y 1 = 0, Eqs. (4.14) and (4.17) are replaced by (4.18) and which gives a standard Gaussian distribution, with a mean that reduces linearly over time, and a standard deviation which increases with √ g.
The solutions (4.19) and (4.17) are similar to the chromosome-level results (2.22), which are plotted in Fig. 3. Before illustrating the behaviour of the cell-level model, we explore the distribution of telomere lengths within each cell so as to determine when senescence occurs.

Determination of time to senescence
A cell will not divide if the telomere length of any of its chromosomes falls below the critical length, here taken to be zero. To determine when a cell reaches senescence, we must determine the length of its shortest telomere. We now develop a method for approximating the length of the shortest telomere in a cell so that our model can predict the time of the onset of senescence.
In each cell there is a subpopulation of N telomeres, whose lengths at generation g we define by {S (g) j } N j=1 . Since daughter telomeres are randomly allocated to daughter cells, the length of each telomere S (g) j is effectively given by a random number R/N where R is drawn from the distribution C If we randomly pick a sample of N such variables, label them {r 1 , . . . , r N } and then order them so that x i = r j for some permutation such that x 1 < x 2 < x 3 . . . < x N , then x 1 = min j {r j }. The minimum of the sample x 1 is distributed according to the first order statistic.
Given the probability density function f (r ) and the cumulative distribution function (cdf) F(r ) for the variables r j , the probability density function f 1 and cumulative distribution function F 1 of the first-order statistic x 1 = min j {r j } (Hogg and Craig 1970) are To illustrate, we consider a Gaussian distribution where the pdf and cdf of a N (μ, σ ) are given by In the left panel of Fig. 6, the narrow line represents the pdf f (r ) for the case of zero mean and unit standard deviation, N (0, 1). The thicker line represents the distribution . Right: similar plots for the log-normal distribution where log X ∼ N (μ, σ ) with μ = 0 ( f = exp(−(log x − μ) 2 /2σ 2 )/xσ √ 2π , which is only defined for x > 0). The pdf (thin line) for the case σ = 1 has a maximum near x = 0.4 and the pdf of the first order statistic (4.20) has a maximum near x = 0.1; the case σ = 0.25 has a pdf with maximum just below x = 1 and the first order statistic has a maximum around x = 0.5 of the first order statistic, f 1 (x 1 ) given by (4.21), that is, x 1 = min j {r j } and {r j } N j=1 are N = 46 sample random variables sampled from the distribution f (r ) (4.20). The distribution of the first order statistic f 1 (x 1 ) is noticeably shifted to the left when compared to f (r ). It also skewed and has a smaller variance.
The mean of the first order statistic can be calculated as E[ f (1) (x 1 )] = −2.216. Whilst the standard deviation of the normalised Gaussian is unity, that of the first order statistic with N = 46 is 0.469-significantly smaller than that for the population as a whole. The mode of the pdf of the first order statistic, f 1 (x 1 ), is given by the maximum of f 1 (x 1 ), which occurs when f 1 (x 1 ) = 0, that is Solving Eq. (4.22) numerically in the case N = 46 yields x 1 = −2.084. The median of the distribution f (1) (x 1 ) is given by F (1) (x 1 ) = 1 2 , which implies erfc(x 1 / √ 2) = 2 1−1/N ; solving this equation numerically, gives the value x 1 = −2.171 when N = 46. Thus, the median and mode are similar, but not identical, to the mean value. To summarise, for a Gaussian distribution of telomere lengths, the length of the shortest of N = 46 telomeres is just over two standard deviations below the mean.
We now return to the full solution for telomere length distributions (4.17) and (4.19). Figures 7 and 8 show that when y 1 = 0, the telomere length decreases over time at a uniform rate; when y 1 > 0 the shortening rate reduces in the later generations. We note also that when y 1 = 0, the distribution spreads considerably, whilst when y 1 > 0, the increase in variance almost ceases at later times. This is particularly noticeable in the behaviour of the first order statistic shown by the narrower lines in Fig. 7. The fraction of dividing and senescent cells over time is shown in Fig. 8, where the curves which decrease from unity as g increases correspond to φ div (g), namely the dividing fraction and those that increase from zero show the senescent proportion, φ sen (g). The chromosome-level definition of senescence is given by Eq. (2.31) and the text immediately following it; for the cell-level description, we define the pdf of total telomere length by f (m, g) = C with the cumulative distribution function being given by F(m, g) = m j=0 f ( j, g). We next compute the pdf of the first order statistic, as f 1 (m, g) = N f (m, g)(1 − F(m, g)) N −1 using (4.20), then the fraction of senescent cells is given by φ sen (g) = m c m=0 f 1 (m, g) where m c = N y 0 /(1 − y 1 ). We denote φ div (g) = 1 − φ sen (g). From Fig. 8 we note that the cell-level definition of senescence gives a much more abrupt transition from dividing population to senescence, and that this occurs at a slightly earlier time than predicted by the chromosome-level models.

Development of continuum cell-level model
In this section we upscale the probabilistic model from Sect. 3 to a cell-level description. We suppose that the amount of telomere lost per replication and the probability of replication depend on telomere length. The probability replication is P div = am + b where a, b are constants chosen to ensure that 0 ≤ P div ≤ 1. The amount of telomere lost per chromosome is y 0 + my 1 , so the discrete cell replication rule can be written as Since, after summing over j, the two terms in square brackets generate identical contributions, for simplicity we take double the first term. In order to determine which parental cells in generation g with total telomere lengths m yield a daughter cell in generation (g +1)with total telomere length m, we replacing m in (5.1) by m + Y j ; Eq. (5.1) can then be simplified to Here, following (3.4) and (4.10), we make use of the scalings m+L Q j( y 0 + y 1 m/N Q)(1+L j y 1 /N ) .
We now introduce new variables, x and τ defined by so that 0 ≤ x ≤ 1 and τ = O(1) is the timescale over which cells transit from their initial state to senescence. The assumed form for the state variable C (g) m separates the behaviour into two components: ξ(τ ) which accounts for the total size of the population and is expected to exhibit rapid growth; and f (x, τ ), which describes the shape of the telomere length distribution. By inserting the ansatz (5.6) into (5.5), we obtain the discrete difference equation To separate the determination of ξ(τ ) and f (x, τ ), we impose a normalisation condition on f and define the mean of the distribution by Then we have that Integrating (5.7) with respect to x and using (5.8)-(5.9), we obtain (5.10) We now expand (5.7) using Taylor series with L 1, and substituting in (5.10), we obtain As with the model considered in Sect. 3, this pde cannot be solved explicitly, so we first consider the first order approximating pde and then the early time asymptotics.

Analysis of probabilistic cell model
Ignoring the O(L) terms in (5.11), we obtain a first order pde (5.15) which can be solved by the method of characteristics to obtain the following approximate equation for the mean Imposing μ(0) = 1, this ode has the implicit solution (5.18) and if a N Q y 0 = b y 1 , then Following Sect. 3.3, we obtain an approximate expression for the shape of the distribution at early times by introducing the rescalings with μ = 1. Then (5.11) supplies (5.21) where the effective diffusivity D is given by We aim to solve (5.21) subject to the conditions As with Eqs. (3.36), (3.38), we note that the leading order terms in (5.21), namely those in the square brackets cancel, leaving (5.24) where λ (T ) represents the (as yet unknown) first correction term from the term in square brackets. If one includes a term of the form λ 2 (T )z f z , then the constraints (5.23) imply that λ 2 = 0, so we omit such a term from the following analysis. Taking Fourier transforms ( f (k, T ) = e ikx f (z, T ) dz), we obtain (5.25) which has the solution where, without loss of generality, we impose λ(0 and f (x, τ ) = π(1 + b + a N Q) (5.28) We note how the generalisation to multiple chromosomes per cell changes the rates at which the distribution of telomere lengths shortens and spreads out. To compare these effects, we introduce two quantities: a loss rate, which gives the rate at which the mean of the telomere length distribution reduces over time, and a 'diffusivity', corresponding to the rate at which the distribution widens in terms of the telomere length variable; note that this is not a literal spreading out in physical space. Firstly, we note that the advection terms are very similar for multiple and single chromosomes (N = 1): defining the loss rate in (3.15)-(3.12) by v/θ and similarly in (5.11) and (5.13)-(5.14), we have which can be made identical under the mapping a → a N . Next, we compare the rate at which the distribution spreads in the cell model with that in the chromosome model. We find this by applying the transformation (5.20) and then examine the coefficient of f zz in he In the cell model, the diffusivity is given by with D, μ given by (5.14), (5.16) respectively, differs from the corresponding quantity in the chromosome model, where D chromo = (D − θμ 2 )/2θ with D, θ, μ given by (3.16), (3.11), (3.22). Evaluating these expressions at x = μ, we have In the case N = 1, the expression for D cell reduces to D chromo . Even after applying the mapping a → a N , we see that for N > 1, D cell < D chromo , due to the mixing of longer and shorter offspring chromosomes into the daughter cell caused by (5.1), in which the dominant terms occur for j ≈ N /2. The effect of having multiple chromosomes in each cell, and randomly allocating daughter chromosomes to each daughter cell is to reduce the relative spread of telomere lengths, which leads to a more abrupt transition from replication to senescence.
In Fig. 9, we use (5.31) to plot both D chromo and D cell against generation number g. The chromosomal case shows little dependence on parameters, and reduces from about 0.8 to approximately 0.2. The diffusivity in the cellular case is much smaller (D cell < 0.1) and, for most of the parameter values illustrated in Fig. 9, D cell increases over the early generations and decreases at later times. The initial increase is due to the term in square brackets in (5.31), which is approximately 1 2 (1 − P div ). Given gen, g D chromo and D cell , we can create a proxy for the standard deviation of the distribution, σ , by noting that (d/dτ )(σ 2 ) = D, hence σ = √ τ 0 D dτ . This expression yields the curves plotted in Fig. 10. On the left, against generation number we plot the average telomere length of a chromosome in base pairs, Qμ. On the right, we plot the approximate standard deviation √ D dτ , to illustrate the significant difference between the cellular and chromosomal models, and smaller differences caused by changing the parameters a, b.

Discussion
In this paper, we have considered the dynamics of telomere loss in a population of independent chromosomes and a population of cells, each containing a population of chromosomes. Our aim has been to understand the kinetics of telomere loss and onset of senescence, which are hallmarks of aging. During each replication event some telomere is lost from one daughter chromosome, and senescence occurs when the telomere length becomes too short for replication to occur. Since telomere length plays an important role in cell division, we have considered a variety of models, which focus on both the amount of telomere lost, and the rate (probability) of cell division depending on telomere length.
We have developed fully discrete models, and approximated them by deriving continuum level pde descriptions. This procedure relies on only small fraction of total telomere length being lost in each generation, so that an asymptotic reduction yields governing equations both at the level of individual chromosomes, and at the cell-level, where each cell contains a subpopulation of N chromosomes. As an example, we consider human cells, where N = 46. We have analysed these models to determine the evolution of the telomere length distribution over many replication cycles (generations). In all cases we observe a shortening of telomeres, until senescence occurs. In cell models senescence occurs when the shortest telomere can no longer undergo replication.
In summary, we have outlined simpler models, which describe a population of individual chromosomes, and more complicated, cell-level, models in which each cell is assumed to contain a subpopulation of N = 46 chromosomes, and upon replication each chromosome divides, giving one offspring to each daughter cell. In both cases, we have considered deterministic models as Case A, where replication occurs each generation, provided senescence has not yet been triggered; and, Case B, a stochastic model in which replication occurs with a probability that may depend on telomere length. In each of these four subcases, there is then the possibility that telomere loss is fixed, or dependent on telomere length. We summarise our results in increasing complexity: -In Case A, when telomere loss is constant, the mean telomere length decreases linearly with generation number, until senescence is reached. This occurs in both the chromosome-level and the cell-level models. In both cases, the distribution of telomere lengths is Gaussian (normal) and with a standard deviation which increases with the square root of the generation number. -In Case A and where the telomere loss is dependent on telomere length, the mean telomere length still reduces, albeit in an exponential rather than a linear fashion being given by (2.24)-(2.26). The distribution has the form of a shifted log-normal (2.22). The two forms of Case A are compared in Fig. 3 for the chromosome level model. -Various parameter values of the stochastic replication model, Case B, are illustrated in Fig. 4. In all cases the reduction in telomere length is not constant, and not given by a simple exponential, rather by the more complicated implicit formula (3.23). In this case considering the simpler case of a size-independent telomere loss rate (y 1 = 0) does not simplify the governing Eq. (3.20). The population size grows subexponentially in this case. The distribution for Case B appears almost Gaussian, however, Fig. 5 shows that it is slightly skewed. In this case the pde describing the shape of the length distribution (3.20) of telomeres does not admit explicit analytical solution. However, an approximation enables us to determine the mean telomere length. Through further use of asymptotic techniques, we show that at very early times, the distribution is Gaussian, however, at later times, it becomes skewed. -The main distinction between the results of the chromosome-level model and the cell-level model is that the latter gives a sharper transition between the whole population being in a replicative state and senescence, as illustrated in Fig. 8.
In all cases, we observe that the distribution of telomere lengths evolves in a manner consistent with corresponding Monte Carlo simulations presented previously (Qi et al. 2014), where fits to the experimental data of Zhang et al. (2000) are illustrated. By varying the telomere loss or replication probability parameters, mean telomere length reduction can be tuned to fit a variety of experimental data. Chromosome-level models similar to our Cases A and B have been considered by other authors. When Levy et al. (1992) studied a special case of Case A, for which telomere loss is constant (independent of telomere length). By assuming that the number of deletions occurring increased with the number of generations following a binomial distribution, they found the mean number of deletions to increase linearly over time, and obtained a sharp transition from replicative to the senescent state. The model of Buijs et al. (2004) is broadly similar to the more general form of Case A that we use. They fitted the model to an experimental distribution of telomere lengths, verifying the length-dependent loss. Our analysis extends this work by showing that the shape of the distribution is log-normal. Portugal et al. (2008) considered the probability of a cell's replication being linearly dependent on telomere length-as in our Case Bbut with constant amount of telomere lost per replication. Our work extends this to the case of length-dependent loss, and we have also considered the fraction of senescent chromosomes.
As far as we are aware, no existing models consider the distribution of telomere lengths within individual cells. We have extended all these models to consider a population of cells, each of which contains a subpopulation of N = 46 chromosomes. In these cell-level models, the longer and shorter daughter chromosomes from replication are randomly allocated to the daughter cells. Whilst some cells inherit predominantly longer chromosomes, and a few shorter chromosomes, the majority of cells inherit approximately equal numbers of longer and shorter daughter chromosomes. The overall effect of this is to reduce the variance of total telomere length in the cell-level model (compared to the individual chromosome description, which can be thought of as a cell-level model with N = 1). Analysis of cell-level continuum models gives some results similar to the individual chromosome model. The mean of the distribution is in good agreement with the Monte Carlo simulations of cell-level models presented previously (Qi et al. 2014). However, the transition between replicative and senescent states is sharper in the cell-level model than in the chromosome model, as noted by Qi et al. (2014). This feature is caused by the length of the shortest telomere in a cell having a smaller variance than the telomere length of overall population, as shown in Fig. 6. Since a cell becomes senescent when its shortest telomere cannot replicate, the total telomere length of senescent cells will increase with the number of chromosomes per cell. In the analysis of the cell-level model we have used the theory of order statistics to determine the expected time at which the shortest telomere reaches the threshold of zero length and hence determined the time at which the onset of senescence occurs. We expect the cell-level results presented here to be the same as if we simulated a population of individual chromosomes and then randomly allocated the telomeres into groups of N = 46. This is because we have assumed a random allocation of long and short daughter telomeres into each cell, as given by (4.7). However, in the case of stem cells, for example, the allocation of long and short telomeres is not random, and so (4.7) would not be valid, and should be replaced with an alternative, which would result in different outcomes.
In future work, we plan to generalise these results to include the effect of telomerase -an enzyme which lengthens telomeres and so allows cells to continue replicating (Qi 2011;Qi et al. 2019); the possibility of varying activity levels of telomerase was discussed by Epel et al. (2004). It may also be possible to generalise the models proposed here to include other mechanisms which influence telomere length, such as the Alternative Lengthening of Telomeres (ALT), telomerase activity, telomere recombination, and Werner's syndrome-a disease in which accelerated aging is caused by the corrupted replication of chromosomes. This is caused by the failure to resolve stalled forks in the DNA replication process. Muraki et al. (2012) discuss the repair mechanism for double-stranded breaks, and other problems that can occur at the microscopic level, for example, the role of deficiencies of the Shelterin proteins, and the interplay between stochastic telomere loss mechanisms and chromosome instability. Oxidative stress is known to accelerate aging at the cellular level, and many other environmental factors have been shown to correlate with telomere length. Starkweather et al. (2013) discuss the effects of various psychosocial factors on telomere length, including chronic psychological stress, sleep quality, socioeconomic status, educational attainment and genetic components. The mathematical modelling of how these macroscopic factors influence telomere loss at the microscale remains an open problem.