Modelling stem cell ageing: a multi-compartment continuum approach

Stem cells are important to generate all specialized tissues at an early life stage, and in some systems, they also have repair functions to replenish the adult tissues. Repeated cell divisions lead to the accumulation of molecular damage in stem cells, which are commonly recognized as drivers of ageing. In this paper, a novel model is proposed to integrate stem cell proliferation and differentiation with damage accumulation in the stem cell ageing process. A system of two structured PDEs is used to model the population densities of stem cells (including all multiple progenitors) and terminally differentiated (TD) cells. In this system, cell cycle progression and damage accumulation are modelled by continuous dynamics, and damage segregation between daughter cells is considered at each division. Analysis and numerical simulations are conducted to study the steady-state populations and stem cell damage distributions under different damage segregation strategies. Our simulations suggest that equal distribution of the damaging substance between stem cells in a symmetric renewal and less damage retention in stem cells in the asymmetric division are favourable strategies, which reduce the death rate of the stem cells and increase the TD cell populations. Moreover, asymmetric damage segregation in stem cells leads to less concentrated damage distribution in the stem cell population, which may be more robust to the stochastic changes in the damage. The feedback regulation from stem cells can reduce oscillations and population overshoot in the process, and improve the fitness of stem cells by increasing the percentage of cells with less damage in the stem cell population.

Stem cells are important to generate all specialized tissues at an early life stage, and in some systems, they also have repair functions to replenish the adult tissues. Repeated cell divisions lead to the accumulation of molecular damage in stem cells, which are commonly recognized as drivers of ageing. In this paper, a novel model is proposed to integrate stem cell proliferation and differentiation with damage accumulation in the stem cell ageing process. A system of two structured PDEs is used to model the population densities of stem cells (including all multiple progenitors) and terminally differentiated (TD) cells.
In this system, cell cycle progression and damage accumulation are modelled by continuous dynamics, and damage segregation between daughter cells is considered at each division. Analysis and numerical simulations are conducted to study the steadystate populations and stem cell damage distributions under different damage segregation strategies. Our simulations suggest that equal distribution of the damaging substance between stem cells in a symmetric renewal and less damage retention in stem cells in the asymmetric division are favourable strategies, which reduce the death rate of the stem cells and increase the TD cell populations. Moreover, asymmetric damage segregation in stem cells leads to less concentrated damage distribution in the stem cell population, which may be more robust to the stochastic changes in the damage. The feedback regulation from stem cells can reduce oscillations and population overshoot in the process, and improve the fitness of stem cells by increasing the percentage of cells with less damage in the stem cell population.
has the potential to remain as a stem cell (self-renewal) or becomes a cell with a more specialized function (differentiation) (figure 1a). If the division produces the same type of cells, such as two stem cells and two differentiated cells, it is called a symmetric division. On the other hand, if the division produces two different types of cells, i.e. a stem cell and a differentiated cell, then it is called an asymmetric division. Available data suggest that most stem cells are able to switch between symmetric and asymmetric divisions, and the balance between these modes is controlled by various internal and external signals to produce appropriate numbers of stem and differentiated cells [1][2][3].
Stem cells can be found in most mammalian tissues, and they participate in tissue repair in response to damage and maintaining tissue homeostasis [4,5]. Research suggests that the decline in adult tissue maintenance and the increase in cancer formation might be a consequence of stem cell ageing [6]. Although age-related manifestation in stem cell population and function differ across tissues and organisms, decline in regenerative capacity due to depletion or dysfunction of stem cells is a hallmark [7]. Protein aggregates, dysfunction organelles and DNA damage are commonly identified as factors of ageing [4,5]. To slow down the accumulation of ageing factors, a hypothesis suggests that mitotic cells (actively dividing cells) might asymmetrically segregate damage away from the cell whose fate is to become a new stem cell [8]. The asymmetric inheritance of cellular components in dividing cells was first observed in yeast, and has been extensively studied over the past decades. In yeast, carbonylated proteins, extrachromosomal ribosomal DNA circles and dysfunction mitochondria are retained by a mother cell during asymmetric division, while its daughter is rejuvenated with little damage [9][10][11]. Recent evidence suggests that stem cells may employ a similar mechanism to protect one progeny from ageing [12][13][14][15][16][17][18] (figure 1b). An asymmetric partition of damaged proteins in stem cell division was observed in adult flies' intestine and germline [12], mammalian stem cells [16], and murine neural stem cells [18]. Despite these emerging findings, it still remains elusive how stem cells cope with damage accumulation and how this is related to stem cell proliferation and differentiation.
Driven by the lack of knowledge of mechanisms regulating stem cell maintenance and differentiation in the ageing process, mathematical models have been employed to address key questions and provide quantitative insights into stem cell renewal and differentiation, as well as the decline of cellular functions in the ageing process. Several mathematical models were proposed to study stem cell population and ageing, which fall into two categories: individual-based modelling and continuous population modelling [19,20].
Individual-based modelling simulates individuals or agents that have a unique set of state variables and usually interact with each other in the local environment [20,21]. The advantages of this approach include that it can take stochastic effects into consideration to describe phenomena on the level of individual cells, and that detailed molecular dynamics and cell-cell interaction can be incorporated. In [22], a stochastic model of stem cell organization was introduced to explain the observed heterogeneity of haematopoietic stem cells by the stochastic switching between the growth environments and the selforganizing process based on within-tissue plasticity properties. Assuming that cell proliferation is negatively affected by telomere shortening, an agent-based stochastic model [23] was proposed to study telomere-dependent stem cell replicative ageing. This model provides a good approximation of the qualitative growth of cultured human mesenchymal stem cells. In [24], mutation accumulation in large populations of stem cells was modelled by a discrete-time branching process where each division produces 0, 1 or 2 stem cell daughters, each of which randomly accumulates a mutation. This model demonstrated that symmetric division could reduce the risk of accumulating phenotypically silent heritable damage in individual stem cells. However, a major drawback of individual-based modelling is the computational inefficiency, especially when the population is large. When a population is large and homogeneous, continuous population model using ODEs [25][26][27][28][29][30][31] or PDEs [32][33][34][35][36][37][38] are more appropriate to describe the population dynamics. In [39], a maturity-structured two-compartment model was proposed to study the regulation of mammalian red blood cell production. The model consists of two transport equations describing population densities of mitotic cells and post-mitotic cells, and an ODE describing hormone dynamics. The model showed that a perturbation of blood-donation type leads to damped oscillatory return to normal status, and that an elevated random peripheral destruction of red blood cells leads to sustained oscillations. In [40], a three-compartment ODE model was applied to study the dynamics of stem cells, transit-amplifying cells and terminally differentiated (TD) cells in the olfactory epithelium of mice. The authors identified conditions on parameters for the stability of the system when negative feedback loops are present either as Hill functions or in more general forms. Their analysis suggested that two factors, autoregulation of the proliferation of transit-amplifying progenitor cells and low death rate of TD cells, enhance the stability of the system. In [41], the authors applied the mean-field approach to approximate an agent-based model for studying heterogeneity within the haematopoietic stem cell population. Their proposed PDE model can capture the key structure of the model including the 'age'-structure of stem cells and improve the efficiency of the numerical algorithms. In [42], a system of PDEs was used to model mutation accumulation hierarchy and differentiation hierarchy of cells with stem cells on the top level and to examine cancer development and growth. In their model, maturity is treated as a continuous variable, while the number of mutation accumulation and telomere shortening are treated as discrete cell classes. The boundary conditions describe transition among different cell classes at division: cells lose telomeres and acquire mutations. The study showed that the more mutation classes and higher proliferation rate are sufficient to explain the faster growth of the cancer cell population.
To understand the effects of different damage segregation strategies on stem cell ageing, we propose a novel model to integrate stem cell proliferation and differentiation with damage accumulation in the stem cell ageing process, and feedback regulation from stem and TD cells. A system of hyperbolic PDEs is constructed to model two compartments in cell lineage: mitotic cells (stem cells) and post-mitotic cells (TD cells). It is assumed that the cell cycle progression of stem cells is a continuous process while stem cell division is discrete. The boundary conditions of the PDEs model the stem cell renewal and differentiation at division when damage segregation takes place. Cell death is modelled as an outcome of damage accumulation. Stem cell proliferation and differentiation are regulated by feedbacks from the population of TD cells and stem cells. Ageing effect is modelled through the inhibition from the damage accumulation on stem cell proliferation and self-renewal. Our analysis and numerical simulations are carried out to compare the effects of different regulation and damage segregation strategies on population dynamics and stem cell fitness, which have not been discussed in the previous studies of stem cell regulation.
Our simulations suggest that equal distribution of the damaging substance between stem cells in symmetric renewal and less damage retention in stem cell in asymmetric division are favourable strategies, which reduce the death rate of stem cells and increase TD cell populations. Also, asymmetric damage segregation in stem cells leads to less concentrated damage distribution in stem cells population, which may be more robust to the stochastic change in damage. Compared to the feedbacks solely from TD cells, the feedback regulation from stem cells (autoregulation) can reduce oscillations and population overshoot in the process, and improve the fitness of stem cells by increasing the percentage of stem cells with less damage in the stem cell population.
This paper is structured as follows. The general description of our model is given in §2. In §3, a simple model without feedback regulations is presented to analyse the relation between population dynamics and various parameters. In §4, two more complex models with feedbacks from TD cells and stem cells are proposed to study different regulation mechanisms and the effect of segregation strategies.

Model description
In our mathematical model, a simplified conceptual model, we consider two types of cells: mitotic cells, which include stem cells and multiple progenitor cells, and post-mitotic cells, which include all TD cells. For simplicity, we call mitotic cells stem cells, and denote its population density by S; and we call post-mitotic cells TD cells, and denote its population density by T. To consider the evolution of cell populations, we define both S(t, p, a) and T(t, p, a) as functions of time t and two other continuous biological state variables: cell cycle progression p and damage level a. To describe cell maturation and ageing processes, we make the following assumptions: -When the amount of damage a accumulates and reaches a certain threshold a Ã , stem cells die by an apoptosis-like process as a result of ageing via damage accumulation; similar to stem cells, when the amount of damaged substance reaches a certain threshold a c , TD cells are removed by an apoptosis-like process. -Cell cycle progression p is an indicator variable, for which a stem cell divides when p increases to a threshold p Ã , unless damage has already reached a threshold a Ã , in which case the stem cell is removed before division. Although TD cells no longer divide, for simplicity, we keep cell cycle progression variable p in T, but assume the upper boundary for p is infinity.
By conservation law, a system of transport equations are derived to describe the evolution of S and T: and where V p and U p are rates of cell cycle progression for stem and TD cells, respectively; V a and U a are rates of damage accumulation for stem and TD cells, respectively. Note that various of feedbacks may regulate cell cycle progression and cell damage accumulations, and these functions will be specified in the following sections.
The boundary conditions at p = 0 describe the reproduction/division process accompanied by the segregation of damage substances. First, we assume that there are three types of cell divisions: (i) The daughter cells after division are two stem cells; this occurs with probability δ 1 .
(ii) The daughter cells after division are two TD cells; this occurs with probability δ 2 . (iii) The daughter cells after division are one stem cell and one TD cell; this occurs with probability δ 3 .
One of these three types of cell division occurs at the end of cell cycle p = p Ã with probability δ 1 , δ 2 or δ 3 , where δ 1 + δ 2 + δ 3 = 1 (figure 2a). These three probabilities may be regulated by various feedbacks and will be specified in later sections.  Upon the completion of a cell cycle, the stem cell cycle progression p will be reset to zero and the damaged proteins are inherited from mother to daughters. The damage inheritance can be described by transition kernels r i,S (a, a 0 ) and r i,T (a, a 0 ). Based on the above assumptions, the boundary conditions for S and T cells at p = 0 are as follows: and The first term of the right-hand side of (2.3) represents the stem cell production process through type (i) cell division; the second term of the right-hand side of (2.3) represents the stem cell process through type (iii) cell division; the first term of the right-hand side of (2.4) represents the TD cell production process through type (ii) cell division; the second term of the right-hand side of (2.4) represents the TD cell production process through type (iii) cell division. The transition kernel r i,S (a, a 0 ) describes how daughter stem cells with damage a come from the mother cells with damage a 0 after the i-th type of division; the transition kernel r i,T (a, a 0 ), describes how daughter TD cells with damage a come from the mother cells with damage a 0 in the i-th type of division. The transition kernels satisfy the following conservation conditions: Since type (i) and type (ii) divisions induce two daughter stem cells and two daughter TD cells, respectively, the integrals of r 1,S (a, a 0 ) and r 2,T (a, a 0 ) equal 2; since type (iii) division induces one daughter stem cell and one daughter TD cell, the integrals of r 3,S (a, a 0 ) and r 3,T (a, a 0 ) both equal one. The last three integrals are based on the assumption that the amount of damage is conserved during cell division. No-flux conditions on the boundary a = 0 are imposed to ensure conservation of population in the direction of a and For TD cells, we assume that the population density of TD is zero when p is very large The populations of stem and TD cells at time t are given by the integrals S(t, p, a) dp da and P T (t) ¼

Model without feedback regulations
In this section, we start with the simple model where no feedbacks are included, and we additionally make the following assumptions to simplify the model (figure 2b): -Cell cycle progression velocities V p = v p and U p = u p are constants, so are the damage accumulation rates V a = v a and U a = u a . -The portions of three types of divisions δ i are constants.
Under these assumptions, the transition kernels in the boundary conditions are Dirac delta functions and the boundary conditions in (2.3) and (2.4) become and Now, we study the role of the damage segregation rules in the long-term behaviour of the stem cell population. Since the segregation rule is fixed, after sufficiently long time, the cellular damage at the end of cell cycle, temporarily assuming no death, converges to a limit damage band where ω 1 and ω 2 are the minimum and the maximum of α 1 , α 2 , γ 1 , respectively. The derivation of damage limit band can be found in appendix A 1.
The population dynamics turns out to depend on the proportions of three division types of stem cells and the position of the lethal threshold a Ã /v a with respect to the limit damage band [ p Ã /v p (1 − ω 1 ), p Ã /v p (1 − ω 2 )]. Before proceeding, we define a key lumped parameter for studying the long-term population behaviour Here, the parameter f r is called the self-renewal fraction. Using the limit damage band and the self-renewal fraction, we can find some conditions for different long-term behaviours.
Proposition 3.1. Stem cells become extinct, if the renewal fraction f r < 1/2 or the limit damage band is completely above the lethal threshold, i.e. a Ã /v a < p Ã /v p (1 − ω 1 ).
Proposition 3.2. Assume that the limit damage band lies completely below the lethal threshold, i.e. a Ã /v a > p Ã /v p (1 − ω 2 ) and on the following assumptions on the renewal fraction and the lethal threshold a Ã /v a lies within the limit damage band : Under the above assumptions in (3.3), there are many situations of damage segregation rules, and more importantly not all of the situations are biologically meaningful. Here, we will focus on one situation: the damage retention in stem cells through asymmetric division is smaller than the damage segregation portions in symmetric stem cell renewal, i.e., This setting is based on the biological observation that asymmetric division is a favourable mechanism for stem cell lineage to remove the damage. For avoiding the case that the lethal threshold a Ã /v a lies below p Ã /v a (1 − α 1 ), i.e. all asymmetrically renewing stem cells are destined to die, we further assume that By the assumptions (3.3)-(3.5), we obtain that : Under the assumption (3.6), the stem cell population may approach zero or blow up in the long-term behaviour. The following proposition provides a condition to guarantee extinction of stem cells, as well as a condition for population blow-up, with the proof given in appendix A 2.
If there exists a combination of δ i such that (2f r ) n À d n 1 , 1, then the stem cell population is eventually zero.
then the stem cell population goes to infinity, or is bounded below to avoid stem cell extinction.
From proposition 3.3, we can obtain several necessary conditions to maintain the stem cell population. Since the term (2f r ) n À d n 1 is increasing with n and goes to infinity when n tends to infinity, n has to be small enough to satisfy the condition (2f r ) n À d n 1 , 1 in proposition 3.3. From the definition of n, we observe that if the difference between a Ã /v a and p Ã /v p (1 − α 2 ) becomes small, or the difference between 1/(1 − α 2 ) and 1/(1 − γ 1 ) becomes large, or the right-hand side p Ã =v p (1 À a 2 ) À ( p Ã =v p )(1=(1 À a 2 ) À (1=(1 À g 1 )))a m 2 becomes large for each m, the value of n will increase and then (2f r ) n À d n 1 , 1 cannot be satisfied in most of the combinations of δ 1 , δ 2 and δ 3 . This result provides some conditions of the parameters to prevent stem cell extinction.
Let us discuss the above three possibilities. To reduce the difference between a Ã /v a and p Ã /v p (1 − α 2 ), the ratio of cell cycle progression speed v p to damage accumulation speed v a should be large enough and close to p Ã /a Ã (1 − α 2 ); to increase the difference between 1/(1 − α 2 ) and 1/(1 − γ 1 ), damage retention in asymmetric division should decrease, i.e. γ 1 should be small; to increase the right-hand side p Ã =v p (1 À a 2 ) À ( p Ã =v p )(1=(1 À a 2 ) À (1=(1 À g 1 )))a m 2 for each m, damage distribution in self-renewal should become more symmetric, i.e. α 1 should increase to close to 0.5. In conclusion, when v p /v a increases, γ 1 decreases or α 1 increases, it may provide a better condition to maintain the stem cell population. These results are supported by the numerical simulations shown in figure 3.
royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191848 In figure 3, we study the long-term population dynamics with different combinations of the division probabilities δ 1 and δ 3 . We uniformly generate 600 pairs of δ 1 and δ 3 with 0 < δ 1 + δ 3 ≤ 1 and 2δ 1 + δ 3 ≥ 1. We skip the region 2δ 1 + δ 3 < 1 as it implies stem cell extinction whatever other parameters are (by proposition 3.1). For each pair of δ 1 and δ 3 , the system is solved by the numerical method described in appendix A 5. The blue regions in figure 3 represent the cases of stem cell extinction and the yellow regions represents the non-extinction cases (the stem cell population may tend to a constant non-zero value, approach infinity, or keep oscillating). When v p /v a increases (in figure 3a, v p /v a = 1.67; in figure 3b, v p /v a = 2), γ 1 decreases (in figure 3c, γ 1 = 0.3; in figure 3d, γ 1 = 0.05) or α 1 increases (in figure 3e, α 1 = 0.1; in figure 3f, α 1 = 0.45), the blue regions become smaller (the yellow regions become larger), which implies that there are more combinations of parameters (δ 1 , δ 3 ) allowing stem cell survival.
For most of the combinations of parameters, the stem cell populations in the models without feedback regulation blow up to infinity or diminish to zero. Although the no-regulation assumption is not realistic, the simple model not only provides us a guidance on the selection of parameters used in the model with feedbacks but also reveals that population dynamics are results of all factors: the cell cycle progression of stem cell, damage accumulation, fractions of divisions and damage segregation rules. In the next section, we will consider the combinations of parameters that guarantee exponential growth in the stem cell population and study feedback regulations that could lead to non-zero steady-state population.

Model with feedback regulations
Biological evidence shows that some mammalian stem cells can switch between symmetric and asymmetric divisions in response to external and internal regulations [1,3]. For example, both epidermal [43] and neural [44] progenitors change from the primarily symmetric division that expands the stem-cell pool during embryonic development to primarily asymmetric in mid to late gestation. It is also observed that nervous [45] and haematopoietic [46] stem cells in adults can divide symmetrically to replace lost cells through injury, although they divide asymmetrically under  steady-state conditions. In particular, two types of feedbacks have been proposed [34,47,48] (also see figure 4a): long-range feedbacks responding to the population of TD cells, and short-range feedbacks acting in an autocrine fashion from stem cells. These two types of feedbacks regulate stem cell population through controlling two types of parameters: (i) the speed of stem cell division, V p (or stem cell cycle progression), (ii) the probabilities of three types of cell division, δ i (differentiation versus renewal).
However, the underlying mechanisms are not fully understood. In addition, the effect of damage segregation rules on cell population is completely unknown. In the following, we use mathematical models with different feedbacks and damage segregation rules to study which types of feedback and damage segregation mechanisms are more plausible given the experimental observations.
For simplicity, we assume that damage segregation rules are fixed, i.e. α i , β i , γ i are constants. Based on the results from the simple model in §3, to ensure that the stem cell population will blow up to infinity without any feedback regulation, additional assumptions are made: -Self-renewal fraction f r . 1 2 . -Damage accumulation speed V a = v a is slower than cell cycle progression of stem cells.

Feedback only from TD cells
The long-range feedback, which acts through signals sent by differentiated cells and inhibits stem cell division and self-renewal, has been biologically observed in numerous tissues, including muscle [49], bone [50], skin  [51], nervous system [52] and haematopoietic systems [53]. Despite this wealth of data, there is less understanding of the exact mechanisms of feedback regulation. A significant number of mathematical models have been developed to explore the possible mechanisms behind feedback regulations [26,27,32,[34][35][36][37][38][39][40]54,55]. The dynamics of signalling molecule s(t) can be described by a simple ODE as follows: where α is a constant synthesis rate, and the degradation is proportional to the level of s and affected by the cell population P of stem or TD cells. Since the dynamics of the signalling molecules take place on a faster time scale than the processes of cell proliferation and differentiation, we assume that the feedback signal can be approximated by a quasi-steady-state solution [25,27,32,[34][35][36][37][38][39][40]54,55]. By properly rescaling s, the quasisteady state of s is given by a Hill function where k is a regulation constant to account for the sensitivity to the cell population and m is the Hill exponent. The function in (4.2) reflects the assumption that the signal intensity achieves its maximum in the absence of cells and decreases asymptotically to zero if the number of cells increases. Note that Hill functions are widely used to describe ligand-receptor interactions, which also makes them natural choices to model the actions of secreted feedback factors [37]. According to biological evidence, we may model the feedback from TD cells in the following ways (figure 4b): when stem cell or TD cell population is small (in the early stage of development or with drastic loss due to injury), symmetric division predominates over asymmetric division; during the stable stage (mid and late gestation or tissue homeostasis), stem cells switch to asymmetric division. Such feedbacks may be added to the model by modifying division fractions δ i as follows: where i ∈ {1, 2}, d 0 i are basal division fractions, k T i are regulation constants, m T is the Hill exponent and P T (t) is the TD cell population. The remaining division fraction δ 3 is defined as 1 − δ 1 − δ 2 .
Recent studies on cell cycles of neural stem cells have shown evidence that TD cells may be a source of signalling molecules that inhibits cell cycle progression of stem cells [56]. Therefore, besides the above regulation on cell fate decisions, negative feedback from TD cells can also be involved in stem cell proliferation V p (cell cycle progression speed), i.e. excessive number of TD cell may slow down the proliferation of stem cell, and hence reduce the population of both stem cells and TD cells. Thus, V p can be modified as where v p is the initial cell cycle progression speed and k T v is a regulation constant. For other velocities, we set V a = v a , U a = u a and U p = u p .
Among all parameters, we are most interested in the effect of segregation rules α i , β i , γ i on the population dynamics. A set of reasonable estimations of the parameters other than α i , β i , γ i are chosen based on appropriate biological and mathematical assumptions, and are shown in table 1. The detailed reasoning of the selection of parameters is provided in appendix A.3.1. In brief, we assume that initially, the stem cells are expanding, the cell cycle progression is fast, and the symmetric renewal predominates in three types of division. In the rest of this subsection, we will focus on the discussion of the effect of different combinations of segregation rules α i , β i , γ i .
In the simulations, we consider the following situations: α 1 varies from 0.1 to 0.5, indicating that the symmetry in damage segregation increases between the two stem daughter cells; β 1 varies from 0.1 to 0.5, indicating that the symmetry in damage segregation increases between TD cells; γ 1 varies from 0.1 to 0.9, indicating that the damage retention in stem cells increases in asymmetric division. Under different segregation rules, the following aspects are studied: the population size of TD cells and the population ratio of TD cells to stem cells at steady state; the death rate of stem cells and the fraction of three types of division at steady state; and the dynamics of population evolution and the damage distribution in stem cell population at steady state.  does not affect population size at steady state. When the parameter β 1 , which determines how damage is distributed between two TD cells in the symmetric differentiation of stem cells, varies from 0.1 to 0.5, neither stem cell population nor TD cell population has a big change. This result is biologically meaningful, since the replenishment of TD cells is efficient compared to the loss of TD cells due to damage accumulation. In the following part, we will focus on the case with β 1 = 0.5. Figure 5 also shows that as damage segregation α i between stem cells becomes more symmetric, TD cell populations increase, while the ratios P T /P S are almost constant with a slight decrease. For fixed α 1 , as damage retention γ 1 in asymmetric division increases, both TD cell populations and the population ratios decrease. Interestingly, when α 1 = 0.5, although the populations of TD cells are very close for the cases γ 1 = 0.1, 0.3 and 0.5, the population ratios differ dramatically. Compared to γ 1 = 0.1, one needs a

Death rate of stem cells and fraction of three types of division
Other than studying the population size, we also compare the death rate of stem cells and the fraction of three types of division at steady state. In our model, the death rates of stem cell (r S ) and TD cell (r T ) are measured in the following way: 0 v a S(t, p, a Ã ) dp P S and r T ¼ Ð 1 0 u a T(t, p, a c ) dp P T , ( 4 :5) where the numerators are population out-fluxes due to death and the denominators are the total populations. Figure 6 shows that more equal distribution of damage in symmetric renewal and less damage retention in asymmetric division result in a lower death rate of stem cells and less symmetric division at steady state. According to our simulations, the segregation rules do not have much impact on the death rate of TD cells, which ranges from 2.130 × 10 −4 to 3.345 × 10 −4 . However, the death rate of stem cells changes dramatically as we vary the segregation rules. From figure 6a, we can see that the death rate of stem cells increases as α 1 decreases or γ 1 increases. The lowest death rate of stem cells is attained when damage is equally distributed between stem cells in symmetric renewal and damage retention is minimal in asymmetric division. To maintain the steady stem cell population and to replenish the short-lived TD cells, a higher death rate should be associated with a fast turnover of stem cells. Indeed, when we examine the fractions of three types of division, we find that the higher death rate of stem cells is always associated with the higher fraction of symmetric renewal (figure 6b). These results suggest that to maintain stabilized populations, tissues may have different mechanisms that involve different damage segregation rules and division rules.

Dynamics of population evolution and damage distribution of stem cell population
The comparison of population dynamics among different combinations of α 1 and γ 1 (figure 7a) shows that higher damage retention results in oscillations in population evolution. Oscillations start to appear as we increase the damage retention γ 1 in stem cells during asymmetric division (γ 1 = 0.7 in figure 7a). Also, the oscillations become more severe when the damage distribution among stem cells in symmetric renewal becomes more symmetric (α 1 = 0.5 and γ 1 = 0.7 in figure 7a). Moreover, population overshoots before steady state are observed in all cases with α 1 ≠ 0.5. The results above suggest that more equal distribution of damage in symmetric renewal and less damage retention in asymmetric division are favourable segregation rules in four aspects: populations converge to steady states faster without oscillations or severe population overshoot; the death rate of stem cells is much smaller than that of TD cells; asymmetric division predominates in three types of division at steady state; not only the size of TD cell population is larger, the population ratio P T /P S is also larger.
However, the asymmetric damage segregation between stem cells in self-renewal may have benefits that cannot be shown from the analysis of total populations as integrated results. Our PDE modelling approach allows us to investigate more details of population density and damage distribution of stem cells. Figure 7b gives the stem cell population densities corresponding to different segregation rules. As we studied in §3, damage segregation rules determine the limit damage band after a sufficiently long time. In particular, the less symmetric the segregation rules are, the wider the limit damage bands will be. Although asymmetric segregation rules may result in a greater percentage of stem cells inheriting more damage and accelerated death, it also leads to a higher percentage of healthier stem cells with less damage. This can be observed in the samples shown in figure 7b. When the symmetry of damage segregation in stem cell population increases, the damage distributions at the end of cell cycle for stem cells become more concentrated, as such symmetry increases. The concentrated damage distribution could be a disadvantage, since it is less resistant to external perturbations. In this sense, we think that some stem cells may sacrifice the lower death rate for more diversified damage distribution, in order to protect the stem cell pool from a possible unfavourable situation that may lead to a sudden increase in damage.

Feedback from both TD and stem cells
The maintenance of the stem cell pool is not only affected by signalling from mature cells, but also by the stem cell pool itself. Stem cells reside in the so-called stem cell niche, where both cellular and non-cellular components interact in order to control stem cell proliferation and differentiation [57][58][59][60]. The regulations from the stem cell population include two aspects: a negative feedback control of stem cell proliferation as a result of contact inhibition; and a self-renewal inhibition factor secreted by the stem cell niche, in which case the more stem cells there are, the less likely stem cells will divide symmetrically. Experiments also show that cells inheriting the majority of damage protein aggregates during asymmetric division have an increased cell cycle length and tendency to differentiate [13,14,18]. In addition to feedback regulations from differentiated cells, we are also interested in how the dynamics of populations would change if regulations from stem cells are included (figure 4c). Many modelling works (e.g. [34,40]) only consider the feedbacks from stem cell populations, due to the simplicity of their models. However, we include damage as a state variable in our model and could simulate feedbacks from both the stem cell population and the stem cell cellular damage level.
Thus, in our model, we assume that the excessive stem cell population and elevated stem cell cellular damage slow down stem cell proliferation and modify V p in (4.4) as is a decreasing function of a with sigmoid shape, in which k a 1 , a 0 1 , a 1 and b 1 are constants. We also assume that the excessive stem cell population inhibits stem cell symmetric division, and that the elevated stem cell cellular damage promotes stem cell differentiation via decreasing stem cell selfrenewal. As a result, we modify δ 1 and δ 2 in (4.3) in the same way as above ( 4 :7) and is another function of a sigmoid shape with constants k a 2 , a 0 2 , a 2 and b 2 . As a continuation of §4.1, our analysis still focuses on the discussion of segregation rules. In addition to the parameters in table 1, the values of more parameters are given in table 2, for which the detailed reasoning can be found in appendix A.3.2.
To examine the change in population dynamics after introducing feedbacks from stem cells, the same combinations of segregation rules as in §4.1 are considered: α 1 varies from 0.1 to 0.5, β 1 varies from 0.1 to 0.5, and γ 1 varies from 0.1 to 0.9. Here, we study the effect of the feedbacks from stem cells on the following two aspects: the population size of TD cells and the population ratio of TD cells to stem cells at steady state; and the dynamics of population evolution and the damage distribution of stem cell population at steady state. for fixed α 1 , as damage retention γ 1 in stem cell increases, both the TD cell population and the population ratio P S /P T decrease; for fixed γ 1 , as damage segregation between stem cells becomes more symmetric, the TD cell population increases, but the population ratio is almost unchanged. It is worth mentioning that the population ratio increases a little bit for each damage segregation configuration, after adding feedbacks from stem cells. This means that the regulation is more efficient, since one needs a smaller stem cell pool to generate the same number of TD cells.

Dynamics of population evolution and damage distribution of stem cell population
A notable consequence of adding feedbacks from stem cells is the reduction of the oscillations in population dynamics. Comparing figure 7a with figure 9a, after considering the feedbacks from stem cells, the oscillations in population dynamics disappear even when γ 1 is large, and that the population overshoot problem resolves when α 1 is small. Since the feedbacks from stem cells described by (4.7) and (4.6) also depend on stem cell cellular damage level, we want to examine the effect of such regulations on damage distribution at the end of cell cycle by comparing the models with/without the feedbacks from stem cells (figures 7b and 9b). Similar to §4.1, the damage distributions at the end of cell cycle become more concentrated, as the degree of symmetry in damage segregation increases. However, compared to the results in §4.1 (figure 7b), the model with the feedback from stem cells (figure 9b) provides better fitness, since under the same segregation rules, there are more cells with less damage at the end of cell cycle at steady state. This effect becomes most significant when the segregation rule is symmetric. This result suggests that slowing down the cell cycle progression of stem cells with a high level of damage and promoting such stem cells to differentiate can indeed improve the overall health of the stem cell population.

Conclusion
In this research, a novel model was developed to integrate stem cell proliferation and differentiation with damage accumulation in stem cell ageing process. A system of two structured PDEs are used to model stem cells (including all multiple progenitors) and TD cells. In our model, cell cycle progression and damage accumulation are continuous while division is discrete, and the damage segregation takes place at each division. Ageing effect is included through the inhibition from damage accumulation on stem cell proliferation and self-renewal. Regulations from TD cell and stem cell populations are incorporated through negative feedbacks on stem cell proliferation and symmetric division.
Our results showed that more equal distribution of damage between stem cells in symmetric renewal and less damage retention in stem cell in asymmetric division are still favourable segregation rules resulting in higher population size and greater population ratio; asymmetric damage segregation in stem cells leads to less concentrated damage distribution in stem cells population, which may be more robust to sudden increase in damage. These two results provide an answer for why some types of stem cells can switch between symmetric and asymmetric divisions in response to external and internal regulations [1,3].
Other than the feedbacks solely from TD cells, adding feedback regulations from stem cells can reduce oscillations and population overshoot in some situations with unfavourable damage segregation rules, say, in the situation where α 1 is small and γ 1 is large. Moreover, the stem cell feedback regulation can slow down the proliferation of stem cells with high level of damage to increase their tendency to differentiate. Overall, the feedback regulation system can improve the fitness of stem cells by increasing the percentage of stem cells with low level of damage. This result provides a new aspect to understand the role of the self-regulations in stem cell lineage [32,47,52,61]. In this study, our model provides a framework to study how the damage segregation rules at division affect the population dynamics. Moreover, our model is an excellent tool to understand the data which may result from observing the dynamics of regeneration with damage distribution in stem cell population after a knock-out process. In the future work, we can extend our study to the system with more than two stages, including progenitor cells [40,61] to get a better understanding for the system of stem cell lineage.
Data accessibility. This article has no additional data. The code for the numerical method can be found in the electronic supplementary material.
Authors' contributions. All authors were involved in the design of the model and the analysis of the numerical data; Y.W. and W.-C.L. carried out the numerical tests, and drafted the manuscript; C.-S.C. critically revised the manuscript. All authors gave final approval for publication and agree to be held accountable for the work performed therein.
Competing interests. We declare we have no competing interests.

Appendix A
A.1. Derivation of the limit damage band Suppose a stem cell currently has damage level a 0 in the beginning of a cell cycle. We want to record the damage in the following cycles. At the end of the first cycle, the damage accumulates to a 0 + p Ã /v p . After division, a stem daughter cell inherits damage (a 0 + p Ã /v p )ℓ 1 with ℓ 1 ∈ {α 1 , α 2 , γ 1 }, whose damage gets to (a 0 + p Ã /v p )ℓ 1 + p Ã /v p = a 0 ℓ 1 + ( p Ã /v p )(ℓ 1 + 1) at the end of the second cycle. After division, a stem daughter cell inherits damage (a 0 ℓ 1 + ( p Ã /v p )(ℓ 1 + 1))ℓ 2 with ℓ 2 ∈ {α 1 , α 2 , γ 1 }, and the damage accumulates to (a 0 ℓ 1 + ( p Ã /v p )(ℓ 1 + 1))ℓ 2 + tp Ã /v p = a 0 ℓ 1 ℓ 2 + ( p Ã /v p )(ℓ 1 ℓ 2 + ℓ 2 + 1) at the end of the third cycle. By induction, we can easily show that the damage in a stem daughter cell (if not dead yet) at the end of the m-th cycle is Let ω 1 and ω 2 be the minimum and the maximum of α 1 , α 2 , γ 1 . In the above, considering a fixed segregation portion ℓ i = ω j for all i ≥ 1, we get the damage level The limits as m → ∞ yield the limit damage band A.2. Proofs of the propositions A.2.1. Analytic solution Before we prove the propositions, we first obtain the analytic solution which can provide us a tool used in the proofs. To simplify the analysis, we make the following change of variables which will be royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191848 used in the proofs: and the domains of new variables are 0 e p t p ¼ p Ã =v p and 0 e a t a ¼ a Ã =v a . Given constant cell cycle progression and damage accumulation for stem cells, the biological meaning of t p and t a can be interpreted as t p is the duration of one cell cycle.
t a is the duration needed to accumulate damage to the lethal threshold from zero damage.
With the new variables, (2.1) and (2.2) become and where e S(t, e p, e a) ¼ S(t, p, a) and e T(t, e p, e a) ¼ T(t, p, a). For simplicity, we will drop all the tildes on the notations. The dynamics of stem cells is independent of TD cells, while the dynamics of TD cells is determined by stem cell through boundary condition (3.2). Once the behaviour of stem cell is determined, the behaviour of TD cell can be easily derived. So we mainly focus on the dynamics of stem cells.
After a change of variable, (A 1) together with the initial condition and boundary condition (3.1) can be solved by the method of characteristics. Due to the complexity of boundary conditions, no closed form of S(t, p, a) can be obtained, and the solution is expressed in a recursive form.
For the first cell cycle, i.e. 0 < t < t p , we have three cases (i) If t ≤ p and t ≤ a, the solution is determined by the initial condition S(t, p, a) ¼ S(0, p À t, a À t); (A 3) (ii) If t > a and a < p, the solution is determined by the boundary condition on a = 0 S(t, p, a) ¼ S(t À a, p À a, 0); (A 4) (iii) If t > p and p < a, the solution is determined by boundary condition on p = 0 S(t, p, a) ¼ S(t À p, 0, a À p): According to these three cases, S(t, p, a) can be solved as S(t, p, a) ¼ S(0, p À t, a À t) i f t p and t a, 0 i f t . a and a , p, þH(t, p, a, g 1 , t a ) d 3 g 1 S 0, t p À (t À p), a À p g 1 À (t À p) if t . p and p , a, where the function H is given by For the time beyond the first cell cycle, i.e. t > t p , by tracing back the process of damage accumulation by royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191848 one cycle, S(t, p, a) is solved as þ H(t p , p, a, g 1 , t a ) d 3 g 1 S t À t p , p, a À p g 1 À (t p À p) : Applying (A 8) iteratively and eventually (A 6), we can obtain the density of the stem cell population at any time.

A.2.2. Proof of proposition (3.1)
First, we suppose f r < 1/2. Consider any fixed time t ≫ t p . By (A 8), we have Then by integration and substitutions By the continuity of P S (t), it attains a maximum P max in a time interval [t 0 , t 0 + t p ] for a fixed t 0 . Repeated application of the above inequality, in view of the arbitrariness of t, implies that for any n ≥ 1, and for any t in [t 0 + n t p , t 0 + (n + 1)t p ], we have P S (t) (2f r ) n P max , which clearly establishes that lim t→∞ P S (t) = 0. For the second part, when the limit damage band exceeds the lethal threshold, i.e. t a < t p /(1 − ω 1 ), we let ɛ = t p /(1 − ω 1 ) − t a > 0. If we assume that (a − p)/ω < t a , where ω ∈ {α 1 , α 2 , γ 1 }, we have a , vt a þ p, where C = (1 − ω 2 )ɛ > 0 which is independent of a and p. Suppose that t ≪ t p , by (A 8), we have þ H(t p , p, a, g 1 , t a ) d 3 g 1 S t À t p , p, g(a, g 1 ) , where g(a, ω) = (a − p)/ω − (t p − p). For each term in the right-hand side, if (a − p)/ω ≥ t a , H is zero; if (a − p)/ω < t a , we have g(a, ω) < a − C. When g(a, ω) is less than or equal to zero, the term S(t − t p , p, g(a, ω)) = 0; when g(a, ω) is larger than zero, we can repeat the process and express the solution at t in terms of the solutions at t − nt p where n is a positive integer. Because g(a, ω) < a − C, the damage level will be reduced by a constant value in each process. Hence, with finite value n, all the terms S(t − nt p , p, g) will become zero as g is less than zero.
royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191848 Intuitively, the stem cell population always diminishes to zero, since after sufficiently long time the offspring of stem cells with initial zero damage, even in the slowest way of accumulating damage, i.e. ℓ i = ω 1 for every i, would gain damage to reach t a and die, and other offspring in the evolution gain more damage and would die sooner. This result is true, no matter what the proportions of three division types are.

A.2.3. Proof of proposition (3.2)
By the assumption that t a > t p /(1 − ω 2 ), no cell will die due to damage accumulation and f r completely determines the evolution of stem cell population. If f r > 1/2, the portion of generated stem cells through divisions is greater than 1, which means that the stem cell population is increasing to infinity. Indeed, for any t, let ω ∈ {α 1 , α 2 , γ 1 }, we define da dp: By the definition of H and change of variables with ω ≤ ω 1 and t a > t p /(1 − ω 2 ), t À t p , p, a da dp: By (A 8) and (A 9), t À t p , p, a da dp: Repeating this procedure, if t ¼ nt p þt andt [ [0, t p ), then we have p, a da dp: With the assumptions f r > 1/2 and ð t p 0 ð taÀt p þp 0 S 0, p, a da dp . 0: Now we proved that the population will blow up when n → ∞. Intuitively, after sufficiently long time, there is no death by damage accumulation since t a ≥ t p /(1 − ω 2 ). Hence in the situation that f r = 1/2, we can use a similar method to show that one of the offspring after stem cell division remains a stem cell and the stem cell population is conserved.

A.2.4. Proof of proposition (3.3)
We assume that (3.6) holds royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191848 20 We first prove (a). For integers m ≥ 1, define 1 À a 2 t p : Let P S (kt p ) denote the stem cell population at the end of the k-th cycle, when the damage levels in all cells are larger than a k ¼ (1 À g k 1 )=(1 À g 1 )t p , from the discussion in A 1. The damage level in stem daughter cells (if not dead yet) after n(k) more cycles would be greater than where n = n(k). If we set ℓ i = α 2 for all 0 ≤ i ≤ n(k), the definition of n(k) can imply that the damage level in (A 10) is larger than t a at the end of the (k + n(k))-th cycle, and the daughter cells in this situation will die. Thus, for the possible situations for the cell to survive at the end of the (k + n(k))-th cycle, there should be at least one ℓ i = γ 1 or α 1 for 1 ≤ i ≤ n(k). Clearly, the proportion of such stem daughter cells in the population is (2f r ) n(k) À d n(k) 1 : Suppose there exists a combination of m and δ i such that (2f r ) n(k) À d n(k) 1 , 1. By the similar method used in A.2.2, we have It is easy to show that the sequence {(2f r ) n − (δ 1 ) n } is increasing in n, under the assumption f r . 1 2 ; the sequence {n(k)} is non-increasing in k, i.e. n(k + 1) ≤ n(k). Then applying the above inequality with k replaced by k + n(k), one gets P S ((k þ n(k) þ n(k þ n(k)))t p ) P S ((k þ n(k))t p )((2f r ) n(kþn(k)) À d n(kþn(k)) 1 ) P S ((k þ n(k))t p )((2f r ) n(k) À d n(k) 1 ) P S (kt p )((2f r ) n(k) À d n(k) 1 ) 2 : Repeating this argument, we have a sequence of population approach to zero as (2f r ) n(k) À d n(k) 1 , 1. Since {(2f r ) n(k) − (δ 1 ) n(k) } is non-increasing in k, we can just consider the condition {(2f r ) n(k) − (δ 1 ) n(k) } < 1 when k is sufficiently large. Then, we can simplify the condition {(2f r ) n − (δ 1 ) n } < 1 with Therefore, part (a) is proven. Next we work with (b). Similar to the method used in the proof of A.2.3 with the assumption , t a , we obtain t À t p , p, a da dp: When δ 2 = 0, δ 1 + δ 3 = 1 and Although there is no concrete biological data in stem cell division types or specific regulation mechanisms to our knowledge, we discuss the selection of parameters in table 1 referring to both the biological observations [43][44][45][46] and the previous mathematical models [25][26][27]32,[34][35][36][37][38][39][40]54,55]. For the simulations we considered here, we assume that damage segregates symmetrically between stem cells and between TD cells, i.e. α 1 = 0.5, β 1 = 0.5, while the damage retention in asymmetric division is low, i.e. γ 1 = 0.25. Cell cycle progression speed and damage accumulation speeds v p , v a and u a . It is generally accepted that the population of mitotic cells is much smaller than that of TD cells, although it is very difficult to identify and count stem cells and progenitors [57][58][59][60]. Due to lack of experimental data, we refer to some existing modelling works for a reasonable range for the ratio of mitotic cells to TD cells. In the simulations of [37,55], populations of post-mitotic cells are about 4-10-fold that of mitotic cells.
On the other hand, mitotic cells are assumed to have a smaller death rate than post-mitotic cells [7]. However, there is no precise description of death rates due to two reasons: the death rate of stem cells is not observable, and the death rate of TD cells is tissue specific. In various modelling works, the assumed death rate of TD cells ranges from 10 −1 to 10 −4 of the total population [25,39]. Analysis and simulations in [25,39,40] suggest that a large death rate of TD cells results in instability of population. In our model, the death rates are measured in the following way: 0 v a S(t, p, a Ã ) dp P S (A 11) and r T ¼ where the numerators are population out-fluxes due to death and the denominators are the total populations. By properly choosing v p , v a , u a , we wish to achieve that P T /P S is no less than 5, r T is about 10 −4 , and r S is reasonably small, compared to r T . Here, we take v p = 0.2, v a = 0.06 and u a = 0.02. Maximum fractions of three types of cell division d 0 1 , d 0 2 and d 0 3 . It was observed that stem cells mainly undergo asymmetric renewal and differentiation when stem cell pool is expanding or when TD cells suffer from great loss, and enter into an inactively dividing stage and then mainly undergo asymmetric division at steady state [43][44][45][46]. Since the feedbacks we consider are negative regulations, d 0 1 provides the maximum fraction of symmetric renewal. By the property of the Hill function, δ 1 (P T ) can switch between large and small values in response to small and large populations of TD cells. But we observed that as we increase the initial d 0 1 , population overshoot and oscillations appear before convergence ( figure 10a). Moreover, the population ratio P T /P S decreases (from 5.664 to 3.622) and death rates of stem cell r S increases (from 1.191 × 10 −7 to 1.365 × 10 −4 ) as d 0 1 increases from 0.6 to 0.8. Hence, we take d 0 1 ¼ 0:6, d 0 2 ¼ 0:3 and d 0 3 ¼ 0:1. Regulation constants k T 1 , k T 2 and k T v . Many works suggest that regulation parameters k T i and k T v are closely related to the TD cell population at steady state [26,27,35,39]. These modelling works mainly focus on the blood system, which is one of the most well-studied systems and serves as a paradigm for understanding stem cells. The regulation parameter k in those models ranges from 10 −7 to 10 −10 , which corresponds to red blood cell population of size 10 7 ∼ 10 10 . However, some other works that do not specify the type of stem cells may choose some other ranges. For example, [ On the other hand, by choosing appropriate regulation parameters k T 1 , k T 2 and k T v in equations (4.3) and (4.4), asymmetrical division will predominate in three types of division of stem cells at steady state. In our simulations shown in figure 10b, we found that equal strengths of feedback on stem cell symmetric renewal δ 1 and cell cycle progression V p would result in oscillations. Stronger regulation on V p relative to δ 1 is needed to obtain stabilized population evolution (figure 10b). We also found that the relative strengths of regulation on symmetric renewal δ 1 and symmetric differentiation δ 2 play a role in determining the death rate of stem cells (r S = 2.408 × 10 −4 when (k T 1 , k T 2 , k T v ) ¼ (10 À8 , 10 À8 , 0:5 Â 10 À8 ); r S = 1.191 × 10 −7 when (k T 1 , k T 2 , k T v ) ¼ (10 À8 , 0:5 Â 10 À8 , 0:5 Â 10 À8 )). To obtain a reasonably small death rate of stem cells, stronger regulation on δ 2 is also needed. Hence, we take k T 1 ¼ 10 À8 , k T 2 ¼ 0:5 Â 10 À8 and k T v ¼ 0:5 Â 10 À8 . Hill exponent for feedbacks m T . The modelling work [40] suggests that a low Hill exponent in the feedback is needed to avoid oscillations, where the authors observed unstable steady states for Hill exponent m ≥ 3 through direct simulations in their ODE models. We also found similar phenomena in our model, where oscillations appear in the simulations with m T > 4. Thus we need m T ∈ {1, 2, 3}. To make the best estimation of m T , we also need to consider all biological observations. Our simulations show that when m T = 1 the death rate of stem cell is only about one third of that of TD cells (r S = 9.586 × 10 −5 when m T = 1; r S = 1.191 × 10 −7 when m T = 2), while when m T = 3 there is population overshoot before convergence to steady state although the death rate of stem cell decreases to 10 −15 (figure 10c). Hence we choose m T = 2 in our model.
After a large number of trials, we arrive at the best choice of the parameters in table 1 with consideration of all biological and mathematical reasons. Under these parameters, the resulting death rates are r S = 1.191 × 10 −7 and r T = 2.357 × 10 −4 . The populations at steady state are P S = 2.496 × 10 7 , P T = 1.414 × 10 8 with the ratio P T /P S = 5.664. But we need to point out that all our assumptions are pieced together through different biological and mathematical research, more precise estimations of parameters need real experimental data.    Regulation constants k S 1 , k S 2 and k S v , and Hill exponent for feedbacks m S . Similar to appendix A.3.1, due to the consideration of reasonable population ratio and death rate of stem cells, we still assume that the feedbacks from stem cell population on cell cycle progression V p and symmetric differentiation δ 2 are stronger than that on symmetric renewal δ 1 . That is, we intentionally choose smaller k S v , k S 2 , compared to k S 1 . And the Hill exponent m S is chosen to be 2.
To determine the magnitude of k S 1 , we simulate the population evolution with k S 1 ranging from 10 −6 to 10 −8 . We found that the relative magnitude of k T 1 and k S 1 has little influence on the population ratio. To effectively model the regulation from the stem cell population, 10 −7 turns out to be the best choice, according to our simulations under the choice of parameters in table 1. Hence, we take k S 1 ¼ 10 À7 , k S 2 ¼ 0:25 Â 10 À7 and k S v ¼ 0:5 Â 10 À7 . Constants for sigmoid regulation. According to experiments, mitotic cells inheriting more damage protein aggregates have an increased cell cycle length and tendency to differentiate [13,18]. To model these observations, we assume that the multipliers f (a) and g(a) decrease in a and have sigmoid shape. Parameters a i and b i in the multipliers determine their maximum and minimum strengths, a 0 i locates the threshold in damage of the transition and k a i controls the stiffness of the transition. Due to lack of experimental data or previous modelling works to refer to, we assume f (a) and g(a) take the same form and range between [0.4, 1.1] with a transition at a 0 1 ¼ 0:75. Hence, we take a 1 = a 2 = 1.1, b 1 = b 2 = −0.7, k a 1 ¼ k a 2 ¼ 20 and a 0 1 ¼ 0:75. With a large number of simulations, we present the best choice of the parameters in table 2 and conclude that this set of parameters is our best estimation. Under these parameters, the resulting death rates are r S = 2.560 × 10 −7 , r T = 2.314 × 10 −4 for stem cells and TD cells, respectively. The populations at steady state are P S = 9.593 × 10 6 , P T = 6.709 × 10 7 with the ratio P T /P S = 6.994.

A.4. Estimation of populations in the models with feedbacks from TD cells
First, we recall the assumptions: a Ã ≤ a c , V a = v a , U p = u p and U a = u a are constants, and δ i and V p have feedbacks from TD cells described by Under the above assumptions, we integrate (2.1) and (2.2) in p and a to get dP S dt ¼ Àv a ð p Ã 0 S(t, p, a Ã ) À S(t, p, 0) dp À V p (P T ) ð a Ã 0 S(t, p Ã , a) À S(t, 0, a) da (A 15) and dP T dt ¼ Àu a S(t, p, a Ã ) dp À V p (P T ) T(t, p, a c ) dp þ V p (P T ) Making a change of variable in a and noticing that a Ã ≤ a c , we obtain a system of ODEs, called system I: S(t, p, a Ã ) dp þ V p (P T )(d 1 (P T ) À d 2 (P T )) ð a Ã 0 S(t, p Ã , a) da (A 17) and dP T dt ¼ Àu a Then under the same initial condition, the solution to system I is not larger than that to system II, by comparison theorems for ODEs. Now we assume that both systems I and II have steady-state solutions, denoted by ( P S , P T ) and (P S ,P T ), respectively. ThenP T can be obtained by solving d 1 (P T ) À d 2 (P T ) ¼ 0: Assuming appropriate values of parameters to guarantee solvability of the last equation, we have the upper-bound estimate : (A 20) Next, we continue to estimate P S . First, we observe that in the steady state both S(t, p, a) and T(t, p, a) are constant along characteristic curves and we refer the reader to figure 11 for an illustration of the situation. Since in system II, there is no outflux of stem cells on the boundary a = a Ã , we havẽ S(t, p Ã , a) da: On the other hand, by the intermediate value theorem for integrals, we have, for some h ∈ [h 1 , h 2 ] that Da and h 2 ¼ a c À h 1 1 À v 1 Da with η 1 = min{β 1 , β 2 , γ 2 }, η 2 = max{β 1 , β 2 , γ 2 } and Da ¼ ( p Ã =V p ( P T ))v a . Using thatP T is a steady-state population, we have by (A 17) u aPT ¼ hV p (P T ) ð a Ã 0 S(t, p Ã , a) da ¼ hV p (P T )P S p Ã and thus the estimate P S P S ¼ u a p Ã hV p (P T )P T u a p Ã V p ( P T )h 1P T , ( In the following simulation, we use third-order WENO scheme and third-order TVD Runge-Kutta time integrator. As U p , U a , V p , V a are all positive, the numerical approximation W i,j to the exact solution W( p i , a j , t) satisfies the following ODE system: whereF iþ1=2,j is called numerical flux, the design of which is the key ingredient to a successful scheme. For the third-order WENO scheme, the numerical fluxF iþ1=2,j is defined as follows: iþ1=2,j , ( whereF (m) iþ1=2,j for m = 1, 2, are the two second-order accurate fluxes on two different stencils given bŷ F (1) iþ1=2,j ¼ À iþ1=2,j ¼ The nonlinear weights ω m are given by