# Stochastic branching models for the telomeres dynamics in a model including telomerase activity

Athanase Benetos, Coralie Fritsch, Emma Horton, Lionel Lenotre, Simon Toupance, Denis Villemonais

► To cite this version:

Athanase Benetos, Coralie Fritsch, Emma Horton, Lionel Lenotre, Simon Toupance, et al.. Stochastic branching models for the telomeres dynamics in a model including telomerase activity. Journal of Mathematical Biology, 2025, 91 (8), 10.1007/s00285-025-02230-z. hal-04648211v2

HAL Id: hal-04648211

https://hal.science/hal-04648211v2

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Stochastic branching models for the telomeres dynamics
in a model including telomerase activity

A. Benetos<fnref n="1" />,<fnref n="2" />, C. Fritsch<fnref n="3" />, E. Horton<fnref n="4" />, L. Lenotre<fnref n="5" />,
S. Toupance<fnref n="1" /> and D. Villemonais<fnref n="6" />,<fnref n="7" />

27th August 2025

Abstract

Telomeres are repetitive sequences of nucleotides at the end of chromosomes, whose evolution over time is intrinsically related to biological ageing. In most cells, with each cell division, telomeres shorten due to the so-called end replication problem, which can lead to replicative senescence and a variety of age-related diseases. On the other hand, in certain cells, the presence of the enzyme telomerase can lead to the lengthening of telomeres, which may delay or prevent the onset of such diseases but can also increase the risk of cancer.

In this article, we propose a stochastic representation of this biological model, which takes into account multiple chromosomes per cell, the effect of telomerase, different cell types and the dependence of the distribution of telomere length on the dynamics of the process. We study theoretical properties of this model, including its long-term behaviour. In addition, we investigate numerically the impact of the model parameters on biologically relevant quantities, such as the Hayflick limit and the Malthusian parameter of the population of cells.

**Keywords:** branching processes, non-conservative semi-groups, many-to-one formula, quasi-stationary distribution, Hayflick limit, telomerase.

MSC: 60J80, 60J85, 60K40.

# 1 Introduction

Telomeres are repetitive sequences of nucleotides located at the ends of linear chromosomes that act as protective caps to ensure the genomic stability. As cells divide, telomeres gradu-

<footnote n="1">¹Université de Lorraine, Inserm, DCAC, F-54000 Nancy, France</footnote>

<footnote n="2">²Université de Lorraine, CHRU-Nancy, Service de Gériatrie, F-54000 Nancy, France</footnote>

<footnote n="3">³Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France</footnote>

<footnote n="4">⁴Department of Statistics, University of Warwick, Coventry, CV4 7AL, UK</footnote>

<footnote n="5">⁵Université de Haute-Alsace, Faculté des Sciences et Techniques, F-68200 Mulhouse, France</footnote>

<footnote n="6">⁶Université de Strasbourg, IRMA, F-67000 Strasbourg, France</footnote>

<footnote n="7">⁷Institut universitaire de France (IUF)</footnote>
ally shorten, ultimately reaching a critical length triggering cellular senescence or apop-
tosis (Olovnikov, 1996; Entringer et al., 2018; Xu and Teixeira, 2019), processes implicated
in a number of age-related diseases such as cardiovascular, metabolic and neurodegenerative
diseases (D'Mello et al., 2015; Haycock et al., 2017; Toupane and Benetos, 2019). In addi-
tion, research has shown that factors such as lifestyle, genetics and oxidative stress can also
impact telomere length dynamics (Starkweather et al., 2014; Shammas, 2011; Louzon et al.,
2019).

On the other hand, some cell types can express an enzyme called telomerase. This
DNA polymerase has the ability to add telomeric repeats to the end of telomeres during
DNA replication, which can compensate the effects of telomere shortening but also plays an
important role in the development and evolution of cancer (Jafri et al., 2016; Blasco, 2005).

In the past 30 years, mathematical modelling has played an important role in under-
standing the long-term behaviour of the distribution of telomere length in cells and the
associated health implications. One of the earliest (deterministic) models for the evolution
of telomere lengths was given by Levy et al. (1992), which mirrored the results seen in *in vitro*
experiments. Shortly after, Arino et al. (1995) and Olofsson and Kimmel (1995) reframed
this model in terms of a branching process to obtain both exact and asymptotic results for
the behaviour of telomeres. Since then, there have been a wide range of stochastic models
for telomere lengths, both with and without the presence of telomerase, for the purpose
of understanding the role of the shortest telomere in senescence (Bourgeron et al., 2015;
Dao Duc and Holcman, 2013), calculating the time until senescence (Eugène et al., 2017),
parameter inference (Lee and Kimmel, 2020) and calculating the stationary distribution (Lee
and Kimmel, 2020). We also refer the reader to the articles of Portugal et al. (2008); Wattis
et al. (2020); Qi et al. (2014) for computational analysis of both deterministic and stochastic
models.

In this article, our first main contribution is to propose a mathematical model that de-
scribes the evolution in continuous time of a population of cells, whose dynamics depend
explicitly on the lengths of their telomeres. This is in accordance with recent empirical evid-
ence (Bourgeron et al., 2015; Dao Duc and Holcman, 2013; Xu et al., 2013), which suggests
a strong link between the length of the shortest telomere and the cell's behaviour. The main
novelty in our approach is that we allow for different cell types, multiple chromosomes per
cell *and* we also take into account the effect of telomerase. Moreover, we model the dynamics
of the *whole population* of cells, which requires one to keep track of the state of all telomeres
in each individual cell within a growing cell population. To the best of our knowledge, this is
the first mathematical model that takes all of these components into account simultaneously.

Among the models in the aforementioned literature, some of them describe the whole cell
population, while others describe a single cell lineage. While both approaches are relevant,
we emphasize that the latter approach does not appropriately describe the distribution of
the whole population when reproduction and/or senescence rates of a cell depend on its
telomere lengths. This is due to an imbalance of mass compared to the whole population of
cells. To counteract this issue, one needs to appropriately weight the single linear process, as
in the *many-to-one formula*. In the setting of telomere dynamics, this is given in Lemma 3.3,
which allows one to represent the first moment of the entire cell population by the average
behaviour of a typical trajectory in the population, appropriately weighted. This result
allows one to more easily analyse the behaviour of the branching process, and provides a
useful tool for numerical simulations.

One of the main mathematical contributions of this paper, Theorem 4.5, is related to the asymptotic behaviour of the first moment of the particle system. In the case that telomerase is present and compensates the attrition of telomeres, so that the population grows indefinitely, we show two things. Firstly, the distribution of an average cell and its telomere lengths converges with time, and secondly, the first moment of the population size grows asymptotically exponentially fast. This analysis is related to the theory of quasi-stationary distributions (see e.g. Collet et al. (2013); van Doorn and Pollett (2013); Méléard and Villemonais (2012)). While the analysis presented in this article is fairly standard in the literature on quasi-stationary distributions, to the best of our knowledge, this is the first time these techniques have been applied in the context of telomere dynamics.

We also present numerical simulations of the model, both with and without telomerase.
In the case where it has very little impact or is not present at all, the population eventually
stabilizes due to the fact that cells can no longer divide when their telomeres are too short.
In this case, the number of times the population doubles is finite; this effect is also called the
*Hayflick limit* (Hayflick and Moorhead, 1961; Shay and Wright, 2000; Trentesaux and Riou,
2010). Thus, we investigate the influence of the model parameters on the Hayflick limit.
On the other hand, when telomerase is present, we study its influence on the asymptotic
growth rate of the population and the stabilization of the distribution of an average cell and
its telomere lengths, illustrating Theorem 4.5 in practice.

The rest of the paper is set out as follows. In Section 2, we present the specific model that we will work with throughout the paper, presenting first the biological mechanism for DNA replication at the level of chromosomes, followed by the stochastic process describing the population of cells. This model takes into account the telomere length distribution of chromosomes, as well as the cell's ability to divide or not. In addition, we allow for the possibility that the dynamics depend on the absence/presence of telomerase. Section 3 is devoted to the mathematical analysis of the first moment of the mathematical model, leading to the *many-to-one formula* (see Lemma 3.3). In Section 4, we establish and prove our main mathematical result, which characterises the long-term behaviour of our process (see Theorem 4.5). In particular, we describe the growth rate of the average number of particles in the system and the corresponding stationary distribution. Section 5 contains the numerical simulations, which illustrate the impact of the model parameters on 1) the Hayflick limit (the number of times the population doubles), 2) the asymptotic growth rate of the population, and 3) the asymptotic distribution of telomere lengths in the cell population. Finally, in Section 6, we discuss our model, results and possible future work and extensions. In addition, we alert the reader to the table of notation in Appendix C, which collects the main parameters used for the simulations in Section 5.
# 2 The model

## 2.1 The biological model

In this section, we describe the biological DNA replication mechanism we adhere to throughout the article. For concreteness, we will describe the replication mechanism that occurs in human cells, however, the mathematical model we consider in the next section is robust enough to apply to other organisms. We will also only describe the relevant aspects of the DNA replication process, and we refer the reader to the articles of Grach (2013); Aubert (2014) and references therein for further details.

Human cells contain 23 pairs of chromosomes, each composed of a DNA strand whose extremities are called telomeres. Each strand of the DNA double helix has a 3' and a 5' end. Due to the antiparallel structure of the double helix, the 3' end of one strand opposes the 5' end of the other strand. Moreover, the 3' end overhangs the opposing 5' end (see Figure 1a).

Starting from the centre of the chromosome (referred to as the origin of replication), the DNA helix is 'unzipped' in the direction of the telomeres, creating two 'Y' shapes, each called a replication fork (see Figure 1b). The two separated strands will then act as templates for an enzyme (DNA polymerase) to make new complementary DNA strands, resulting in two copies of the original chromosome.

The orientation of the DNA helix means that when it is unzipped, in one direction it unzips from the 3' to the 5' end, which will form the *leading* strand, while in the other direction the helix is unzipped in the 5' to 3' direction, and will form the *lagging* strand. Due to the way DNA polymerase replicates DNA and the antiparallel nature of a DNA helix, the mechanism is different for each of the two strands. Indeed, as the DNA polymerase synthesises the leading strand, DNA fragments are added continuously to form the new complementary (antiparallel) strand from the origin to the end of the 5' end of the template. Note that the new complementary strand will be shorter than the original one since DNA polymerase can only use the 3' to 5' strand as a template, resulting in the loss of the original overhang. In order to re-establish the overhang, an enzyme called nuclease removes the end of the template 5' end so that the new 3' is longer. On the other hand, when DNA polymerase synthesises the lagging strand, it does so in separated segments, called Okazaki fragments, using a 'back-stitch' type method. Due to the discontinuous nature of the replication along the lagging strand, DNA cannot be replicated all the way to the end of the template strand, hence also resulting in an overhang of the 3' end. We refer the reader to Figures 1a, 1b and 1c for a diagram of these steps, as well as the articles of Grach (2013) and Muraki et al. (2012) for further details and references.

Thus, as a result of the mechanism described above, telomeres become shorter with each cell division. This is known as the *end replication problem*. When the shortest telomere in a cell becomes too short, the cell is unable to divide any further (see Bourgeron et al. (2015); Hemann et al. (2001)), since another division would risk damaging the DNA. When this occurs, the cell becomes senescent, which we refer to as deactivation in the rest of the paper. Thus, if the shortest telomere in each cell falls below a certain threshold, the number of cells in the population ceases to increase. This leads to a concept called the Hayflick limit, which is defined as the number of times a population of cells can double (before cell division is no
longer possible).

In certain cells, such as stem cells, the majority of cancer cells and some somatic cells, an enzyme called telomerase is present, which provides a mechanism by which telomeres can lengthen. In these cells, after DNA replication has occurred, telomerase has the ability to add new nucleotides to the 3' end. After a certain number of nucleotides have been added, the complementary sequence of nucleotides is then added to the corresponding 5' end (see Figure 1d, and Blackburn and Collins (2011); Rhyu (1995) for further details).

As discussed by Teixeira et al. (2004), the number of nucleotides added varies between a few to more than 100 nucleotides, and the number added is independent of the length of the strand. However, when telomerase is present, it does not lengthen every telomere in the cell (as illustrated in Figure 1d); it is more likely to elongate a telomere if the telomere is shorter (we also refer the reader to the telomere length regulation model proposed by Greider (2016)).

At the end of the replication process (including the possible elongation of telomeres), the resulting object is a pair of chromosomes, with one being given to each of the daughter cells.

## 2.2 Branching model

We now propose a stochastic model for the evolution of a population of cells and their chromosomes based on the biological mechanism presented in the previous section. In particular, we will define a branching process that mimics the behaviour of a population of dividing cells as a function of their telomeres.

We first introduce some notation in order to keep track of the telomere lengths in each chromosome. Recall that the 3' end overhangs the opposing 5' end of a DNA strand. We denote by $\omega$ the length of this overhang and assume that it is the same for all telomeres. We will also write $(m, n)$ to denote a chromosome whose 3' ends are of length $m$ and $n$ (classical units are either nucleotides (nt) for single strand regions or base pair (bp) and kilobase (kb) for double stranded regions, and, for sake of simplicity we will use kilobases throughout the paper). Using this notation, it follows that at one end of the chromosome the 3' end and its opposing 5' end are of length $m$ and $m - \omega$, respectively, while at the other end of the chromosome, they are of length $n$ and $n - \omega$, see Figure 1a.

If telomerase is not expressed then, after replication, the chromosome $(m, n)$ gives rise to two daughter chromosomes $(m, n - \omega)$ and $(m - \omega, n)$ (see Figure 1c). However, when telomerase is present, each of the four new telomeres can be elongated. Recall that if the 3' end of one strand is elongated then the complementary 5' end will also be elongated, keeping an overhang of size $\omega$. We will assume that the four telomeres (two per daughter chromosome) are elongated independently of the other telomeres, and with a probability that depends on its length. Moreover, the size of the elongation is assumed to be independent of the length of the telomere.

To represent this mathematically, let $q(\ell)$ denote the probability that a telomere with length $\ell$ will be elongated, and let $\chi_\ell$ denote a Bernoulli random variable with parameter $q(\ell)$. Then, given $\chi_\ell = 1$, we let $L$ denote the (random) quantity added to the telomere, whose law is given by a probability distribution $\mu_f$ on $\mathbb{N} = \{0, 1, 2, \dots\}$. Writing $\Lambda_\ell = \chi_\ell L$, the two resulting daughter chromosomes are thus given by $(m + \Lambda_m, n - \omega + \Lambda'_{n-\omega})$ and
<figure><img src="image_2.png" /><figcaption>(a) Chromosome before cell division.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Replication bubble and replication forks</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) The two daughter chromosomes before lengthening by telomerase</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(d) The two daughter chromosomes after lengthening by telomerase</figcaption></figure>

Figure 1: Complete reproduction of a chromosome with the new overhangs at each end.
$(m - \omega + \Lambda''_{m-\omega}, n + \Lambda'''_{n})$ (see Figure 1d), where $(\Lambda_{\ell})_{\ell \in \mathbb{N}}$, $(\Lambda'_{\ell})_{\ell \in \mathbb{N}}$, $(\Lambda''_{\ell})_{\ell \in \mathbb{N}}$ and $(\Lambda'''_{\ell})_{\ell \in \mathbb{N}}$ are mutually independent and identically distributed. In the example of Figure 1d, $\Lambda'''_{n} = 0$.

*Remark 2.1.* According to the empirical evidence presented by Teixeira et al. (2004), we choose the distribution $\mu_f$ to be independent of the length of the telomere $\ell$. However, the mathematical results and proofs developed in the next section can be easily adapted to the more general setting where the distribution does depend on $\ell$.

Apart from their telomere lengths, cells are also represented by their state: the cell can be in an active or a non-active state. In an active state, it can undergo replication, deactivation (i.e. cellular senescence, defined as an irreversible arrest of cell proliferation (Hayat, 2017)), or be removed from the system (including in particular apoptosis). In a non-active state (i.e. after senescence), the cell cannot replicate, and can thus only be removed from the system.

With this in mind, we now define our model at the level of the cells. Let $N_t$ be the number of cells in the system at time $t \ge 0$. The collection of cells will be denoted by $\{(c_i(t), x_i(t)) : i = 1, \dots, N_t\}$, where $c_i(t)$ represents the lengths of the telomeres in the $i$-th cell at time $t$ and $x_i(t)$ represents its state (active or non-active). More precisely, letting $K$ denote the number of chromosomes in a cell, we write $c_i(t) = (m_{i,j}(t), n_{i,j}(t))_{j=1}^K \in (\mathbb{N} \times \mathbb{N})^K$, where $(m_{i,j}(t), n_{i,j}(t))$ represents the $j^{th}$ chromosome in cell $i$ at time $t \ge 0$. Moreover, if the $i^{th}$ cell is active we set $x_i(t) = \mathbf{a}$, and otherwise $x_i(t) = \mathbf{na}$. Note that $N_t$ also counts the number of non-active cells. We assume that if the minimum telomere length of a cell is smaller than a given value $L_{\min} \ge \omega$, then it is inactive.

For a given element $c = (m_j, n_j)_{j=1}^K \in (\mathbb{N} \times \mathbb{N})^K$, we let $\min c := \min\{m_j, n_j, 1 \le j \le K\}$ denote the minimum telomere length in the cell and $\max c := \max\{m_j, n_j, 1 \le j \le K\}$ the maximum telomere length in the cell. We thus denote by $E = E_{\mathbf{a}} \cup E_{\mathbf{na}}$ the set of possible ‘values’ a cell can take, with $E_{\mathbf{a}}$ and $E_{\mathbf{na}}$ denoting the set of values of the active cells and non-active cells, respectively, i.e.

$$E_{\mathbf{a}} := \{c \in (\mathbb{N} \times \mathbb{N})^K : \min c \ge L_{\min}\} \times \{\mathbf{a}\} \text{ and } E_{\mathbf{na}} := (\mathbb{N} \times \mathbb{N})^K \times \{\mathbf{na}\}.$$

Hence, a cell is represented by an element of the form $(c, x) \in E$.

Let

$$X_t := \sum_{i=1}^{N_t} \delta_{(c_i(t), x_i(t))}, \quad t \ge 0,$$

denote the population of cells at time $t$, which is an element of $\mathcal{M}(E)$, the set of non-negative discrete measures on $E$. The system evolves as a branching process so that, given their point of creation, cells evolve independently of each other according to the following dynamics. In what follows, and throughout the rest of the article, we use the term *rate* for a function $r_x(s, c)$, $x \in \{\mathbf{na}, \mathbf{a}\}$, $s \ge 0$, $c \in E$, to mean there is a random time, $T$, whose distribution is given by $\text{Pr}_{u,(c,x)}(T > t) = e^{-\int_u^t r_x(s,c)ds}$, $0 \le u \le t$.

An active cell $(c, \mathbf{a}) = ((m_j, n_j)_{j=1}^K, \mathbf{a}) \in E_{\mathbf{a}}$ will remain in the same state until one of the following events occur.
• At rate $d_a(t, c)$, at time $t \ge 0$, the cell is removed from the system:

$$X_t = X_{t-} - \delta_{(c,a)}.$$

Biologically, this represents the event that the cell is removed from the system via e.g. apoptosis.

• At rate $s_a(t, c)$, at time $t \ge 0$, the cell becomes non-active, which means that $x$ switches from a to na:

$$X_t = X_{t-} - \delta_{(c,a)} + \delta_{(c,na)},$$

which represents the physical phenomena of cellular senescence.

• at rate $b_a(t, c)$, at time $t \ge 0$, the cell divides into two daughter cells given by,

$$c^1 := c^1(c) = \{(m_j - \omega B_j + \Lambda_{m_j-\omega B_j,j}, n_j - \omega(1-B_j) + \Lambda'_{n_j-\omega(1-B_j),j})^K_{j=1}\}$$

and

$$c^2 := c^2(c) = \{(m_j - \omega(1 - B_j) + \Lambda''_{m_j-\omega(1-B_j),j}, n_j - \omega B_j + \Lambda'''_{n_j-\omega B_j,j})^K_{j=1}\},$$

where $B_j, j = 1, \dots, K$ are a collection of independent Bernoulli random variables with parameter 1/2, describing the allocation of daughter chromosomes in both daughter cells, and the $\Lambda_{\ell,j}, \Lambda'_{\ell,j}, \Lambda''_{\ell,j}$ and $\Lambda'''_{\ell,j}$ are independent (of each other and the $B_j$), with law described in the previous subsection. If $\min c^i < L_{\min}$, for $i \in \{1, 2\}$, then the cell becomes non-active immediately, and if $\min c^i \ge L_{\min}$, we assume that the cell is, as its mother cell, active. Hence, the term $\delta_{(c,a)}$ of $X_{t-}$ is replaced by

$$1_{\min c^1 \ge L_{\min}} \times \delta_{(c^1,a)} + 1_{\min c^1 < L_{\min}} \times \delta_{(c^1,na)} + 1_{\min c^2 \ge L_{\min}} \times \delta_{(c^2,a)} + 1_{\min c^2 < L_{\min}} \times \delta_{(c^2,na)}.$$

Physically, this describes the event of cell replication and the possible configurations of the daughter cells.

A non-active cell $(c, na) \in E_{na}$ remains in the same state until it is removed from the system, at rate $d_{na}(t, c)$, at time $t \ge 0$:

$$X_t = X_{t-} - \delta_{(c,na)}.$$

Again, this is the mathematical representation of apoptosis or other mechanisms that may cause cells to 'die'.

# 3 Mean semigroup of the branching process

## 3.1 Evolution equation of the mean semigroup

In this section, we consider the dynamics of the expectation of the branching process $X_t$. We thus define the expectation semigroup associated to this process. For $(c, x) \in E$ and a bounded measurable function $g: E \to \mathbb{R}$, set

$$\psi_{u,t}[g](c, x) := \mathbb{E}_{u,(c,x)}[\langle g, X_t \rangle] := \mathbb{E}_{u,(c,x)}\left[\sum_{i=1}^{N_t} g(c_i(t), x_i(t))\right], \quad t \ge u \ge 0, \quad (1)$$
where $\mathbb{E}_{u,(c,x)}$ is the expectation operator associated to the law $\mathbb{P}_{u,(c,x)}(\cdot) := \mathbb{P}(\cdot | X_u = \delta_{(c,x)})$, i.e. the law of the process whose initial population is composed of a single cell, $(c, x) \in E$ at time $u$.

In order to understand the evolution of the mean semigroup $(\psi_{u,t})_{0\le u\le t}$, we introduce the following assumption that ensures only a finite number of events can occur in finite time intervals.

**Assumption 3.1.** The birth, deactivation and removal rates, $b_a$, $s_a$, $d_a$ and $d_{na}$, are bounded.

Based on the dynamics described in the previous section, we have the following proposition, which gives the evolution equation for $(\psi_{u,t})_{0\le u\le t}$. In what follows, we say that the semigroup $\psi$ is bounded in time if $\|\psi_{u,t}[g]\|_{\infty} < \infty$ for all positive bounded measurable functions $g$ and for all $0 \le u \le t$.

**Proposition 3.2.** Let $g: E \to \mathbb{R}$ be a bounded measurable function, $0 \le u \le t$, $(c, x) = ((m_j, n_j)_{j=1}^K, x) \in E$. Under Assumption 3.1, the expectation semigroup $(\psi_{u,t})_{0\le u\le t}$ is the unique solution, that is bounded in time, to

$$ \psi_{u,t}[g](c, x) = g(c, x) + \mathbf{1}_{\{x=a\}} \int_{u}^{t} b_{a}(v, c) [\mathcal{F}[\psi_{v,t}[g]](c, a) - \psi_{v,t}[g](c, a)] dv \\ + \mathbf{1}_{\{x=a\}} \int_{u}^{t} s_{a}(v, c) [\psi_{v,t}[g](c, na) - \psi_{v,t}[g](c, a)] - \int_{u}^{t} d_{x}(v, c) \psi_{v,t}[g](c, x) dv, \quad (2) $$

where the operator $\mathcal{F}$ is defined by

$$ \mathcal{F}[g](c, a) = \mathbb{E} \left[ \mathbf{1}_{\min c^1 \ge L_{\min}} g(c^1, a) + \mathbf{1}_{\min c^1 < L_{\min}} g(c^1, na) \right. \\ \left. \qquad + \mathbf{1}_{\min c^2 \ge L_{\min}} g(c^2, a) + \mathbf{1}_{\min c^2 < L_{\min}} g(c^2, na) \right], \quad (3) $$

with

$$ c^1 = (m_j - \omega B_j + \Lambda_{m_j - \omega B_j, j}, n_j - \omega(1 - B_j) + \Lambda'_{n_j - \omega(1 - B_j), j})_{j=1}^K, $$

$$ c^2 = (m_j - \omega(1 - B_j) + \Lambda''_{m_j - \omega(1 - B_j), j}, n_j - \omega B_j + \Lambda'''_{n_j - \omega B_j, j})_{j=1}^K, $$

and where the expectation $\mathbb{E}$ is taken with respect to the law of the $B_j$, $\Lambda_{\ell,j}$, $\Lambda'_{\ell,j}$, $\Lambda''_{\ell,j}$ and $\Lambda'''_{\ell,j}$.

*Proof.* The proof is similar to (Cox et al., 2019, §6). Let us recall the main arguments. By Assumption 3.1, $\psi$ is bounded in time. By conditioning $\psi_{u,t+t'}$ at time $t$ and applying the Markov property, it is a straightforward exercise to show that $(\psi_{u,t})_{t\ge u\ge 0}$ is a semigroup. To show that it satisfies (2) first consider the case when $x = a$. Splitting the expectation semigroup on the first event (branching, deactivation or removal), we have

$$ \psi_{u,t}[g](c, a) = g(c, a)e^{-\int_{u}^{t} b_{a}(v,c) + d_{a}(v,c) + s_{a}(v,c)dv} \\ + \int_{u}^{t} b_{a}(v,c)e^{-\int_{u}^{v} b_{a}(w,c) + d_{a}(w,c) + s_{a}(w,c)dw} \mathcal{F}[\psi_{v,t}[g]](c, a)dv \\ + \int_{u}^{t} s_{a}(v,c)e^{-\int_{u}^{v} b_{a}(w,c) + d_{a}(w,c) + s_{a}(w,c)dw} \psi_{v,t}[g](c, na)dv, \quad (4) $$
where the second term follows from an application of the branching property and the strong Markov property. Applying A.1, we obtain equation (2) for the case $x = a$.

For the case where $x = na$, the cell can only remain non-active or be removed from the system. Thus splitting $\psi_{u,t}[g](c, na)$ on the first time the cell is removed from the system and again, applying A.1, we obtain

$$ \psi_{u,t}[g](c, na) = g(c, na) - \int_{u}^{t} d_{na}(v, c) \psi_{v,t}[g](c, na) dv, \quad (5) $$

that is equation (2) for the case $x = na$.

We now show that (2) admits a unique solution. Suppose that we have two solutions bounded in time, $\psi^{(1)}$ and $\psi^{(2)}$, to (2). Then, denoting $\chi_{u,t}[g] = \psi_{u,t}^{(1)}[g] - \psi_{u,t}^{(2)}[g]$ and $\bar{\chi}_{u,t}[g] = \sup_{(c,x)} |\chi_{u,t}[g](c, x)|$, due to Assumption 3.1, we have

$$ \begin{aligned} |\chi_{u,t}[g](c, x)| &= |\mathbf{1}_{\{x=a\}} \int_{u}^{t} b_{a}(v, c) [\mathcal{F}[\chi_{v,t}[g]](c, a) - \chi_{v,t}[g](c, a)] dv \\ &\quad + \mathbf{1}_{\{x=a\}} \int_{u}^{t} s_{a}(v, c) [\chi_{v,t}[g](c, na) - \chi_{v,t}[g](c, a)] - \int_{u}^{t} d_{x}(v, c) \chi_{v,t}[g](c, x) dv| \\ &\le (3\|b_{a}\|_{\infty} + 2\|s_{a}\|_{\infty} + \|d_{a}\|_{\infty} + \|d_{na}\|_{\infty}) \int_{u}^{t} \bar{\chi}_{v,t}[g] dv. \end{aligned} $$

Fixing $t$ and setting $f(r) = \bar{\chi}_{t-r,t}[g]$, we obtain, for all $r \in [0, t]$

$$ \begin{aligned} f(r) &\le (3\|b_a\|_{\infty} + 2\|s_a\|_{\infty} + \|d_a\|_{\infty} + \|d_{na}\|_{\infty}) \int_{t-r}^{t} \bar{\chi}_{v,t}[g] dv \\ &= (3\|b_a\|_{\infty} + 2\|s_a\|_{\infty} + \|d_a\|_{\infty} + \|d_{na}\|_{\infty}) \int_{0}^{r} \bar{\chi}_{t-w,t}[g] dw \\ &= (3\|b_a\|_{\infty} + 2\|s_a\|_{\infty} + \|d_a\|_{\infty} + \|d_{na}\|_{\infty}) \int_{0}^{r} f(w) dw. \end{aligned} $$

Uniqueness now follows easily from Grönwall's inequality and from $f(0) = \bar{\chi}_{t,t}[g] = 0$. $\square$

## 3.2 Many-to-one

We now consider a many-to-one representation for the semigroup $(\psi_{u,t})_{t \ge u \ge 0}$. The many-to-one formula offers another representation for the first moment of the system of cells via a weighted jump process. The motivation for considering this second representation is two-fold. Firstly, as we shall see in the next section, it more easily allows us to characterise the long-term behaviour of the branching process: indeed, considering the long-term behaviour of a single (weighted) trajectory is much more tractable than that of the entire system of cells. Secondly, we will later simulate our model using interacting particle systems that are based on multiple copies of single trajectories.

To this end, consider the process $(C_t, \mathcal{X}_t)_{t \ge 0}$ that evolves as a pure jump process in $E \cup \{\partial\}$, where $(C_t, \mathcal{X}_t) = \partial$ means that the cell is removed from the system and does not evolve further (in particular, the point $\partial$ is absorbing). The dynamics of this process are as follows.

* When in configuration $(c, a) = ((m_j, n_j)_{j=1}^K, a) \in E_a$, the process
- jumps with rate $2b_a(t, c)$, at time $t$, to $(c', x')$ with

$$c' = (m_j - \omega B_j + \Lambda_{m_j - \omega B_j, j}, n_j - \omega(1 - B_j) + \Lambda'_{n_j - \omega(1 - B_j), j})_{j=1}^K,$$

$$x' = \begin{cases} \mathbf{a} & \text{if } \min c' \ge L_{\min}, \\ \mathbf{na} & \text{if } \min c' < L_{\min}; \end{cases}$$

- switches with rate $s_a(t, c)$ from $(c, \mathbf{a})$ to $(c, \mathbf{na})$;

- switches with rate $d_a(t, c)$ from $(c, \mathbf{a})$ to $\partial$.

• The process in a configuration $(c, \mathbf{na})$ switches to $\partial$ with rate $d_{\mathbf{na}}(t, c)$.

Note that for active cells, when a jump occurs at rate $2b_a$, the cell jumps from $C_{t-}$ to the cell $c^1$ defined as in (3). Note that as the law of $c^1$ and $c^2$ are the same, it is equivalent to choose the cell $c^1$ with probability $1/2$ and the cell $c^2$ with probability $1/2$. Then, if $\min c' \ge L_{\min}$, the cell remains active ($\mathcal{X}_t = \mathbf{a}$), and otherwise the cell becomes inactive ($\mathcal{X}_t = \mathbf{na}$).

Let $\mathbf{P}_{u,(c,x)}$ denote the law of this process when started from a single cell with configuration $(c, x)$ at time $u$, with corresponding expectation operator $\mathbf{E}_{u,(c,x)}$. Similarly to the proof of Proposition 3.2, it follows that the linear semigroup $(\varphi_{u,t})_{t \ge u \ge 0}$ associated to $(C_t, \mathcal{X}_t)_{t \ge 0}$ satisfies

$$\begin{aligned} \varphi_{u,t}[g](c, x) &:= \mathbf{E}_{u,(c,x)} [g(C_t, \mathcal{X}_t) \mathbf{1}_{t<\tau}] \\ &= g(c, x) + \mathbf{1}_{x=\mathbf{a}} \int_u^t 2b_{\mathbf{a}}(v, c) \left[ \frac{1}{2} \mathcal{F}[\varphi_{v,t}[g]](c, \mathbf{a}) - \varphi_{v,t}[g](c, \mathbf{a}) \right] dv \\ &\quad + \mathbf{1}_{x=\mathbf{a}} \int_u^t s_{\mathbf{a}}(v, c) [\varphi_{v,t}[g](c, \mathbf{na}) - \varphi_{v,t}[g](c, \mathbf{a})] dv \\ &\quad - \int_u^t d_x(v, c) \varphi_{v,t}[g](c, x) dv, \end{aligned} \quad (6)$$

where $(c, x) \in E$, $g: E \to \mathbb{R}$ is a measurable, bounded function, $0 \le u \le t$ and $\tau$ denotes the hitting time of $\partial$ by the process $(C, \mathcal{X})$.

The above semigroup describes the average behaviour of a typical particle in the branching process, $X$. However, this semigroup does not take into account mass creation and loss, as in the branching process. The following many-to-one formula shows one how to deal with this.

**Lemma 3.3.** Under the assumptions of Proposition 3.2, the semigroup, defined for all bounded measurable functions $g: E \to \mathbb{R}$, initial times $u$, initial cell configurations $(c, x) \in E$ and times $t \ge u$ by

$$\phi_{u,t}[g](c, x) := \mathbf{E}_{u,(c,x)} \left[ \exp \left( \int_u^t \mathbf{1}_{\{\mathcal{X}_v=\mathbf{a}\}} b_{\mathbf{a}}(v, C_v) dv \right) g(C_t, \mathcal{X}_t) \mathbf{1}_{t<\tau} \right], \quad (7)$$

also solves equation (2), and hence $\psi_{u,t} = \phi_{u,t}$ for each $t \ge u \ge 0$.
*Sketch proof.* The method behind the proof is almost identical to that of Proposition 3.2 and so we outline the steps and leave the details as an exercise for the reader.

Conditioning on the first event time, as outlined in the dynamics described above, yields an equation similar to that of (6) except that the first term will be weighted by $e\int_u^t 1_{\{x=a\}} b_a(s,c)ds$, and each of the integrands will be weighted by $e\int_v^t 1_{\{x=a\}} b_a(s,c)ds$. The conclusion now follows from an application of Theorem A.1. □

The main idea behind the many-to-one lemma is that since we are only interested in the average behaviour of the branching process, rather than the pathwise behaviour, we can interpret the first moment of the whole branching particle system as the average behaviour of a typical trajectory in the whole particle system, weighted by the number of particles we expect to see in the system. In this case, the process $(C_t, \mathcal{X}_t)_{t\ge0}$ follows the dynamics of a 'typical trajectory' and the expectation of the exponential term in (7) corresponds to the expected number of particles in the branching process at time $t$.

# 4 Long-term behaviour

This section is devoted to studying the asymptotic stability of the telomere length profile, as stated in Theorem 4.5 below. For simplicity, we only consider the time homogeneous dynamics (hence replacing $b_a(t, c)$ by $b_a(c)$, $d_a(t, c)$ by $d_a(c)$, $\psi_{s,t}$ by $\psi_{t-s}$ and so on). This will allow us to make a connection to the theory of quasi-stationary distributions. We refer the reader to the book of Del Moral (2004) and the articles of Del Moral and Miclo (2002); Del Moral (2013); Bansaye et al. (2019); Champagnat and Villemonais (2018); Velleret (2019) for time inhomogeneous versions of this theory.

In order to state our main result, we need some additional assumptions. We first make a technical assumption (Assumption 4.1) on the telomere lengthening caused by telomerase, which will ensure that the process is irreducible in $E_a$ (see Lemma 4.6 below). The second assumption (Assumption 4.2) then imposes that the process is supercritical (in the usual sense, made precise below). Finally we make assumptions which guarantee quasi-compactness (Assumptions 4.3 and 4.4). We also remind the reader that the table of notation included in Appendix C may be of use.

**Assumption 4.1.** We assume that

1. there exists $b_{tel} > \omega$ in $\mathbb{N}$ such that $\mu_f\{b_{tel}\} > 0$ and $\mu_f\{b_{tel} - 1\} > 0$,

2. $b_a(c) > 0$ for all $(c, \mathbf{a}) \in E_a$,

3. $0 < q(i) < 1$ for all $i \ge L_{\min}$.

We will soon see that under Assumptions 3.1 and 4.1, the process is irreducible on $E_a$ (see Lemma 4.6) and, by construction, its complement is absorbing. Hence we can define

$$ \lambda_0 = \inf \left\{ \lambda \in \mathbb{R}, \text{ such that } \liminf_{t \to \infty} e^{-\lambda t} \psi_t[\mathbf{1}_F](c, \mathbf{a}) < +\infty \right\}, \quad (8) $$
which does not depend on $(c, a) \in E_a$ nor on the (arbitrary) non-empty finite set $F \subset E_a$. The informal interpretation of $\lambda_0$ is that it measures the asymptotic exponential growth of the expected number of active cells in the population. We say that the process is *supercritical* when $\lambda_0 > 0$, *subcritical* when $\lambda_0 < 0$ and *critical* when $\lambda_0 = 0$. The value $\lambda_0$ is sometimes referred to as the Malthusian parameter.

In the rest of this section, we focus on the supercritical case, which corresponds to exponential growth of the average number of active cells (such as expanding cancerous tumour (Blasco, 2005) or germinal cells (Hiyama and Hiyama, 2007, Fig. 1)).

**Assumption 4.2.** We have $\lambda_0 > 0$.

Next, we introduce an assumption on the moments of $\mu_f$, the law of the length added to telomeres when telomerase acts.

**Assumption 4.3.** We assume that there exists $\alpha_0 > \frac{1}{\omega} \ln \frac{\|b_a\|_{\infty}}{\lambda_0}$ such that $\mu_f$ admits an exponential moment of order $\alpha_0 > 0$, i.e.

$$ \sum_{n \ge 0} \mu_f\{n\} \exp(\alpha_0 n) < +\infty. $$

We conclude with an assumption on the probability that telomerase acts on a telomere, depending on its length.

**Assumption 4.4.** We assume that the probability of telomerase activity, $q(i)$, decreases to 0 when $i \to +\infty$.

Note that Assumption 4.4 is satisfied, in particular, when $q(i)$ decreases geometrically fast in $i$, as suggested in the telomere length regulation model (Greider, 2016). Indeed, the replication fork model described therein posits that telomerase travels with the lagging strand replication machinery, requiring it to stay bound to the fork through nucleosomes and telomere proteins, either of which can cause dissociation of telomerase, until reaching the chromosome terminus for telomere extension. By linking telomerase activity to fork progression subject to the renewal of succeeding non-dissociation events along the telomere, this model suggests that the activity of telomerase occurs with a probability that decreases geometrically fast in the length of the telomere.

We are now ready to state our main result.

**Theorem 4.5.** Suppose Assumptions 3.1, 4.1, 4.2, 4.3 and 4.4 hold. Then there exists a function $\eta: E \to \mathbb{R}$, positive on $E_a$ and vanishing on $E_{na}$, a probability measure $\nu$ on E, and a function $V: E \to [1, \infty)$ such that, for all $g: E \to \mathbb{R}$ satisfying $|g| \le V$, we have

$$ |e^{-\lambda_0 t} \psi_t[g](c, x) - \eta(c, x) \nu[g]| \le C e^{-\gamma t} V(c, x), \quad t \ge 0, (c, x) \in E, $$

where $C, \gamma$ are positive constants. In addition, one can choose $V$ such that, for some $k \ge 1, \alpha > 0$,

$$ V(c, na) = 1 \quad \text{and} \quad V(c, a) = \exp \left[ \alpha \sum_{j=1}^{K} \left( (m_j - k + \omega)_+ + (n_j - k + \omega)_+ \right) \right]. $$
We prove this Theorem in three steps. The first takes the form of a Lemma that shows that the process is irreducible on $E_a$. The second gives the asymptotic behavior for the process restricted to the active particles (Proposition 4.7). Finally the proof of Theorem 4.5 follows from Theorem 3.1 of Champagnat and Villemonais (2022).

**Lemma 4.6.** Under Assumptions 3.1 and 4.1, the set $E_a$ is irreducible for the process $(C_t, \mathcal{X}_t)_{t \ge 0}$ defined in Section 3.2.

*Proof.* Let $(c, x) = ((m_j, n_j)_{j=1}^K, \mathbf{a}) \in E_a$. Let us denote by $\tau_1 < \tau_2$ the two first jump times of the process (where $\tau_2 = +\infty$ if $(C_{\tau_1}, \mathcal{X}_{\tau_1}) = \partial$). Then, for any fixed $t_0 > 0$ and any $(A_1 \times \cdots \times A_K) \times \{\mathbf{a}\} \in E_a$, we have

$$ 
\begin{aligned}
\mathbf{P}_{(c,x)}((C_{t_0}, \mathcal{X}_{t_0}) \in (A_1 \times \cdots \times A_K) \times \{\mathbf{a}\}) & \\
& \ge \mathbf{P}_{(c,x)}(\tau_1 < t_0 < \tau_2, (C_{t_0}, \mathcal{X}_{t_0}) \in (A_1 \times \cdots \times A_K) \times \{\mathbf{a}\}) \\
& \ge (1 - e^{-t_0 \inf(2b_a+s_a+d_a)})e^{-t_0 \sup(2b_a+s_a+d_a)} \frac{\inf 2b_a}{\sup(2b_a + s_a + d_a)} \\
& \times \mu_{(m_1,n_1)}(A_1) \times \cdots \times \mu_{(m_K,n_K)}(A_K),
\end{aligned} 
\quad (9) $$

where

$$ \mu_{(m_j,n_j)}(A_j) = \mathbb{P}((m_j - \omega B_j + \Lambda_{m_j-\omega B_j,j}, n_j - \omega(1 - B_j) + \Lambda'_{n_j-\omega(1-B_j),j}) \in A_j). $$

Note that due to the second part of Assumption 4.1, $\inf 2b_a > 0$.

We observe that, restricting to the events $\{B_j = 1, \Lambda'_{n_j,j} = 0\}$ (i.e., we have attrition $-\omega$ on the left-hand side of the chromosome and telomerase does not act on the right-hand side) and $\{B_j = 0, \Lambda_{m_j,j} = 0\}$, we obtain

$$ 
\begin{aligned}
\mu_{(m_j,n_j)}(\cdot) & \ge \mathbb{P}(B_j = 1, \Lambda'_{n_j,j} = 0) \mathbb{P}((m_j - \omega + \Lambda_{m_j,j}, n_j) \in \cdot) \\
& \quad + \mathbb{P}(B_j = 0, \Lambda_{m_j,j} = 0) \mathbb{P}((m_j, n_j - \omega + \Lambda'_{n_j,j}) \in \cdot) \\
& \ge C \left( \delta_{(m_j+b_{tel}-\omega,n_j)} + \delta_{(m_j+b_{tel}-\omega-1,n_j)} + \delta_{(m_j-\omega,n_j)} \right) \\
& \quad + C \left( \delta_{(m_j,n_j+b_{tel}-\omega)} + \delta_{(m_j,n_j+b_{tel}-\omega-1)} + \delta_{(m_j,n_j-\omega)} \right),
\end{aligned} $$

where Assumption 4.1 ensures that $C$ can be chosen to be positive.

Now consider the discrete time Markov process $Z^j$ evolving in $\mathbb{N}^2$ and with transition kernel

$$ \frac{1}{6} \left( \delta_{(m_j+b_{tel}-\omega,n_j)} + \delta_{(m_j+b_{tel}-\omega-1,n_j)} + \delta_{(m_j-\omega,n_j)} + \delta_{(m_j,n_j+b_{tel}-\omega)} + \delta_{(m_j,n_j+b_{tel}-\omega-1)} + \delta_{(m_j,n_j-\omega)} \right). $$

This process can jump with positive probability from $(m_j, n_j)$ to $(m_j + \omega(b_{tel} - \omega), n_j)$ in $\omega$ steps and from $(m_j + \omega(b_{tel} - \omega), n_j)$ to $(m_j, n_j)$ in $b_{tel} - \omega$ steps, so that there is path of length $b_{tel}$ linking $(m_j, n_j)$ to itself. Similarly, there exists a path of length $b_{tel} - 1$ linking $(m_j, n_j)$ to itself. This implies that the process is aperiodic. Moreover, the process can jump from $(m_j, n_j)$ to $(m_j + \omega(b_{tel} - \omega - 1) + 1, n_j)$ in $\omega$ steps, and then come back to $(m_j + 1, n_j)$
in $b_{tel} - \omega - 1$ steps. This shows that the process can reach $(m_j + 1, n_j)$ from $(m_j, n_j)$. Similarly, one shows that the process can reach $(m_j - 1, n_j)$, $(m_j, n_j + 1)$, $(m_j, n_j - 1)$ from $(m_j, n_j)$ in less than $b_{tel}$ steps. In particular, we deduce that, for all $(m_1, n_1), \dots, (m_K, n_K)$ and $(m'_1, n'_1), \dots, (m'_K, n'_K)$, there exists $n_0$ such that

$$ \mathbb{P}(Z_{n_0}^1 = (m'_1, n'_1), \dots, Z_{n_0}^K = (m'_K, n'_K) \mid Z_0^1 = (m_1, n_1), \dots, Z_0^K = (m_K, n_K)) > 0, $$

where the $Z^i$ are chosen independent. Using the definition of $\mu_{(m_j, n_j)}$ and the inequality (9), we deduce that, for all $(c', \mathbf{a}) \in E_{\mathbf{a}}$,

$$ \mathbf{P}_{(c,x)}(C_{n_0 t_0} = c', \mathcal{X}_{n_0 t_0} = \mathbf{a}) > 0. $$

□

In the next proposition, we use the notation

$$ G := \{c \in (\mathbb{N} \times \mathbb{N})^K, \min c > L_{min}\} $$

to denote the set of possible values of active cells.

**Proposition 4.7.** Suppose Assumptions 3.1, 4.1, 4.2, 4.3 and 4.4 hold. Then there exists a positive function $\eta_{\mathbf{a}}: G \to \mathbb{R}$, a probability measure $\nu_{\mathbf{a}}$ on $G$, and a function $V_{\mathbf{a}}: G \to [1, \infty)$ such that, for all $g: G \to \mathbb{R}$ satisfying $|g| \le V_{\mathbf{a}}$,

$$ |e^{-\lambda_0 t} \psi_t[g_0](c, \mathbf{a}) - \eta_{\mathbf{a}}(c) \nu_{\mathbf{a}}[g]| \le C e^{-\gamma t} V_{\mathbf{a}}(c), \quad \forall t \ge 0, \forall c \in G, $$

where $g_0(c, x) := g(c)\mathbf{1}_{c \in G, x=a}$, and where $C, \gamma$ are positive constants. In addition, one can choose $V_{\mathbf{a}}$ such that, for some $k \ge 1, \alpha > 0$,

$$ V_{\mathbf{a}}(c) = \exp \left[ \alpha \sum_{j=1}^{K} \left( (m_j - k + \omega)_+ + (n_j - k + \omega)_+ \right) \right]. $$

**Remark 4.8.** The above proposition is stated for $\lambda_0 > 0$ which corresponds to exponentially growing populations, which is our main focus. However, a straightforward adaptation of the proof shows that the result also holds if $\inf(d_{\mathbf{a}} + s_{\mathbf{a}} + \lambda_0) > 0$. In addition, numerical simulations suggest that this result may hold in a more general context, at least for some choices of the model parameters and without restriction on $\lambda_0$ (see Remark B.1).

**Proof.** In order to prove the proposition, we consider the semigroup on $L^\infty(G)$ (the space of bounded measurable functions on $G$), defined, for all $c \in G$ and all bounded measurable function $g: G \to \mathbb{R}$, by

$$ \begin{aligned} \tilde{\varphi}_t[g](c) &:= e^{-\|b_{\mathbf{a}}\|_{\infty} t} \phi_t[g_0](c, \mathbf{a}) \\ &= \mathbf{E}_{(c, \mathbf{a})} \left[ \exp \left( -\int_0^t (\|b_{\mathbf{a}}\|_{\infty} - b_{\mathbf{a}}(C_v)) dv \right) g_0(C_t, \mathbf{a}) \mathbf{1}_{t<\tau} \right] \\ &=: \mathbf{E}_c \left[ g(\tilde{C}_t) \mathbf{1}_{t<\tilde{\tau}} \right], \end{aligned} \quad (10) $$
where $(\tilde{C}_t)_{t \ge 0}$ is a pure jump sub-Markov process on $G$ evolving as $\mathcal{C}$ under $\mathbf{E}_{c,a}$ but sent to a cemetery point $\dagger \notin G$ at time

$$ \tilde{\tau} := \inf \{t > 0 : \mathcal{C}_t \in (E_{\mathbf{a}})^c\} \wedge \inf \{t > 0 : \int_0^t (\|b_{\mathbf{a}}\|_{\infty} - b_{\mathbf{a}}(C_v))dv > \mathbf{e}\}, $$

where $\mathbf{e}$ is an independent rate 1 exponential random variable (i.e. $(\tilde{C}_t)_{t \ge 0}$ is sent to $\dagger$ from when it reaches $(E_{\mathbf{a}})^c = E_{na} \cup \{\partial\}$ and at an additional rate $\|b_{\mathbf{a}}\|_{\infty} - b_{\mathbf{a}}$). More formally, its infinitesimal generator for bounded measurable functions $g: G \to \mathbb{R}$ is given by

$$ \begin{aligned} \mathcal{L}g(c) &= 2b_{\mathbf{a}}(c) \left( \frac{1}{2} \mathcal{F}[g_0](c, \mathbf{a}) - g(c) \right) - (d_{\mathbf{a}}(c) + s_{\mathbf{a}}(c) + \|b_{\mathbf{a}}\|_{\infty} - b_{\mathbf{a}}(c))g(c) \\ &= b_{\mathbf{a}}(c) \mathcal{F}[g_0](c, \mathbf{a}) - (d_{\mathbf{a}}(c) + s_{\mathbf{a}}(c) + \|b_{\mathbf{a}}\|_{\infty} + b_{\mathbf{a}}(c))g(c), \end{aligned} \quad (11) $$

where $\mathcal{F}$, defined in (3), satisfies

$$ \mathcal{F}[g_0](c, \mathbf{a}) = 2\mathbb{E} \left[ \mathbf{1}_{\min c^1 \ge L_{\min}} g \left( (m_j - \omega B_j + \Lambda_{m_j - \omega B_j, j}, n_j - \omega(1 - B_j) + \Lambda'_{n_j - \omega(1 - B_j), j})^{K}_{j=1} \right) \right]. $$

The irreducibility of $\tilde{\mathcal{C}}$ is a direct consequence of the irreducibility of $\mathcal{C}$ proved in Lemma 4.6. In particular, we can define

$$ \tilde{\lambda}_0 = \inf \left\{ \lambda \in \mathbb{R}, \text{ such that } \liminf_{t \to \infty} e^{-\lambda t} \tilde{\varphi}_t[\mathbf{1}_H](c) < +\infty \right\}, $$

independently of $c \in G$ and of the finite set $H \subset G$. We observe, by Lemma 3.3, (8) and (10), that

$$ \tilde{\lambda}_0 = \lambda_0 - \|b_{\mathbf{a}}\|_{\infty}. \quad (12) $$

Our aim is to apply Theorem 5.1 of Champagnat and Villemonais (2023). In order to do so, it is sufficient to find a Lyapunov type function, $V_{\mathbf{a}}: G \to [1, +\infty)$, such that

$$ \mathcal{L}V_{\mathbf{a}} \le (\tilde{\lambda}_0 - \varepsilon)V_{\mathbf{a}} + C\mathbf{1}_H, \quad (13) $$

for some $\varepsilon > 0$ and finite $H \subset G$. Beware that the convention for the sign of $\lambda_0$ is not the same in the reference. Indeed, the constant $\lambda_0$ therein refers to a *death rate*, while in our case it refers to a *growth rate*. Hence, one need to replace $\lambda_0$ by $-\lambda_0$ in Theorem 5.1 of Champagnat and Villemonais (2023) to reconcile with our setting.

In order to find $V_{\mathbf{a}}$, we set $\alpha \in (\frac{1}{\omega} \ln \frac{\|b_{\mathbf{a}}\|_{\infty}}{\lambda_0}, \alpha_0)$, where $\alpha_0$ is given by Assumption 4.3. We look at functionals of the type

$$ V_k(c) := \exp \left[ \alpha \sum_{j=1}^{K} \left( (m_j - k + \omega)_+ + (n_j - k + \omega)_+ \right) \right], \quad k \ge 0, $$

with $c = (m_j, n_j)_{j=1}^K \in G$. We introduce the event

$$ A_k(c) = \left\{ \begin{array}{l} \Lambda_{m_j - \omega B_j, j} = 0, \forall j \text{ s.t. } m_j > k/2 \\ \cap \left\{ \Lambda'_{n_j - \omega(1 - B_j), j} = 0, \forall j \text{ s.t. } n_j > k/2 \right\} \\ \cap \left\{ \Lambda_{m_j - \omega B_j, j} \lor \Lambda'_{n_j - \omega(1 - B_j), j} < k/2 - \omega, \forall j \right\} \end{array} \right\}. $$
and its probability

$$p_k(c) = \mathbb{P}(A_k(c)).$$

It follows from Assumption 4.4 that for all $c \in E$ the probability of the complementary of $A_k(c)$ satisfies

$$ \mathbb{P}(A_k(c)^c) \le \sum_{j=1}^{K} [q(m_j - \omega B_j)1_{m_j > k/2} + q(n_j - \omega(1 - B_j))1_{n_j > k/2} + 2\mathbb{P}(L \ge k/2 - \omega)] \\ \xrightarrow{k \to +\infty} 0, $$

then

$$ \inf_{c \in E} p_k(c) \xrightarrow{k \to +\infty} 1. $$

Denoting by $c' = (m_j - \omega B_j + \Lambda_{m_j - \omega B_j, j}, n_j - \omega(1 - B_j) + \Lambda'_{n_j - \omega(1 - B_j), j})_{j=1}^K$ the configuration of the process after a jump that doesn't lead to the removal of the particle from the system, we easily check that

$$ V_k(c') \le \exp \left( \alpha \sum_j (\Lambda_{m_j - \omega B_j, j} + \Lambda'_{n_j - \omega(1 - B_j), j}) \right) V_k(c). \quad (14) $$

On the one hand, by Hölder's inequality, for any fixed $1 < p < \frac{\alpha_0}{\alpha}$, we have

$$ \begin{aligned} & \mathbb{E}(\mathbf{1}_{A_k(c)^c} \exp(\alpha \sum_j (\Lambda_{m_j - \omega B_j, j} + \Lambda'_{n_j - \omega(1 - B_j), j}))) \\ & \le (1 - p_k(c))^{1-1/p} \left( \mathbb{E} \left( \exp(p\alpha \sum_j (\Lambda_{m_j - \omega B_j, j} + \Lambda'_{n_j - \omega(1 - B_j), j})) \right) \right)^{1/p} \\ & \le (1 - p_k(c))^{1-1/p} (\theta_{p\alpha})^{2K/p}, \end{aligned} \quad (15) $$

where $\theta_{p\alpha} := \mathbb{E}(\exp(p\alpha L)) = \sum_{n \ge 0} \mu_f\{n\} \exp(p\alpha n) < +\infty$ by Assumption 4.3. On the other hand, conditioning on the event $A_k(c)$, the left-hand telomere length $m'_j$ of the j-th chromosome of $c'$ satisfies

$$ m'_j - k + \omega = m_j - \omega B_j + \Lambda_{m_j - \omega B_j, j} - k + \omega \le \begin{cases} m_j - k + \omega - \omega B_j & \text{if } m_j > k/2 \\ 0 & \text{if } m_j \le k/2 \end{cases} $$

and similarly for the right-hand side. Then $V_k(c') \le V_k(c)$, and if $c$ is such that $\max c \ge k$, then with probability at least $1/2$ (for at least one chromosome side which leads to $\max c \ge k$), $V_k(c') \le e^{-\alpha\omega} V_k(c)$. We then deduce from (11) that, for all $c \in G$ such that $\max c \ge k$,

$$ \begin{aligned} \mathcal{L}V_k(c) & \le 2b_{\mathbf{a}}(c) \left( \frac{p_k(c)}{2} + \frac{p_k(c)}{2}e^{-\alpha\omega} + (1-p_k(c))^{1-1/p}(\theta_{p\alpha})^{2K/p} \right) V_k(c) \\ & \quad - (d_{\mathbf{a}}(c) + s_{\mathbf{a}}(c) + \|b_{\mathbf{a}}\|_{\infty} + b_{\mathbf{a}}(c))V_k(c) \\ & \le (b_{\mathbf{a}}(c)e^{-\alpha\omega} + 2b_{\mathbf{a}}(c)(1-p_k(c))^{1-1/p}(\theta_{p\alpha})^{2K/p} - d_{\mathbf{a}}(c) - s_{\mathbf{a}}(c) - \|b_{\mathbf{a}}\|_{\infty}) V_k(c) \\ & \le (b_{\mathbf{a}}(c)e^{-\alpha\omega} + 2b_{\mathbf{a}}(c)(1-p_k(c))^{1-1/p}(\theta_{p\alpha})^{2K/p} - \|b_{\mathbf{a}}\|_{\infty}) V_k(c). \end{aligned} $$
As $\alpha$ has been chosen such that $\alpha > \frac{1}{\omega} \ln \frac{\|b_a\|_{\infty}}{\lambda_0}$ then

$$ \varepsilon := \frac{1}{2}(\|b_a\|_{\infty} + \tilde{\lambda}_0 - \|b_a\|_{\infty}e^{-\alpha\omega}) = \frac{1}{2}(\lambda_0 - \|b_a\|_{\infty}e^{-\alpha\omega}) > 0, $$

where we used (12), so that $\|b_a\|_{\infty} - b_a(c)e^{-\alpha\omega} \ge 2\varepsilon - \tilde{\lambda}_0$. Moreover, as $\theta_{p\alpha} < \infty$, then choosing $k$ large enough such that $2b_a(c)(1-p_k(c))^{1-1/p}(\theta_{p\alpha})^{2K/p} \le \varepsilon$ for all $c \in G$, we deduce that, for all $c \in G$ such that $\max c \ge k$,

$$ \mathcal{L}V_k(c) \le (\tilde{\lambda}_0 - \varepsilon)V_k(c). $$

Now, choosing the finite set $H := \{c \in G, \max c \le k\}$, from (14) and (11), $\mathcal{L}V_k \le V_k$ which is bounded in the compact set $H$. We deduce that there exists a constant $C > 0$ such that (13) holds true for $V_a = V_k$, which concludes the proof. $\square$

**Proof of Theorem 4.5.** In a similar manner to (10), consider the pure jump sub-Markov process $(\bar{C}_t, \bar{\chi}_t)_{t \ge 0}$ on $E$ evolving as $(\mathcal{C}, \mathcal{X})$ under $E_{(c,x)}$, but is removed from the system at an additional rate $\|b_a\|_{\infty} - \mathbf{1}_{x=a}b_a(c, x)$. Denote by $E_{(c,x)}$ the expectation operator associated to this process, and by $\bar{\varphi}_t$ the semigroup on $L^{\infty}(E)$ (the space of bounded measurable functions on $E$), defined, for all $(c, x) \in E$ and all bounded measurable functions $g: E \to \mathbb{R}$, by

$$ \bar{\varphi}_t[g](c, x) := e^{-\|b_a t\|_{\infty}} \phi_t[g](c, x) = E_{(c,x)} [g(\bar{C}_t, \bar{\chi}_t) \mathbf{1}_{t<\bar{\tau}}], $$

where $\bar{\tau}$ denotes the time at which the particle is removed from the system $(\bar{C}_t, \bar{\chi}_t)_{t \ge 0}$. Note that, before time $T_{na} := \inf\{t \ge 0, \bar{\chi}_t = \mathbf{na}\}$, the process $\bar{C}_t$ evolves as the process $\tilde{C}_t$ introduced in the proof of Proposition 4.7. Combining this fact with (10), (12) and Proposition 4.7, we have

$$ \left| e^{-\tilde{\lambda}_0 t} E_{(c,x)} [g(\bar{C}_t, \bar{\chi}_t) \mathbf{1}_{t<\bar{\tau} \wedge T_{na}}] - \eta_a(c) \nu_a[g(\cdot, \mathbf{a})] \right| \le C e^{-\gamma t} V_a(c), \quad \forall t \ge 0, \forall (c, x) \in E_a. $$

Next, we make use of Theorem 3.1 of Champagnat and Villemonais (2022) for the discrete time process $(\bar{C}_n, \bar{\chi}_n)_{n \in \mathbb{N}}$, in the situation corresponding to Assumption A1 therein. In order to do so, we set $D_1 = E_a$, $D_2 = E_{na}$, $j_{0,P} = 0$, $\theta_{0,P} = e^{-\tilde{\lambda}_0}$, $\gamma = e^{-\|b_a\|_{\infty}} < \theta_{0,P}$, $c_1 = 1$, $W_R \equiv \mathbf{1}_{(c,x) \in D_2}$ and $W_P = V_a \mathbf{1}_{(c,x) \in D_1}$. We deduce that for $g: E \to \mathbb{R}$ such that $|g| \le V := W_P + W_R$,

$$ \left| e^{-\tilde{\lambda}_{0n}} E_{(c,x)} [g(\bar{C}_n, \bar{\chi}_n) \mathbf{1}_{n<\bar{\tau}}] - \eta(c,x) \nu[g] \right| \le C \beta^n V(c,x), \quad \forall (c,x) \in E, \forall n \in \mathbb{N}, \quad (16) $$

for some $\beta \in (0, 1)$. Now fix $h \ge 0$ and consider

$$ g_h : (c, x) \mapsto E_{(c,x)} [g(\bar{C}_h, \bar{\chi}_h) \mathbf{1}_{h<\bar{\tau}}]. $$

Applying (16) to $g_h$ and using the Markov property at time $h$, we obtain

$$ \left| e^{-\tilde{\lambda}_{0n}} E_{(c,x)} [g(\bar{C}_{n+h}, \bar{\chi}_{n+h}) \mathbf{1}_{n+h<\bar{\tau}}] - \eta(c,x) \nu[g_h] \right| \le C \beta^n V(c,x) \xrightarrow{n \to +\infty} 0. $$
Taking $(c, x) \in E_a$, we thus observe that $\nu$ is a quasi-limiting distribution and hence a quasi-stationary distribution for the process (this is a classical result from the theory of quasi-stationary distributions, see for instance the works of Méléard and Villemonais (2012); Collet et al. (2013); van Doorn and Pollett (2013)). In particular, there exists $\lambda \le 0$ such that $P_\nu((\bar{C}_t, \bar{\chi}_t) \in \cdot, t < \bar{\tau}) = e^{\lambda t} \nu(\cdot)$, for all $t \ge 0$. The above convergence shows that $\lambda = \tilde{\lambda}_0$ and we deduce that $\nu[g_h] = e^{\tilde{\lambda}_0 h} \nu[g]$. Finally, we proved that

$$ \left| e^{-\tilde{\lambda}_0 n} \mathrm{E}_{(c,x)} [g(\bar{C}_{n+h}, \bar{\chi}_{n+h}) \mathbf{1}_{n+h<\bar{\tau}}] - \eta(c,x) e^{\tilde{\lambda}_0 h} \nu[g] \right| \le C \beta^n V(c,x). $$

Up to a change in the constant $C$ and setting $\gamma = \ln^{1/\beta}$, this concludes the proof of the theorem. $\square$

# 5 Numerical simulations

This section is dedicated to the numerical exploration of the theoretical results of the previous section. In Section 5.1, we consider the model without telomerase (i.e. with $q \equiv 0$) and study its limiting population size, through the so-called Hayflick limit, as a function of the parameters of the model. In Section 5.2, we consider the model with telomerase and study, first the values of the Malthusian parameter $\lambda_0$, defined by (8), as a function of the parameters of the model, and second the convergence of the process as stated by Theorem 4.5 when $\lambda_0 > 0$. Notice that, contrary to the Hayflick limit which takes transitory effects into account, $\lambda_0$ and the limit of the process only describe the asymptotic average behaviour, and do not depend on the initial condition.

Note that under Assumption 4.2, the expected number of particles grows exponentially and thus a naïve numerical simulation of the process is not appropriate. Instead, we use an interacting particle approximation scheme based on genetic algorithms, described in Algorithm 1 in appendix, which allows us to replace the mean semigroup by the average behaviour of a fixed size population. We refer the reader to the articles of Del Moral (2004); Del Moral (2013) for detailed results and methods allowing to deduce the average dynamics of the population (including in particular $t \mapsto \mathrm{E}_{(c_0,a)}(N_t)$ as well as the Malthusian parameter $\lambda_0$) from simulations of a population with fixed size $N$.

For simulations presented in Sections 5.1 and 5.2, we compute the expectation of 100 runs of Algorithm 1 with $N = 10000$ particles. Where possible, we have chosen parameter values that are consistent with empirical data or existing literature. In particular, the experiments described by Chow et al. (2012) suggest that $\omega$ should be chosen in the range [12b, 300b]. However, when information on possible parameter values is not available, we have chosen the values arbitrarily.

We remind the reader that we have included a table of notation in Appendix C containing a description of the model parameters and a reference to where they were first introduced in this article.
## 5.1 Without telomerase

In this section, we consider the model without telomerase. As such, this model describes the dynamics of a population of cells whose telomere lengths can only decrease with time. More precisely, we choose $q \equiv 0$ and by default

$$L_{min} = 2 \text{ kb}, \omega = 0.2 \text{ kb}, K = 46, c_0 = \prod_{i=1}^{K}(8\text{kb}, 8\text{kb}), \qquad (17)$$

so that all the telomeres in the population have initial length 8kb (the unit kb refers to 1 kilobase, sometimes written as kbp for kilobase pair), and

$$s_a(c) = e^{-(\min c - L_{min})/r_s} \text{ with } r_s = 100, b_a(c) = 1, d_{na}(c) = d_a(c) = 0. \qquad (18)$$

In particular, when all the cells have (at least) one telomere shorter than $L_{min}$, they are all in a non-active state, and the population can no longer evolve.

One quantity of interest is the final size of the population, which we denote by $N_\infty$. The Hayflick limit (Hayflick and Moorhead, 1961; Shay and Wright, 2000) of a cell population describes the total number of doubling of the population. Mathematically, we define this limit as

$$H = \log_2 \mathbb{E}_{(c_0, \mathbf{a})} [N_\infty], \qquad (19)$$

where $c_0$ describes the telomere lengths of an initial cell. In what follows, we will investigate (numerically) the dependence of the Hayflick limit on certain model parameters, namely the number of chromosomes $K$, the minimal telomere length for active cells $L_{min}$, the overhang length $\omega$, and the parameter $r_s$ of the “deactivation” rate $s_a$.

*Remark 5.1.* In biological experiments, the Hayflick limit is better described by $\log_2 \frac{N_\infty}{N_0}$. However, under mild assumptions, $\frac{N_\infty}{N_0}$ can be approximated by $\mathbb{E}(\frac{N_\infty}{N_0})$ when the initial size of the population $N_0$ is large (this is a consequence of the law of large numbers and the branching property in our model). In addition, if all the initial cells are in state $(c_0, \mathbf{a})$, $\mathbb{E}(\frac{N_\infty}{N_0}) = \mathbb{E}_{(c_0, \mathbf{a})} [N_\infty]$. This justifies our mathematical definition of the Hayflick limit $H$.

We first present the evolution of the expected population size for the above choice of parameters, see Figure 2. In the initial population, all telomere lengths are sufficiently large so that the non-active population is negligible. From (18), the expected population size grows exponentially with rate $b_a(c) = 1$ (see the second figure in Figure 2 where the total expected population size follows the red dotted line until $t \approx 19$). Then the growth slows down until the population stabilises at $\mathbb{E}_{(c_0, \mathbf{a})}(N_\infty)$ (see figure on the left), which corresponds to the time when the population of active cells vanishes (see figure on the right).

In Figure 3, we represent the Hayflick limit as a function of the number of chromosomes $K$; the minimal length, $L_{min}$, of telomeres in active cells; the overhang $\omega$, and the parameter $r_s$ associated with the “deactivation” rate $s_a$. We make the following observations.

* In Figure 3a, we observe that the Hayflick limit is decreasing with respect to the number of chromosomes $K$. This is not surprising since one expects that the larger the number
<figure><img src="image_6.png" /><figcaption>Figure 2: Expected population size as a function of time $t \ge 0$, for the model without telomerase ($q \equiv 0$) with parameters as in (17) and (18). Left: the expected total population dynamics $E_{(c_0,a)}(N_t)$, with the proportion of active cells $E_{(c_0,a)}(X_t(E_a))$ (red) and non active cells $E_{(c_0,a)}(X_t(E_{na}))$ (blue). Centre: zoom of the start of the dynamics in the left figure. Right: the dynamics of the expected number of active cells $E_{(c_0,a)}(X_t(E_a))$. In each of the three graphs, the dotted red line is the curve of the function $t \mapsto e^t$.</figcaption></figure>

of chromosomes, the smaller (in law) the minimal telomere length in the cell. Therefore increasing the number of chromosomes increases the probability of deactivation (with our choice of $s_a$ given by (18)), as well as the probability that a cell is ‘non-active’ after a division.

With our choice for $L_{min}$, $\omega$ and $c_0$ given in (17), the minimal theoretical value of the Hayflick limit (when $K \to \infty$) is $(8000 - 2000)/200 = 30$, however we observe that it remains around 45 for a large but realistic number of chromosomes. For the case with 46 chromosomes, which corresponds to the setting of Figure 2, we observe a Hayflick limit of 45.9, which is comparable to the experimentally measured interval [40, 60] (see Hayflick and Moorhead (1961) and Verma et al. (2020)).

* In Figure 3b, we observe that the Hayflick limit decreases linearly with respect to the minimal length $L_{min}$ of active cells. This is due to the fact that increasing the minimal telomere length of active cells increases the number of non-active cells and hence decreases the number of times the population doubles.

* Similarly, the Hayflick limit is also decreasing with respect to the overhang $\omega$ (see Figure 3c) as well as the parameter $r_s$ of the deactivation rate, $s_a$ (see Figure 3d).

## 5.2 With telomerase

We consider now the model with telomerase and we choose default parameters satisfying Assumptions 4.1, 4.3 and 4.4. Our values of $q(\ell)$ are derived from the empirical measures described by Teixeira et al. (2004) and the mechanism described by Greider (2016). We choose the following geometric distribution for the probability that telomerase lengthens a telomere, and the following uniform distribution for the distribution of the size of the increase:

$$q(\ell) = q_0 2^{-(\ell - L_{min})/L_{min}} \quad \text{and} \quad \mu_f = \mathcal{U}([0, M]), \qquad (20)$$
<figure><img src="image_7.png" /><figcaption>(a) Hayflick limit as a function of K.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Hayflick limit as a function of $L_{min}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Hayflick limit as a function of $\omega$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Hayflick limit as a function of $r_s$.</figcaption></figure>

Figure 3: Hayflick limit (19) for different values of K, $L_{min}$, $\omega$ and $r_s$, all other parameters being as in (17) and (18). Note that, the numerical approximation method used in Algorithm 1 implies some uncertainties in the computation of the Hayflick limit, which are slightly visible, mainly in Figures 3a and 3d.
<figure><img src="image_11.png" /><figcaption>Figure 4: Malthusian parameter $\lambda_0$ (8) for different values of $q_0$ (left), $M$ (center) and $r_s$ (right), for the model with the telomerase mechanism given by (20) with the default values $q_0 = 0.7$, $M = 600$ and all other parameters being as in (17) and (18).</figcaption></figure>

where $q_0 \in [0, 1]$ and $M > 0$. By default, we choose $q_0 = 0.7$, $M = 600$ and the other parameters remain as in (17) and (18).

If $\lambda_0 > 0$ then the expected population size goes to infinity at rate $\lambda_0$ and hence the Hayflick limit is equal to infinity. On the other hand, if $\lambda_0 < 0$, then the asymptotic expected population size of active cells decreases exponentially fast and eventually goes extinct. Finally, if $\lambda_0 = 0$, then the population of active cells is critical and does not asymptotically increase nor decrease exponentially fast.

In order to understand the effect of the biological parameters on the behaviour of the population (survival, growth rate, Hayflick limit), we compute $\lambda_0$ for different values of $s_a$, $q$ and $\mu_f$ (equivalently, $r_s$, $q_0$ and $M$, respectively). When $\lambda_0 < 0$, we compute the Hayflick limit and, when $\lambda_0 > 0$, we illustrate the convergence of the telomere length distribution in the population of cells to the limiting distribution $\nu$, as stated in Theorem 4.5.

In Figure 4, we plot $\lambda_0$ as a function of the parameters $q_0$, $M$ and $r_s$. We make the following observations.

* As expected, the eigenvalue $\lambda_0$ increases when the parameter $q_0$ increases, since this increases then the probability of telomerase lengthening a telomere.

* Similarly, when $M$ increases, so does $\lambda_0$ since this increases the upper bound on the amount a telomere can be lengthened by.

* In the left-hand figure, we see that $\lambda_0$ decreases monotonically with $r_s$, the parameter that appears in the deactivation rate $s_a$. The reasoning is the same as for the Hayflick limit. In fact, as we see in Figure 3d for the model without telomerase, increasing the parameters $r_s$ is not beneficial for the growth of the population (seen in terms of the number of times the population doubles). Similarly, in this case, it is not beneficial when the telomerase is active (seen here in terms of the population size growth rate $\lambda_0$).

As previously mentioned, we also consider parameter values which yield $\lambda_0 < 0$, i.e. a subcritical population. In particular, in Figure 4, we see that certain values of $q_0$ and $M$ yield a subcritical population, despite the presence of telomerase. In this case, the asymptotic
<figure><img src="image_12.png" /><figcaption>Figure 5: Hayflick limits (19) for different values of $q_0$ (left) and $M$ (right) leading to $\lambda_0 < 0$, for the model with the telomerase mechanism given by (20) with the default values $M = 600$ (left picture) and $q_0 = 0.7$ (right picture) and all other parameters being as in (17) and (18).</figcaption></figure>

average population size of active cells decays exponentially until extinction, so that the total population size eventually plateaus out, and thus the Hayflick limit is finite.

To study this in more detail, in Figure 5, we represent the Hayflick limit as a function of the parameters $q_0$ and $M$ for these values. As expected, the monotonic behaviour is the same as in Figure 4.

On the other hand, when $\lambda_0$ is positive, the Hayflick limit is infinite. This is the case, for example, for the default set of parameters given at the start of the subsection. Figure 6 represents the evolution of certain functions of the distribution of telomere lengths in cells for this default set of parameters. For convenience, we have plotted the distribution of the minimum, maximum and mean (in each cell) of the telomere lengths. We observe the convergence of the distribution toward the distribution $\nu$, as stated in Theorem 4.5.

Finally, in Figure 7, we have plotted these distributions at the final time of the simulation for different values of $K$. We see that, the higher the number of chromosomes $K$, the higher the variance of the telomere length distribution. In particular, this figure shows that as the number of chromosomes increases, the average of the distribution of the minimum telomere length decreases. Since in addition, the deactivation and removal of cells from the system are determined by the length of the shortest telomeres in the cells, this entails that $\lambda_0$ decreases with $K$, as demonstrated by Table 1.

<table><caption>Table 1: Malthusian parameter $\lambda_0$ defined in (8) for the different values of $K$ used in Figure 7.</caption><thead><tr><td>Chromosomes number K</td><td>Malthusian parameter λ<sub>0</sub></td></tr></thead><tbody><tr><td>1</td><td>0.9992918</td></tr><tr><td>2</td><td>0.9986076</td></tr><tr><td>20</td><td>0.9863393</td></tr><tr><td>46</td><td>0.9692310</td></tr><tr><td>90</td><td>0.9416240</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 6: Dynamics of certain characteristics of the (normalized) population density for the model with the telomerase mechanism given by (20) with the default values $q_0 = 0.7$ and $M = 600$ and all other parameters being as in (17) and (18) (in particular with $K = 46$). More precisely, we plot the normalization of the distributions $E_{(c_0,a)}(\sum_{i=1}^{N_t} \mathbf{1}_{g(c_i(t),x_i(t)) \in a})$, with $g(c_i(t), x_i(t)) = \mathbf{1}_{x_i(t)=a} \min c_i(t)$ (blue), $g(c_i(t), x_i(t)) = \mathbf{1}_{x_i(t)=a} \max c_i(t)$ (red) and $g(c_i(t), x_i(t)) = \mathbf{1}_{x_i(t)=a} \text{mean } c_i(t)$ (green), where mean $c_i(t)$ denotes the average telomere length in $c_i(t) \in (\mathbb{N} \times \mathbb{N})^K$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 7: Distribution of the minimal (blue), mean (green) and maximal (red) lengths of telomere for the limit distribution $\nu$ given by Theorem 4.5, for different values of the number of chromosomes $K$ in cells for the model with the telomerase mechanism given by (20) with the default values $q_0 = 0.7$ and $M = 600$ and all other parameters being as in (17) and (18).</figcaption></figure>
# 6 Discussion

We built a probabilistic individual based model for the telomere dynamics in a population of cells. We studied theoretically and numerically the dynamics of this model, with and without telomerase. In particular, we estimated numerically the effect of several parameters on the Hayflick limit and the growth of the size of the cell population.

Although the probabilistic model was built to mimic the biological mechanisms involved in telomere length dynamics, we made several assumptions in order to make the model more tractable for theoretical and numerical analysis, and in particular, to keep the number of parameters in the model sufficiently small. However, in some situations, one may wish to consider the following mechanisms: abrupt telomere shortening (Špoljarić et al., 2019), single strand breaks (von Zglinicki et al., 2000) and oxidative stress (von Zglinicki, 2002; Arkus, 2005). This would require one to modify our model to allow $\omega$ to be random and to include rare abrupt events. In our model, we assumed that chromatids are divided randomly between daughter cells during mitosis and there is no biological evidence against it. In some organism stem cells, biased DNA segregation during mitosis has been observed in vitro (Karpowicz et al., 2009; Ranjan and Chen, 2022). These observations only concern stem cells and still need to be confirmed in vivo, but may also be included in our model by removing the assumption that the random variables $B_j$ are not independent.

Additional mechanisms may also be studied mathematically and numerically. Namely, one may take into account the reactivation of telomerase in inactive cells (e.g. Telomerase reverse transcriptase in cancer cells (Dratwa et al., 2020)), fluctuation of the telomerase intensity through time or the population (depending on the environment or on the state of the population, leading to probabilistic models in random environments or with density dependence) or on/off commutator mechanisms, typically modelled by a Piecewise Deterministic Markov Process representing the density of a commutator enzyme (see Tucey and Lundblad (2014) for a biological description of commutators). To model logistic constraints, one may also wish to consider a spatial model (see e.g. David and Perthame (2021)) or include the dynamics of a shared limited resource in the model (as e.g. for chemostat models, see Fritsch et al. (2015)). Finally, one may study the influence of the dependence of the overhang $\omega$ on the length of the telomere, empirically observed by Huffman et al. (2000); Rahman et al. (2008).

We finally mention the possibility of modelling an age-dependent process. At the expense of additional technicalities and an increase in model complexity, one may choose to study an age-dependent model, where the different rates described in our model depend on the age of the cell, see e.g. the work of Olayé and Tomasevic (2024). From a probabilistic perspective, this lies in the setting of Crump-Mode-Jagers branching processes (Jagers, 1975).

In our numerical simulations, we used parameter values inspired by empirical measures from the biological literature. Another natural approach, which will be the subject of further research, is to estimate the parameters from empirical data. Since most experimental data concerning human cells are available for differentiated cells, such as leukocytes (see e.g. Toupane et al. (2019) where the distribution of telomere length in humans is studied), in future work, we will study population processes to model multi-tissue organisms where stem cells differentiate into specialized cells.
# Statements and Declarations

## Acknowledgements

The authors would like to thank Anne Gégout-Petit for her support and useful discussions throughout the project. This work was supported partly by the french PIA project “Lorraine Université d'Excellence”, ANR-15-IDEX-04-LUE and received government funding managed by the Agence Nationale de la Recherche under the France 2030 program, reference ANR-23-EXMA-0005. Simulations are run on the babycluster of the Institut Élie Cartan de Lorraine.

## Competing interest

The authors have no relevant financial or non-financial interests to disclose.

## Compliance with Ethical Standards

The authors report there are no conflict of interest to declare.

## Data availability statement

No datasets were generated or analysed during the current study.

# A Dynkin's lemma

In this section, we discuss a useful result that is used in several places in the article. The full version of this result and further details can be found in (Horton and Kyprianou, 2023, Theorem 2.1). We note that this result is sometimes referred to as Dynkin's lemma due to the fact that a slightly less general version of this result was first introduced by Dynkin (2002).

Let us first introduce some necessary notation. Let $E$ be a Polish space (it is possible to make $E$ even more general but it will not be necessary for this article), which will denote the state space. We will also append an extra point $\dagger \notin E$ to $E$, that will denote the “cemetery” state, i.e. where particles are sent when they die. We let $\xi = (\xi_t, t \ge 0)$ denote a Markov process on $E$ with lifetime $\zeta := \inf\{t \ge 0 : \xi_t = \dagger\}$. Note that we consider $\dagger$ to be an absorbing state meaning that $\xi_t \in \{\dagger\}$ for all $t \ge \zeta$. We also let $\mathbf{P}_x$ denote the law of $\xi$ conditional on $\xi_0 = x$ with corresponding expectation operator $\mathbf{E}_x$.

In what follows, we will let $B(E)$ denote the space of bounded, measurable functions on $E$ such that $f(\dagger) = 0$ for all $f \in B(E)$. We will also let $B^+(E)$ denote the subset of $B(E)$ consisting only of non-negative functions. Then, for $f \in B^+(E)$, $t \ge 0$ and $x \in E$, we define

$$ \mathbb{Q}_t[f](x) := \mathbf{E}_x[f(\xi_t)] $$

to be the expectation semigroup associated to $(\xi, \mathbf{P}_x)$. We now extend the notion of expectation semigroups. Suppose that $\gamma \in B(E)$ and define, for $x \in E$, $t \ge 0$ and $f \in B(E)$,

$$ \mathbb{Q}_t^{\gamma}[f](x) := \mathbf{E}_x \left[ e^{\int_0^t \gamma(\xi_s) ds} f(\xi_t) \right]. $$
The weight $e^{\int_0^t \gamma(\xi_s)ds}$ is often referred to as a multiplicative potential. This potential can be thought of as weighting the Markov process, consequently giving path a notion of importance. For example, if $\sup_{x \in E} \gamma(x) \le 0$, then the multiplicative potential acts as a penalisation term, since this weight will be less than or equal to one.

Often, the expectation semigroup $(\mathbb{Q}_t^\gamma, t \ge 0)$ forms the basis of an evolution equation:

$$ \chi_t(x) := \mathbb{Q}_t^\gamma[g](x) + \int_0^t \mathbb{Q}_s^\gamma[h_{t-s}](x)ds, \quad t \ge 0, x \in E, g \in B^+(E), \quad (21) $$

and where $h : [0, \infty) \times E \to [0, \infty)$ satisfies $\sup_{s \le t} |h_s| \in B^+(E)$.

We are now in a position to state a result that will be useful throughout the article, which offers an alternative form of (21).

**Theorem A.1.** Suppose that $|\gamma| \in B^+(E)$ and that $\sup_{s \le t} |h_s| \in B^+(E)$ for all $t \ge 0$. If $(\chi_t, t \ge 0)$ is represented by (21), then it also solves

$$ \chi_t(x) = \mathbb{Q}_t[g](x) + \int_0^t \mathbb{Q}_s[h_{t-s} + \gamma\chi_{t-s}](x)ds, \quad t \ge 0, x \in E. \quad (22) $$

Note that under an extra condition, the converse of this statement is true however, we will not state the details since we only use the result given above in this article.

The full proof of the above result is essentially integration by parts and we refer the interested reader to Horton and Kyprianou (2023). A quick calculation that can help the reader to convince themselves of this result is to compute the generator associated with (21) to obtain a differential version of (22).

The above result can be thought of as a form of calculus by which the multiplicative potential in $\mathbb{Q}^\gamma$ is removed and appears instead as an additive potential in the integral term. Heuristically, it shows us how to transfer mass within the evolution equation, resulting in two different stochastic interpretations of the solution to (21) or equivalently, (22). For simplicity, let us suppose that $h \equiv 0$. Then, on the one hand, (22) suggests a solution given by the expectation of a branching process whose particles move according to the semigroup $\mathbb{Q}$ and mass is created according to $\gamma$. On the other hand, (21) offers another interpretation of the solution as the expectation of a single particle, weighted by the potential $\gamma$, that is, the solution is given by the semigroup $\mathbb{Q}^\gamma$. The inclusion of $h_{t-s}$ in the integral term simply allows for a more general setting, which we now illustrate with an example. Further examples can be found in (Horton and Kyprianou, 2023, Chapter 2).

**Example A.2.** Consider the case where $\xi$ is a Brownian motion in $E = \mathbb{R}$ such that, at rate $\alpha \in B^+(\mathbb{R})$, the particle jumps to a new position in $\mathbb{R}$ according to the law $\mu$. Hence, if $T$ is the time of the first such jump, then

$$ \mathrm{Pr}(T > t|\sigma(B_s, s \le t)) = e^{-\int_0^t \alpha(B_s)ds}, \quad t \ge 0. \quad (23) $$

Let $(X_t, t \ge 0)$ denote this process with law $P_x$ conditional on $X_0 = x \in \mathbb{R}$ and set $\chi_t(x) = E_x[f(X_t)]$ for $f \in B^+(\mathbb{R})$. By conditioning on whether $T \le t$ or $T > t$, it can be shown (see p.23 of Horton and Kyprianou (2023) for the details) that $\chi_t(x)$ satisfies (21) with

$$ \mathbb{Q}_t^{-\alpha}[f](x) = \mathbf{E}_x \left[ e^{-\int_0^t \alpha(B_s)ds} f(B_t) \right] $$
and

$$h_{t-s} = \alpha(\cdot) \int_{\mathbb{R}} \chi_{t-s}(y) \mu(dy).$$

Theorem A.1 then shows that, equally, $\chi_t(x)$ satisfies

$$\chi_t(x) = \mathbf{E}_x[f(B_t)] + \int_0^t \mathbf{E}_x \left[ \alpha(B_s) \left( \int_{\mathbb{R}} \chi_{t-s}(y) \mu(dy) - \chi_{t-s}(x) \right) \right] ds.$$

This second equation can be interpreted as follows. At time $t$, either no jump has occurred and so the particle is at position $B_t$, or with rate $\alpha$, the ‘parent’ particle is replaced by a ‘new’ particle whose position is given by $\mu$.

## B Algorithm

Here we present the interacting particle approximation scheme used to produce the simulations in Section 5. The following algorithm allows one to simulate a population with fixed size, $N$ say. Initially, each particle evolves according to an independent copy of a (sub)Markov process. When a particle is killed, it is resampled from the surviving population and the particles then continue to evolve independently. In this case, the process evolves according to $(\bar{C}_t, \bar{X}_t)_{t \ge 0}$ defined in the proof of Theorem 4.5. We refer the reader to the article of Del Moral (2004); Del Moral (2013) for further details of algorithms of this type.

*Remark B.1.* The estimation of $\lambda_0$ from Algorithm 1 by the methods of Del Moral (2004); Del Moral (2013) is valid under the convergence stated in Theorem 4.5. Despite the fact that Theorem 4.5 was proved under restrictive assumptions (including $\lambda_0 > 0$), numerical simulations suggest that the convergence also holds true for all the parameters choices of Section 5. Based on these considerations, we allow ourselves to use the same numerical methods to estimate the parameter $\lambda_0$ across Section 5. In particular, we consider certain parameter regimes where $\lambda_0 \le 0$, in order to illustrate that our result should hold under less restrictive assumptions, and to provide the reader with a more complete picture. We leave the question of the possible generalization of Theorem 4.5 to these parameters open.
sample $(c_i(0), x_i(0))_{i \in \llbracket 1, N \rrbracket} = (c_0, \mathbf{a})_{i \in \llbracket 1, N \rrbracket} \{\text{initialization}\}

t \leftarrow 0

while $t \le T$ do

$\tau \leftarrow \sum_{i=1}^{N} 1_{x_i(t)=\mathbf{a}}(2b_{\mathbf{a}}(c_i(t)) + s_{\mathbf{a}}(c_i(t)) + d_{\mathbf{a}}(c_i(t)) + 1_{x_i(t)=\mathbf{na}}d_{\mathbf{na}}(c_i(t)))$
$\Delta t \sim \text{Exp}(\tau)$

$i \sim \sum_{j=1}^{N} \frac{1_{x_i(t)=\mathbf{a}}(2b_{\mathbf{a}}(c_i(t)) + s_{\mathbf{a}}(c_i(t)) + d_{\mathbf{a}}(c_i(t)) + 1_{x_i(t)=\mathbf{na}}d_{\mathbf{na}}(c_i(t))}{\tau}$
$(c_j(t + \Delta t), x_j(t + \Delta t))_{j \in \llbracket 1, N \rrbracket, j \neq i} = (c_j(t), x_j(t))_{j \in \llbracket 1, N \rrbracket, j \neq i}$

$u \sim \mathcal{U}[0, 1]$

if $x_i(t) = \mathbf{a}$ and $u \le \frac{2b_{\mathbf{a}}(c_i(t))}{2b_{\mathbf{a}}(c_i(t)) + s_{\mathbf{a}}(c_i(t)) + d_{\mathbf{a}}(c_i(t))}$ then
{cell division}

$(m_j, n_j)_{j \in \llbracket 1, K \rrbracket} = c_i(t)$

$B_j \sim \text{Ber}(1/2), j = 1, \dots, K$

$L_j \sim \mu_f, L'_j \sim \mu_f, j = 1, \dots, K$

$\chi_j \sim \text{Ber}(q(m_j - \omega B_j)), j = 1, \dots, K$

$\chi'_j \sim \text{Ber}(q(n_j - \omega(1 - B_j))), j = 1, \dots, K$

$m_j \leftarrow m_j - \omega B_j + L_j \chi_j, j = 1, \dots, K$

$n_j \leftarrow n_j - \omega(1 - B_j) + L'_j \chi'_j, j = 1, \dots, K$

if $\min_{1 \le j \le K}\{m_j, n_j\} \le L_{\min}$ then
{non-viable cell}

$(c_i(t + \Delta t), x_i(t + \Delta t)) \leftarrow ((m_j, n_j)_{j \in \llbracket 1, K \rrbracket}, \mathbf{na})$

else

{viable cell}

$(c_i(t + \Delta t), x_i(t + \Delta t)) \leftarrow ((m_j, n_j)_{j \in \llbracket 1, K \rrbracket}, \mathbf{a})$

end if

else

if $x_i(t) = \mathbf{a}$ and $u \le \frac{2b_{\mathbf{a}}(c_i(t)) + s_{\mathbf{a}}(c_i(t))}{2b_{\mathbf{a}}(c_i(t)) + s_{\mathbf{a}}(c_i(t)) + d_{\mathbf{a}}(c_i(t))}$ then
{cell deactivation}

$(c_i(t + \Delta t), x_i(t + \Delta t)) \leftarrow (c_i(t), \mathbf{na})$

else

{cell removal: the cell is uniformly sampled on existing cells}

$(c_i(t + \Delta t), x_i(t + \Delta t)) \sim \mathcal{U}((c_j(t + \Delta t), x_j(t + \Delta t))_{j \in \llbracket 1, K \rrbracket, j \neq i})$

end if

end if

$t \leftarrow t + \Delta t$

end while

Algorithm 1: Fleming-Viot type scheme: simulation of $N$ processes $(C_t, \mathcal{X}_t)_{t \ge 0}$ defined in Section 3.2. Each process is simulated by a Gillespie Algorithm (Gillespie, 1977). When one process goes extinct, it is replaced by one alive cell, drawn uniformly among the alive processes.

C Notation

Here we provide the reader with a table of notation for the parameters that are used in the numerical simulations in Section 5. We give the notation used, a description and the page number where it was first introduced.
<table><thead><tr><th>Notation</th><th>Description</th><th>Reference</th></tr></thead><tbody><tr><td>(m, n)</td><td>Chromosome indexed by its telomeres</td><td>p.5</td></tr><tr><td>ω</td><td>Overhang</td><td>p.5</td></tr><tr><td>K</td><td>Number of chromosomes per cell</td><td>p.7</td></tr><tr><td>L<sub>min</sub></td><td>Minimal telomere length of active cells</td><td>p.7</td></tr><tr><td>(c, a) = ((m<sub>j</sub>, n<sub>j</sub>)<sup>K</sup><sub>j=1, a</sub>)</td><td>Active cell</td><td>p.7</td></tr><tr><td>(c, na) = ((m<sub>j</sub>, n<sub>j</sub>)<sup>K</sup><sub>j=1, na</sub>)</td><td>Non active cell</td><td>p.7</td></tr><tr><td>c<sub>0</sub></td><td>Telomere lengths of each initial cell</td><td>p.20</td></tr><tr><td>q(l)</td><td>Probability that telomerase acts</td><td>p.5</td></tr><tr><td>q<sub>0</sub></td><td>Parameter of q(l)</td><td>p.21</td></tr><tr><td>μ<sub>f</sub></td><td>Telomere lengthen law</td><td>p.5</td></tr><tr><td>M</td><td>Maximal telomere lengthen</td><td>p.21</td></tr><tr><td>b<sub>a</sub></td><td>Division rate</td><td>p.8</td></tr><tr><td>s<sub>a</sub></td><td>Deactivation rate</td><td>p.8</td></tr><tr><td>r<sub>s</sub></td><td>Parameter of s<sub>a</sub></td><td>p.20</td></tr><tr><td>d<sub>x</sub></td><td>Removal rate</td><td>p.8</td></tr><tr><td>N<sub>t</sub></td><td>Cell number in the population</td><td>p.7</td></tr><tr><td>N<sub>∞</sub></td><td>Final population size</td><td>p.20</td></tr><tr><td>H</td><td>Hayflick limit</td><td>p.20</td></tr><tr><td>λ<sub>0</sub></td><td>Malthusian parameter</td><td>p.12</td></tr></tbody></table>

References

Arino, O., Kimmel, M. and Webb, G. F. (1995). Mathematical modeling of the loss of telomere sequences. *Journal of Theoretical Biology*, 177(1):45–57.

Arkus, N. (2005). A mathematical model of cellular apoptosis and senescence through the dynamics of telomere loss. *Journal of Theoretical Biology*, 235(1):13–32.

Aubert, G. (2014). Telomere dynamics and aging. *Progress in Molecular Biology and Translational Science* 125, 89–111.

Bansaye, V., Cloez, B. and Gabriel, P. (2019). Ergodic behavior of non-conservative semigroups via generalized Doeblin's conditions. *Acta Applicandae Mathematicae*, 1–44.

Blackburn, E. H. and Collins, K. (2011). Telomerase: an RNP enzyme synthesizes DNA. *Cold Spring Harbor Perspectives in Biology*, 3(5):a003558.

Blasco, M. A. (2005). Telomeres and human disease: ageing, cancer and beyond. *Nature Reviews Genetics*, 6(8):611–622.

Bourgeron, T., Xu, Z., Doumic, M. and Teixeira, M. T. (2015). The asymmetry of telomere replication contributes to replicative senescence heterogeneity. *Scientific Reports*, 5(1):1–11.
Champagnat, N. and Villemonais, D. (2023). General criteria for the study of quasi-stationarity. *Electronic Journal of Probability*, 28:1–84.

Champagnat, N. and Villemonais, D. (2018). Uniform convergence of penalized time-inhomogeneous markov processes. *ESAIM Probability and Statistics*, 22:129–162.

Champagnat, N. and Villemonais, D. (2022). Quasi-stationary distributions in reducible state spaces. *ArXiv preprint*, arXiv:2201.10151.

Chow, T. T., Zhao, Y., Mak, S. S., Shay, J. W. and Wright, W. E. (2012). Early and late steps in telomere overhang processing in normal human cells: the position of the final RNA primer drives telomere shortening. *Genes & Development*, 26:1167–1178.

Collet, P., Martínez, S. and San Martín, J. (2013). *Quasi-stationary distributions*. Probability and its Applications (New York). Springer, Heidelberg. Markov chains, diffusions and dynamical systems.

Cox, A. M., Harris, S. C., Horton, E. L. and Kyprianou, A. E. (2019). Multi-species neutron transport equation. *Journal of Statistical Physics*, 176(2):425–455.

Dao Duc, K. and Holcman, D. (2013). Computing the length of the shortest telomere in the nucleus. *Physical Review Letters*, 111:228104.

David, N. and Perthame, B. (2021). Free boundary limit of a tumor growth model with nutrient. *Journal de Mathématiques Pures et Appliquées*, 155:62–82.

Del Moral, P. (2004). *Feynman-Kac formulae*. Probability and its Applications (New York). Springer-Verlag, New York. *Genealogical and interacting particle systems with applications*.

Del Moral, P. (2013). *Mean field simulation for Monte Carlo integration*, volume 126 of *Monographs on Statistics and Applied Probability*. CRC Press, Boca Raton, FL.

Del Moral, P. and Miclo, L. (2002). On the stability of nonlinear Feynman-Kac semigroups. *Annales de la Faculté des Sciences de Toulouse : Mathématiques, Série 6*, 11(2):135–175.

D'Mello, M. J., Ross, S. A., Briel, M., Anand, S. S., Gerstein, H. and Paré, G. (2015). Association between shortened leukocyte telomere length and cardiometabolic outcomes: systematic review and meta-analysis. *Circulation: Cardiovascular Genetics*, 8(1):82–90.

Dratwa, M., Wysoczańska, B., Łacina, P., Kubik, T. and Bogunia-Kubik, K. (2020). TERT–regulation and roles in cancer formation. *Frontiers in Immunology*, 11:589929.

Dynkin, E.B. (2002). *Diffusions, Superdiffusions and Partial Differential Equations*. *American Mathematical Society Colloquium Publications (American Mathematical Society, Providence)*, vol. 50.
Entringer, S., de Punder, K., Buss, C. and Wadhwa, P. D. (2018). The fetal programming of telomere biology hypothesis: an update. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 373(1741):20170151.

Eugène, S., Bourgeron, T. and Xu, Z. (2017). Effects of initial telomere length distribution on senescence onset and heterogeneity. *Journal of Theoretical Biology*, 413:58-65.

Fritsch, C., Harmand, J. and Campillo, F. (2015). A modeling approach of the chemostat. *Ecological Modelling*, 299:1-13.

Gillespie, D. (1977). Exact stochastic simulation of coupled chemical reactions. *Journal of Physical Chemistry*, 81(25):2340-2361.

Grach, A. (2013). Telomere shortening mechanisms. *The Mechanisms of DNA Replication*, 10:55244.

Greider, C. W. (2016). Regulating telomere length from the inside out: the replication fork model. *Genes & Development*, 30(13):1483-1491.

Hayat, M. (2017). Chapter 1 - overview of autophagy. In Hayat, M., editor, *Autophagy: Cancer, Other Pathologies, Inflammation, Immunity, Infection and Aging*, pages 3-90. Academic Press.

Haycock, P. C., Burgess, S., Nounu, A., Zheng, J., Okoli, G. N., Bowden, J., Wade, K. H., Timpson, N. J., Evans, D. M., Willeit, P., et al. (2017). Association between telomere length and risk of cancer and non-neoplastic diseases: a mendelian randomization study. *JAMA Oncology*, 3(5):636-651.

Hayflick, L. and Moorhead, P. S. (1961). The serial cultivation of human diploid cell strains. *Experimental cell research*, 25(3):585-621.

Hemann, M. T., Strong, M. A., Hao, L.-Y. and Greider, C. W. (2001). The shortest telomere, not average telomere length, is critical for cell viability and chromosome stability. *Cell*, 107(1):67-77.

Hiyama, E. and Hiyama, K. (2007). Telomere and telomerase in stem cells. *British Journal of Cancer*, 96(7):1020-1024.

Horton, E. and Kyprianou, A. E. (2023). *Stochastic Neutron Transport: And Non-local Branching Markov Processes*. Springer Nature.

Huffman, K. E., Levene, S. D., Tesmer, V. M., Shay, J. W. and Wright, W. E. (2000). Telomere shortening is proportional to the size of the G-rich telomeric 3'-overhang. *Journal of Biological Chemistry*, 275(26):19719-19722.

Jafri, M. A., Ansari, S. A., Alqahtani, M. H. and Shay, J. W. (2016). Roles of telomeres and telomerase in cancer, and advances in telomerase-targeted therapies. *Genome Medicine*, 8:1-18.
Jagers, P. (1975). *Branching processes with biological applications*. John Wiley & Sons, London.

Karpowicz, P., Pellikka, M., Chea, E., Godt, D., Tepass, U., and van der Kooy, D. (2009). The germline stem cells of Drosophila melanogaster partition DNA non-randomly. *European journal of cell biology*, 88(7):397-408.

Lee, K. H. and Kimmel, M. (2020). Stationary distribution of telomere lengths in cells with telomere length maintenance and its parametric inference. *Bulletin of Mathematical Biology*, 82(12):1-16.

Levy, M. Z., Allsopp, R. C., Futcher, A. B., Greider, C. W. and Harley, C. B. (1992). Telomere end-replication problem and cell aging. *Journal of Molecular Biology*, 225(4):951-960.

Louzon, M., Coeurdassier, M., Gimbert, F., Pauget, B. and de Vaufléury, A. (2019). Telomere dynamic in humans and animals: Review and perspectives in environmental toxicology. *Environment International*, 131.

Méléard, S. and Villemonais, D. (2012). Quasi-stationary distributions and population processes. *Probability Surveys*, 9:340-410.

Muraki, K., Nyhan, K., Han L., and Murnane, J.P. (2012). Mechanisms of telomere loss and their consequences for chromosome instability. *Frontiers in Oncology*, 2(135).

Olayé, J. and Tomasevic, M. (2024). Long-time behaviour of a multidimensional age-dependent branching process with a singular jump kernel. *ArXiv preprint*, arXiv:2201.10151.

Olofsson, P. and Kimmel, M. (1995). Stochastic models of telomere shortening. *Mathematical Biosciences*, 158(1):75-92.

Olovnikov, A. M. (1996). Telomeres, telomerase, and aging: origin of the theory. *Experimental Gerontology*, 31(4):443-448.

Portugal, R., Land, M. and Svaiter, B. (2008). A computational model for telomere-dependent cell-replicative aging. *Biosystems*, 91(1):262-267.

Qi, Q., Wattis, J. A. and Byrne, H. M. (2014). Stochastic simulations of normal aging and Werner's syndrome. *Bulletin of Mathematical Biology*, 76(6):1241-1269.

Rahman, R., Forsyth, N. R. and Cui, W. (2008). Telomeric 3'-overhang length is associated with the size of telomeres. *Experimental Gerontology*, 43(4):258-265.

Ranjan, R. and Chen, X. (2022). Mitotic drive in asymmetric epigenetic inheritance. *Biochemical Society transactions*, 50(2):675-688.

Rhyu, M. S. (1995). Telomeres, telomerase, and immortality. *Journal of the National Cancer Institute*, 87(12):884-894.
Shammas, M. A. (2011). Telomeres, lifestyle, cancer, and aging. *Current Opinion in Clinical Nutrition and Metabolic Care*, 14(1):28.

Shay, J. W. and Wright, W. E. (2000). Hayflick, his limit, and cellular ageing. *Nature Reviews Molecular Cell Biology*, 1(1):72–76.

Špoljarić, A. M., Rubelj, I. and Huzak, M. (2019). Mathematical model and computer simulations of telomere loss. *Journal of Theoretical Biology*, 465:78–89.

Starkweather, A. R., Alhaeeri, A. A., Montpetit, A., Brumelle, J., Filler, K., Montpetit, M., Mohanraj, L., Lyon, D. E. and Jackson-Cook, C. K. (2014). An integrative review of factors associated with telomere length and implications for biobehavioral research. *Nursing Research*, 63(1):36.

Teixeira, M. T., Arneric, M., Sperisen, P. and Lingner, J. (2004). Telomere length homeostasis is achieved via a switch between telomerase-extendible and-nonextendible states. *Cell*, 117(3):323–335.

Toupance, S. and Benetos, A. (2019). *Telomere Dynamics and Aging Related Diseases*, 66–90. Cambridge Handbooks in Psychology. Cambridge University Press.

Toupance, S., Villemonais, D., Germain, D., Gegout-Petit, A., Albuisson, E. and Benetos, A. (2019). The individual's signature of telomere length distribution. *Scientific Reports*, 9(1):685.

Trentesaux, C. and Riou, J.-F. (2010). *Sénescence et immortalité cellulaire*. *Bulletin du Cancer*, 97(11):1275-1283.

Tucey, T. M. and Lundblad, V. (2014). Regulated assembly and disassembly of the yeast telomerase quaternary complex. *Genes & Development*, 28(19):2077–2089.

van Doorn, E. A. and Pollett, P. K. (2013). Quasi-stationary distributions for discrete-state models. *European Journal Operational Research*, 230(1):1–14.

Velleret, A. (2023). Exponential quasi-ergodicity for processes with discontinuous trajectories. *ESAIM: Probability and Statistics*, 27:867–912.

Verma, A., Verma, M. and Singh, A. (2020). Animal tissue culture principles and applications. In *Animal Biotechnology*, pages 269–293. Elsevier.

von Zglinicki, T. (2002). Oxidative stress shortens telomeres. *Trends in Biochemical Sciences*, 27(7):339–344.

von Zglinicki, T., Pilger, R. and Sitte, N. (2000). Accumulation of single-strand breaks is the major cause of telomere shortening in human fibroblasts. *Free Radical Biology and Medicine*, 28(1):64–74.

Wattis, J. A., Qi, Q. and Byrne, H. M. (2020). Mathematical modelling of telomere length dynamics. *Journal of Mathematical Biology*, 80(4):1039–1076.
Xu, Z., Duc, K. D., Holcman, D. and Teixeira, M. T. (2013). The length of the shortest telomere as the major determinant of the onset of replicative senescence. *Genetics*, 194(4):847-857.

Xu, Z. and Teixeira, M. T. (2019). The many types of heterogeneity in replicative senescence. *Yeast*, 36(11):637-648.