## Abstract

Stem cells (SCs) perform the task of maintaining tissue homeostasis by both self-renewal and differentiation. While it has been argued that SCs divide asymmetrically, there is also evidence that SCs undergo symmetric division. Symmetric SC division has been speculated to be key for expanding cell numbers in development and regeneration after injury. However, it might lead to uncontrolled growth and malignancies such as cancer. In order to explore the role of symmetric SC division, we propose a mathematical model of the effect of symmetric SC division on the robustness of a population regulated by a serial differentiation cascade and we show that this may lead to extinction of such population. We examine how the extinction likelihood depends on defining characteristics of the population such as the number of intermediate cell compartments. We show that longer differentiation cascades are more prone to extinction than systems with less intermediate compartments. Furthermore, we have analysed the possibility of mixed symmetric and asymmetric cell division against invasions by mutant invaders in order to find optimal architecture. Our results show that more robust populations are those with unfrequent symmetric behaviour.

## 1. Introduction

Tissues in higher organisms such as mammals maintain homeostasis by means of cascades of serial cell differentiation [1]. These cascades are composed of hierarchically organized compartments of different cell types. The cornerstone of this system is stem cells (SCs) which sustain the tissue by producing both more SCs, a property known as self-renewal and lineage-specific differentiated cells. The serial differentiation cascade proceeds through a number of intermediate stages or compartments of transient amplifying cells (TACs) until it terminates with the compartment of fully-differentiated mature cells (MCs).

There has been much discussion around the issue of how SCs accomplish the feat of undergoing both self-renewal and differentiation [2–6]. One strategy by which SCs can perform these two tasks is asymmetric cell division, whereby, upon SC division, one daughter cell retains the cellular identity of its mother, whereas its sister differentiates. Asymmetric cell division seems to proceed through segregation of the material that determines cell fate so that, when the SC divides, one of the daughters keeps the SC-fate determinants. Lineage-specific determinants are thus passed onto its sister cell [2,7].

Although asymmetric SC division is a simple and elegant solution to the SC puzzle, there has been much discussion and controversy around what is perceived to be a fundamental flaw of this model, namely, under asymmetric division conditions the SC compartment is unable to expand [4,5]. The argument runs that in situations such as injury, where the rate of cellular turnover must be largely increased in order to regenerate the affected tissue, the size of the SC pool should undergo a marked expansion. This hints that asymmetric division might not be the full solution.

Experimental evidence points to the fact that SCs in some tissues, such as the central neural system [8] and the epidermis [9], can divide symmetrically. Symmetric SC division consists of both daughter cells sharing the same fate: SC or differentiated. In this way, the number of cells is regulated by the frequencies of (symmetric) divisions producing SCs or differentiated cells. As this division modality allows for expansion of the SC compartment, an alternative model arises [4] where most SCs can divide both symmetrically and asymmetrically with the balance between these two proliferation modes controlled by environmental signals to produce the appropriate number of SCs and differentiated cells. However, it has been recognized that, while symmetric division confers the capability of enhanced growth and increased regenerative capacity, it also increases the likelihood of cancer. This view is supported by findings according to which the cellular machinery involved in asymmetric division has been evolutionary conserved with a role in tumour suppression, the gene products that promote symmetric SC division also function as oncogenes [4].

The aim of this paper is to address some of the issues arising from the consideration of symmetric SC division in serial differentiation cascades, in particular those regarding the evolution of uncontrolled growth and the risk of cancer. To do so, we consider a stochastic mathematical model of the serial differentiation cascade and use it to address two problems (i) the long-term viability (stability) of a population of cells maintained by a differentiation cascade whose SCs undergo symmetric division only and (ii) the likelihood of invasion of a mutant population with increased symmetric SC division in competition with a population with asymmetric SC division (robustness). Regarding the former, we show that serial differentiation cascades with symmetric SC differentiation exhibit an intrinsic instability mechanism, which induces extinctions with high likelihood. Thus, these cascades seem to have an in-built fail-safe mechanism against excessive symmetric SC division. The competition model allows us to estimate *optimality* conditions, in the sense of maximum robustness, for the rate of symmetric division in terms of maximizing its chance to take over another population with less symmetric division capability while minimizing its chances of undergoing extinction.

There is large body of literature examining the properties of differentiation cascades, in particular in the haematopoietic system [10–12] and in relation to several kinds of leukaemias and haematological diseases [13–16]. Marciniak *et al*. [12] have studied a deterministic model in which they explore the role of regulation in a system with asymmetric SC division. They have found that for their model to exhibit efficient repopulation, regulation by environmental signals must occur at the level of the fraction of self-renewal. Mackey and co-workers [13–15] have done extensive modelling and analysis of a number of haematological diseases, in which circulating cellular blood components exhibit oscillatory behaviour. By means of delay differential equation models and exhaustive parameter sensitivity analysis, they have proposed which processes in the haematopoietic differentiation cascade are more likely to be involved in the development of disease. They have also found that delays in the differentiation cascade are fundamentally linked to the onset of the oscillations, a mechanism akin to the extinction mechanism we explore here. Dingli *et al.* [17] have studied a model of the architecture and dynamics of haematopoiesis to provide estimates of the number of compartments, the size of each compartment and the corresponding replication rates.

There have also been a number of stochastic models of serial differentiation cascades that tackle evolutionary issues. The model proposed by Pepper *et al*. [18] deals with a number questions arising when considering which evolutionary purpose can such a costly system (from the energetic point of view) possibly serve. The conclusion of their study is that differentiation cascades provide intrinsic protection against somatic mutations leading to malignancies such as cancer. Evolution towards malignant behaviour and cancer has been studied, for example, by Sun & Komarova [19] and Rodriguez-Brenes *et al*. [20,21], where the mechanisms of regulation of normal tissues are studied in detail to ascertain how can they be subverted in order to produce cancer. Other related work is the symmetric SC division models proposed by Dingli *et al*. [7], where a Moran process is proposed to study the competition between a population of normal SCs and a population of cancer SCs. They observed that a larger duplication rate in the invader (cancer) SCs increases the probability of invasion.

Michor and co-workers [16,22–25] have modelled the dynamics of myleoid leukaemia using models consisting of hierarchical differentiation cascades. It has been found that cancer initiation requires only a single mutation. Their models suggest that imatinib is a potent inhibitor of the production of differentiated leukaemic cells, but does not deplete leukaemic SCs. Lenaerts *et al.* [26] have shown that although imatinib does not eradicate chronic myleoid leukaemic SCs, because of the stochastic nature of haematopoiesis, leukaemic SCs undergo extinction. This result shows the importance of stochastic effects on the dynamics of hierarchical cell populations.

Zhao *et al*. [27] have proposed a linear model for colonic epithelial cells that explicitly takes into account the proliferation kinetics of a cell as a function cell position within the crypt. They have observe that those cells with mitotic activities concentrated near the SCs delay the rate of mutation accumulation in colonic SCs.

Mutations in hierarchically structured populations have been previously studied by Werner *et al*. [28], they modelled a general hierarchically organized multi-compartment tissue, allowing any number of mutations in a cell. They showed that these tissues strongly suppress the cells that are carrying multiple mutations. Traulsen *et al.* [29] have studied the effect of different kinds of mutations and showed that although neutral mutations often lead to clonal extinction, disease is still possible, and in this case, it may become difficult to eliminate neutral mutations with therapy.

Also related to the context of this paper, a mechanism for stochastic evolutionary escape and its consequences for drug resistance based on the presence of a quiescent, drug-resistant population has been recently described [30]. As quiescent cancer SCs have been recently found in intestinal crypts [31], this escape mechanism can be relevant to ascertaining the efficacy of therapies targeting cancer SCs.

The organization of the paper is as follows. Section 2 is devoted to presenting our model of a differentiation cascade with symmetric SC division in detail and discussing its underlying hypothesis. In §3, we present simulation results for this model and discuss the emergence of a mechanism of extinction induced by symmetric SC division. In §4, we analyse the robustness of the different division strategies by formulating a model for the competition between populations using different such strategies and extract conclusions regarding their optimality. Finally, in §5, we summarize and discuss our results as well as directions for future work.

## 2. Stochastic modelling of a serial differentiation cascade with symmetric stem cell division

In this section, we proceed to present and discuss our stochastic model of a hierarchically organized cell population. Our aim is to formulate a regulated stochastic model which, by incorporating symmetric SC division, can lead to changes in the size of the SC compartment. Regulation is achieved by means of a cytokine, which is produced by the MC compartment, which controls the rate of apoptosis of the SCs.

The stochastic approach is justified by the small number of SCs, typically present in their host organism. For example, according to Morrison & Kimble [4], in the haematopoietic system only 0.003% of the cells are SCs. Given these small numbers in the SC compartment, we expect random effects to be important if not dominant.

### 2.1. Model description

We consider a compartmental model of the serial differentiation cascade with symmetric SC division. As shown in figure 1, each compartment corresponds to a stage in the cascade: from the SC compartment to the MC compartment with a number of intermediate TAC compartments. We further introduce in our model the secretion of a cytokine by the MCs which regulates the size of the SC compartment. We study the dynamics of the population of cells in each compartment. In particular, we examine the behaviour of a differentiation cascade of length *n* + 1, with SCs, *n* − 1 partially differentiated, TACs, and MCs: *x*_{0} represents the number of SCs, *x _{i}* with

*i*= 1, … ,

*n*− 1, the number of TACs at differentiation stage

*i*and

*x*, the number of MCs. Furthermore,

_{n}*x*

_{n}_{+}

_{1}represents the amount of a cytokine that regulates the size of the SC compartment through a negative feedback [13,14,32,33]. This cytokine is secreted at a rate that is proportional to the number of fully-MCs. Therefore, this negative feedback accounts for the MC regulation of the size of the SC compartment which has been incorporated, for example, in models of haematopoiesis [12].

Our stochastic model is formulated in terms of the corresponding master equation [34,35]
2.1where *X* = (*x*_{0}, … , *x _{n}*,

*x*

_{n}_{+}

_{1}) and

*P*(

*X*,

*t*) are the probability of the population vector to have value

*X*at time

*t*. For full specification of our stochastic models, we need to prescribe the transition rates,

*W*(

_{j}*X*,

*t*), corresponding to the probability per unit time of the process

*j*. These rates are modelled in terms of the law of mass action [36].

*r*is the change in

_{j}*X*when elementary process

*j*occurs, i.e.

*X*(

*t*+ Δ

*t*) =

*X*(

*t*) +

*r*with probability

_{j}*W*Δ

_{j}*t*.

We will consider two different scenarios schematically represented in figure 1. The first one, corresponding to figure 1*a* and which will be referred to as the well-stirred scenario, the cytokine that controls the size of SC compartment is assumed to be instantaneously transported to it. By contrast, the second scenario, corresponding to figure 1*b* and to be referred to as the delayed scenario, the cytokine requires a finite amount of time to reach the SC compartment. We have modelled this by introducing a compartment where, upon secretion, the cytokine is transported to at a certain rate. The elementary processes involved in our model are:

— SCs divide symmetrically, so that they can undergo:

(1) Self-renewal. Both daughter cells share the SC fate of the mother with transition rate

*W*_{1}(table 1).(2) Differentiation. Both daughter cells differentiate to first-stage TAC with transition rate

*W*_{2}(table 1).(3) Apoptosis with cytokine-regulated transition rate

*W*_{3}if we are considering the well-stirred scenario or*W*_{7}_{+}_{2n}, otherwise (table 1).

— Stage-

*i*,*i*= 1,*…*,*n*− 1, TACs can undergo:— MCs can undergo:

— Cytokines can undergo:

(1) Clearance with transition rate

*W*_{4}_{+}_{2n}if we are considering the well-stirred scenario, or*W*_{6}_{+}_{2n}, otherwise (table 1).(2) Transport to the intermediate compartment (

*x*_{n}_{+}_{2}in figure 1*b*) with transition rate*W*_{5}_{+}_{2n}if we are considering the delayed scenario (table 1).

We will analyse this system by means of numerical Monte Carlo simulations using the Gillespie stochastic simulation algorithm [36,37]. The mean-field limit of the model also provides useful information. For the model described by equation (2.1) and table 1, corresponding to the well-stirred scenario, the mean-field model is given by the following system of ODEs:
2.2
2.3
2.4
2.5which has two fixed points. The trivial (unstable) equilibrium (0, 0, … , 0), which corresponds to the extinction of the system, a positive equilibrium fixed point, *P _{S}*
The stability of this fixed point is studied in appendix A. We will fix the number of

*x*cells at the metastable state. We will denote as , the number of cells in the

_{n}*i*compartment at the metastable state.

Our model assumes that the signal is produced by the MCs, we make this assumption instead of considering that the signal is produced by, for example, the TACs in position *k*. In this case, the signal-producing rate *s* should be of the form:
2.7

The signal-producing rate is necessary to the cells in each compartment in order to keep the amount of signal necessary to control the population. We see that the only efficient possibility is that the signal is generated by the MCs. However, this possibility of TAC cells producing signal is not out of our study. If the signal is produced by the *x _{k}* cells, the new system is equivalent to a system described above with length

*k*taking as a new (the death rate of the new MCs) the value of

*λ*+

_{k}*d*, since the

_{k}*x*

_{k}_{+}

_{1}, …

*x*cells has no effect on the system.

_{n}In the scenario (b), that is, the cytokine is not transported immediately to the SCs compartment, we consider another compartment, *x _{n}*

_{+}

_{2}, and the reactions and , changing by . In §3, we will see the effect of this change in function of the transport rate

*s*

_{2}.

## 3. Extinction induced by symmetric stem cell division

We now proceed to analyse the stochastic model described in §2.1 and describe an extinction mechanism intrinsic to the structure of a serial differentiation cascade with symmetric cell division. We interpret this mechanism as an anti-cancer mechanism inherent to the hierarchical structure of serial differentiation cascades.

### 3.1. Delay-induced extinction in serial differentiation cascades with symmetric stem cell division

In order to proceed with our analysis, we perform numerical simulations of our model by means of the Gillespie stochastic simulation algorithm [36,37]. This algorithm is based on an exact reformulation of the stochastic process described by master equation (2.1) that allows us to generate exact sample paths or realizations of the underlying process. An example of such a realization, which illustrates the extinction mechanism in our model with symmetric SC division, is illustrated in figure 2.

Figure 2 shows simulation results where the population has been initially set to be equal to the (stable) positive equilibrium of the mean-field limit equations (2.2)–(2.5). This particular example is illustrative of the mechanism through which extinction occurs. Initially, a fluctuation in the system induces an increase in the number of SCs. This increase in SC population propagates through the differentiation cascade and it eventually reaches the MC compartment which starts expanding. Such expansion in MC numbers induces an increase in cytokine concentration, which, in turn, leads to decay in the size of SC compartment. As a consequence, oscillatory behaviour ensues. However, since the system, specially the SC compartment, is subjected to fluctuations and these are amplified by the length of the differentiation cascade (an increase of one cell in the size of the SC compartment induces an increase *O*(2* ^{n}*) in the MC compartment), the amplitude of the oscillations grow until the SC compartment size hits the absorbing barrier at which point the population undergoes extinction.

As delays appear to be instrumental for this extinction pathway, we are interested in the effects of altering the balance between the *x _{n}*

_{+}

_{1}/

*x*cells' lifetime and the replenish rate via

_{n}*x*

_{0}on the average extinction time. To carry on this analysis, we compute the extinction time as a function of the SC proliferation rate,

*p*, and the death rate of the MCs,

*λ*. Our results are shown in figure 3.

_{n}Figure 3 shows that there is a region in the *λ _{n}* −

*p*-plane where the average extinction time,

*T*

_{E}, is such that

*T*

_{E}> 10

^{4}days. For all other values of these two parameters, the waiting time for extinction to occur becomes much shorter. It is relatively straightforward to see that the effect of these two parameters within the region of the parameter space where

*T*

_{E}> 10

^{4}is to cancel the effect of the delay between the SC and MC compartments owing to the length of the differentiation cascade. Large values of the MCs death rate,

*λ*, stabilizes the system via a dissipative-like effect, namely, it impairs the ability of the MC compartment to expand in response to fluctuations in the SC compartment size. It also weakens the effect of the feedback between the MC and SC compartments, as it effectively reduces the concentration of cytokine.

_{n}The effect of varying the value of the self-renewal rate, *p*, is slightly different. When , the SC compartment cannot grow, on the contrary, under the effect of SC apoptosis and random fluctuations, it will dwindle and the population will become extinct. As *p* grows, its effect is strongly coupled to *λ _{n}*. If

*λ*is small, i.e. the MC life expectancy is increased, the effects of fluctuations and delay are reinforced, and the population undergoes rapid extinction. If, on the contrary,

_{n}*λ*is large, we observe two regimes: a regime of intermediate values of

_{n}*p*where the effects of fluctuations and delay are cancelled by the rapid death of the MCs (see above), and a regime of high values of

*p*where MC death is not fast enough to compensate for the increased effects of fluctuations in the SC compartment which take over and produce short survival times.

Another model parameter that has an obvious effect on the delay is the length of the differentiation cascade, *n*. Longer chains promote the effect of the delay and therefore favours extinctions. This is illustrated in figure 4, where we show sample paths of our model for different values of *n*. We see that indeed longer differentiation cascades destabilize the system and favour extinction. This is further explored in figure 5, where we show the extinction probability within the time window of time of duration (0, *T*_{E}), *P*_{E}(*T*_{E}), as a function of *n*. Figure 5 shows a sharp increase in this quantity as *n* increases.

### 3.2. Delayed cytokine scenario

In this section, we consider the delayed cytokine scenario where:

— the secreted cytokine,

*x*_{n}_{+}_{1}, is transported to*x*_{n}_{+}_{2}, that is, , with transition rate*W*_{5}_{+}_{2n}=*s*_{2}*x*_{n}_{+}_{1},— the cytokine in

*x*_{n}_{+}_{2}is degraded, that is, , with transition rate*W*_{6}_{+}_{2n}=*λ*_{n}_{+}_{2}*x*_{n}_{+}_{1}and— the size of the SC compartment is controlled by

*x*_{n}_{+}_{2}, that is , with transition rate*W*_{7}_{+}_{2n}=*λ*_{0}*x*_{0}*x*_{n}_{+}_{2.}

We start by analysing the mean-field behaviour of the system
3.1
3.2
3.3
3.4
3.5which implies that the secretion rate, *s*, is now determined by
3.6and the positive equilibrium is of the form:
Note that if *λ _{n}*

_{+}

_{2}=

*λ*

_{n}_{+}

_{1}, we recover the previous model when

*s*

_{2}→

*∞*.

We are interested in how the dependence of the behaviour of population on the cytokine delay, which is controlled by *s*_{2}. Figure 6 shows the extinction time as a function of *s*_{2}. In figure 7, we have plotted *s* as a function of *s*_{2}, as given by equation (3.6). Taken together, these two results show that a decrease in *s*_{2} contributes to destabilize the population, as this contributes to increase the delay between cytokine secretion and SC response. Moreover, more cytokine needs to be produced to maintain the SC compartment size.

## 4. Robustness of populations with asymmetric stem cell division

So far, we have focused on the study of differentiation cascades with symmetric SC division. We have uncovered that symmetric SC division may induce extinction of the population, especially in long differentiation cascades where the delays in the regulation of SCs by MCs are more important. This result points out that utilization of symmetric SC alone may not be an optimal strategy.

In normal situations, SC division is commonly believed to be mostly asymmetric with symmetric division being a rare event [2]. However, it has also been put forward that symmetric SC division could also play an important role, especially in situations, such as wound repair, where the SC compartment must be allowed to expand in order to compensate for an increased demand in cellular turnover [9]. In the view of our results regarding the induction of instabilities by symmetric SC division, one could ask the question of how frequent symmetric SC division can be before anomalies in the population start to appear.

The aim of this section is to address these two issues, namely, whether symmetric SC division is an optimal strategy and the frequency of symmetric SC division, as well as ascertaining which are the generic features of an optim, i.e. maximally robust, differentiation cascade.

### 4.1. Stochastic modelling of a differentiation cascade with asymmetric stem cell division

Our first step is to present our stochastic model for a differentiation cascade with asymmetric SC division. We formulate our model in terms of the transition rates that go into the corresponding master equation (2.1), which is shown in table 2. Model formulation is very similar to the one for the differentiation cascade with symmetric SC division: we consider the same *n* + 2 different cell and chemical species, i.e. SC, *n* − 1 TACs, MCs and the MC-secreted cytokine.

The dynamics of this system is essentially the same as in the symmetric SC division case (table 1) except for the SC cell dynamics. Following Knoblich [2], we assume that SCs divide asymmetrically with occasional symmetric cell division being a rare event. We thus consider that SCs divide asymmetrically with probability 1 − *e*. We further consider than the rate of asymmetric division is regulated (inhibited) by a cytokine produced by the MCs [12]. These two factors are summarized in the transition rate *W*_{1}, table 2. With probability *e*, SCs behave as a symmetrically dividing cell, whose behaviour is characterized by transition rates *W*_{1}, *W*_{2} and *W*_{3} (table 2), i.e. an SC that behaves as a symmetrically dividing cell can (i) self-renew, (ii) differentiate or (iii) undergo apoptosis, with probability rates *W*_{1}, *W*_{2} and *W*_{3}, respectively. In our study, the parameters *a* and *b* will be chosen to have the same division rate as in the symmetric division model at the metastable state, that is, .

The mean-field behaviour of the model is described by the following system of law-of-mass action ODEs:
4.1
4.2
4.3
4.4
4.5where *y _{i}* is the number of cells in the compartment

*i*. The variables and parameters are described in table 2. Equations (4.1)–(4.5) have a positive equilibrium given by 4.6The stability of the fixed point is studied in appendix A. We will denote as the number of cells in the

*i*compartment at the metastable state. We will study the system in the region of the parameter space for which the positive equilibrium is a stable attractor.

Figure 8 shows the average stationary size of the SC compartment as predicted by equations (4.1)–(4.6). We observe that, while an increase in the size of the SC compartment is obtained as *e* is allowed to grow, the magnitude of such increase, for the chosen set of parameter values, is of the order of 10%.

Regarding the behaviour of the stochastic model, the main difference with respect to the symmetric cell division model (table 1) involves the lack of delay-induced extinction behaviour exhibited by the differentiation cascade model with pure symmetric SC division. This behaviour is illustrated in figure 9 in which two typical sample paths are shown for different values of the probability of symmetric division, *e*. If the probability of symmetric division, *e*, is small, the variation in size of the SC compartment is bounded, which renders the system resilient to the effect of the delay-induced oscillations that, in the pure symmetric SC division system, eventually lead to extinction of SC compartment followed by the extinction of the whole population. In figure 10, we have plotted the extinction probability of a differentiation with asymmetric SC division as a function of *e* and show that, indeed, differentiation cascades with small probability of symmetric division are impervious to the delay-induced extinction mechanism.

### 4.2. Competition between populations

We now move on to the study of the robustness of the asymmetric population, that is, whether it can withstand attempts of invasion by other populations using different strategies. In what follows, such strategies will refer to variations of two parameters: the probability of symmetric division, *e* and the rate of symmetric self-renewal *p*. Our aim is to ascertain criteria of evolutionary optimality based on our results regarding the outcome of competition between populations characterized by different values of these parameters.

#### 4.2.1. Asymmetric versus symmetric population

We start by studying the competition between a resident population generated by differentiation cascade with asymmetric SC division, hereafter referred to as *A*, and a symmetric-SC-division invader, referred to as *S*, whose population dynamics is given by the model described in §2.

*Mean-field behaviour.* In order to proceed further, we start by giving some results regarding the behaviour of the mean-field limit using qualitative theory of ODEs. For simplicity, we will assume that both populations of MCs are producing the same cytokine.

Let *y*_{0} (*x*_{0}) be the size of the SC compartment of the asymmetric (symmetric) division population. The mean-field behaviour of the model of competition between the populations *A* and *S* is described by the following system of law-of-mass action ODEs:
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15

Let the quantities *r* and *r’* be defined as
4.16i.e. *r* and *r’* are the ratios between the net SC growth rate and the SC differentiation rate for each of the competing populations, which coincide with steady-state value of the mean-field cytokine concentration for each of the populations.

From equation (4.7), we see that , i.e. *y*_{0} grows in time, if . On the contrary, *y*_{0} decreases as time progresses if . Similarly, equation (4.11) implies that *x*_{0} is increasing with time when and decreasing when . Moreover, according to equations (4.7) and (4.11), when *t* → *∞*, the stationary level of cytokine is given by *x _{n}*

_{+}

_{1}→ max (

*r*,

*r’*). We thus consider three possible scenarios for the mean-field limit:

—

*r*>*r’*. The*S*population becomes extinct: as*t*→*∞, x*_{n}_{+}_{1}→*r*which implies that, as*t*→*∞*, and with the equality in the latter holding for*x*_{0}= 0 only.—

*r*=*r’*. Both populations coexist: as*t*→*∞*,*x*_{n}_{+}_{1}→*r*=*r’*which implies that, as*t*→*∞*, both and hold.—

*r*<*r’*. By the same reasoning as in the first case, the*A*population becomes extinct.

Regarding the behaviour of the stochastic system, there are some important differences with respect to its mean-field counterpart. In figure 11, we show that the probability of invasion as a function of the self-renewal rate, *p’*, and death rate , of the invader SCs. Instead of the three scenarios that characterize the mean-field limit, we have

—

*r*≥*r’*. In this case, extinction of*S*is the most likely scenario: the population*S*is forced to evolve under higher on-average cytokine levels, which induce upregulation of apoptosis of the*S*-SCs. As*r*′ →*r*, this competitive advantage of*A*over*S*weakens. However, at that point, the delay-induced extinction mechanism for*S*takes over. As we are considering*e*to be small, delays have little effect on the population*A*, which leads to eventual extinction of*S*. This behaviour is illustrated in figure 11, which shows that the probability of*S*to invade*A*vanishes in this regime.—

*r*<*r’*. Under these conditions, the population*A*is the one that evolves under*excess of cytokine*conditions, as the average steady-state concentration of cytokine is higher than the one the*A*population would have in the absence of the population*S*. As a consequence, the probability of*S*invading*A*is now positive, as shown in figure 11. We also observe that the probability of invasion as a function of*r′*reaches a maximum value and then starts decreasing. This behaviour is again due to delay-induced extinctions: as the relative net growth rate of the*S*-SCs increases the probability of delay-induced extinction of the*S*population also increases (as shown in figure 3).

In order to further explore the behaviour of the stochastic competition model in the *r* < *r’* regime, we have done simulations to compute the probability of invasion and the probability of extinction of both populations as a function of *e*, i.e. the probability of asymmetric behaviour of *A*-SCs, and the rate of (symmetric) self-renewal of *S*-SCs. The results are shown in figure 12. We observed (figure 12*a*) that there is a threshold value of *e* below which the resident *A*-type population cannot be invaded by an *S*-type invader, regardless of the value of *p′*. If the *A*-type resident has an *e* above this threshold, then two outcomes are possible: invasion, for moderately high values of *p′*, or extinction of the symmetric invader (see also figure 11), for even higher values of *p′*. Furthermore, as *e* grows, figure 12*b* shows that the more probable outcome is both populations be extinct. This is owing to the fact that for larger values of *e* the symmetric component becomes dominant in population *A* as well as *S*, and consequently the *A*-population starts to becoming affected by the effects of the corresponding delays. The probability of long-term coexistence is negligible, according to our simulations.

#### 4.2.2 Optimal differentiation cascades

Finally, we tackle the issue of examining the optimality properties of differentiation cascades. In this context, we define optimality of a differentiation cascade in terms of its resilience to invasion, namely, optimal differentiation cascades are those which are maximally resilient to invasion.

In order to address this issue, we again consider the competition between two populations generated by differentiation cascades characterized by different values of the probability of symmetric behaviour, *e*, and the self-renewal rate, *p*. We study the ability of the resident population to withstand invasion as a function of the values of these parameters for both populations. We have chosen these two particular quantities as the focus of our analysis, as they are directly related to the ability for adaptive behaviour of differentiation cascades, which has been argued to be at the root of symmetric cell division [4].

We start by analysing the optimal strategy regarding the probability of symmetric behaviour, *e*. To this end, we have run simulations of the competition between a resident with probability of symmetric behaviour, *e*, and an invader with probability of symmetric behaviour, *e′*. All the other parameters are taken to be equal for both populations. Results are shown in figure 13.

In the light of our results regarding delay-induced extinction for higher values of *e*, i.e. when the symmetric part of the dynamics dominates, we expect that there exists a trade-off between resilience to invasion and long-term survival and elevated frequency of SC symmetric behaviour (which favours adaptive responses). This expectation is confirmed by the simulations results shown in figure 13*a*, which show the probability of invasion. We can observe that resident populations with are very stable and resilient to invasion. As *e* increases, so that a faster and more efficient adaptive response be possible, an invader with *e’* ≤ *e* is likely to wipe out the resident population (figure 13*a*).

We have also investigated the effect of the rate of self-renewal on the resilience against invasion of an *A*-type population. Our results are reported in figure 14 where we have plotted the probability of invasion (figure 14*a*), the probability of extinction (figure 14*b*) and the probability of coexistence after a time *T _{F}* has elapsed (figure 14

*c*). These simulations have been done for

*e*= 0.1 and

*e’*= 0.05, which correspond to the region of parameter values where, according to figure 13, the probability of invasion is positive. The results shown in figure 14 indicate that if

*p*>

*p’*, i.e. if the rate self-renewal of the resident is larger than that of the invader, the resident population is resilient to invasion.

Taken together, the results shown in figures 13 and 14 show that the optimal trade-off between resilience to invasion and capability for adaptive behaviour is realized by populations with small probability of symmetric behaviour, *e*, and large value of the self-renewal rate.

## 5. Conclusion

We have presented and analysed stochastic models of hierarchical populations governed by differentiation cascades to study the effect of symmetric SC self-renewal on their long-time stability properties as well as their resilience to invasion. The behaviour of these models has given us some insight on intrinsic protection mechanisms, in particular against invasions by malignant, potentially tumorigenic, populations.

Our model of a hierarchically organized tissue consists of a stochastic compartmental model, where each compartment corresponds to a differentiation stage, with negative feedback regulating the size of the SC compartment. This feedback is assumed to be mediated by a cytokine whose secretion rate is, in turn, modulated by the number of fully-MCs. The general framework we have used in our model formulation, i.e. in terms of cellular compartments with feedback between the SC and the MC compartments, has been used in a wide variety of studies, ranging from general models of robustness against mutations leading to cancer [18–21] to more specific models of particular tissues such as the haematopoietic system [10–14], epidermis [38] or neurogenesis [39]. Our model formulation is therefore rather generic and the resulting model should be applicable to a wide variety of tissues where hierarchical structure is present.

However, the issue of the cytokine-mediated feedback has more supporting evidence in the haematopoietic system, where several of these cytokines have been identified which, upon variation in the size of the differentiated cell compartment, regulate the size of the SC compartment by means of negative feedback loops. For example, in the production of white blood cells, granulocyte colony-stimulating factor (G-CSF) has been demonstrated to mediate negative feedback [13]. Another example is thrombopoietin, which is involved in negative feedback regulation in the production of platelets [13]. Erythropoiesis, i.e. the production of red blood cells, is yet another example in which negative feedback is mediated by a cytokine, namely, erythropoietin [32,33].

In order to investigate the effects of symmetric self-renewal on the robustness of hierarchical populations, we have first formulated a model in which SCs undergo only symmetric self-renewal. Regarding this model, we have identified a new extinction mechanism in which the coupling between fluctuations that increase the pool SCs and the delay in SC regulation leads to the extinction of the SC population and, therefore, to the eventual extinction of the whole: fluctuations in the size of the SC compartment propagate along and are amplified by the differentiation cascade in an amount that typically scales as 2* ^{n}*, where

*n*is the length of the cascade. As the regulation of SC proliferation by the fully-differentiated cells has a delay proportional to

*n*, the negative feedback between both populations may lead to extinction of the SCs and therefore of the whole population.

The extinction time, *T*_{E}, in such symmetric populations has been shown to be sensitive to variations of the death rate of the MCs, *λ _{n}*, and the rate of self-renewal,

*p*. We have shown that increasing

*λ*leads to longer extinction times, whereas populations with larger values of the self-renewal rate

_{n}*p*exhibit shorter extinction times. Biologically, these results imply that symmetric SC self-renewal is more likely to be observed in those with shorter delays between the SC and MCs and those where the lifetime of the MCs is shorter and therefore larger cellular turnover is necessary to maintain tissue homeostasis.

We have also explored the possibility that the cytokine is secreted in response to variations in the size of some of the TAC compartments. We observe that the behaviour is similar to that of shorter differentiation cascades: if the cytokine is produced at a rate proportional to the number of cells in compartment *i*, the behaviour of the SCs is, for all intents and purposes, independent of compartments *i* + 1, … *m* which therefore do not affect the system. We have also seen that if the secretion rate of the cytokine were proportional to the number of some of the TACs, then the corresponding secretion rate would have to be much larger to maintain the same population of SCs, as the number of TACs is smaller than the number of MCs (results not shown). These two results suggest that it is more likely that the cytokine is produced in response to variations in the size of the MCs compartment.

In order to proceed further, we have turned our attention to a model that corresponds to a more biologically realistic situation: a stochastic model of differentiation cascades where SCs can proliferate both symmetrically or asymmetrically. In our model, the choice between these two division modes is a random event. Here, we have focused on the robustness of the population generated by a mixed (i.e. with both symmetric and asymmetric SC division) differentiation cascade in terms of its resilience against invasion (for example, invasion by a potential tumorigenic differentiation cascade). In particular, we have focused on optimal strategies for a differentiation cascade to be robust against invasions as defined by the values of two parameters: the probability of symmetric behaviour, *e*, and the rate of symmetric self-renewal, *p*. Our results show that there exists a trade-off between robustness and the capability of the differentiation cascade for adaptive behaviour in response to some sort of stressful situation (e.g. healing response to injury): whereas the smallest the value of *e*, the more robust the population, we have also found that upregulation of symmetric self-renewal, i.e. increased value of *p*, increases the likelihood of the population to withstand invasion attempts. However, we found that populations with big values of *p* are not robust, in the sense that lead the population to extinction due fluctuations and the mechanisms studied in this paper.

The issue of the optimal (i.e. *maximally anti-tumour*) architecture of has been recently addressed by Rodríguez-Brenes *et al*. [21]. Comparisons between our results and their's are difficult as the models are based on slightly different premises. For example, the model considered in [21] does not take regulation into account. In turn, our models do not account for self-renewal of the TACs. In spite of these differences, some of the results of Rodríguez-Brenes *et al*. [21] are compatible with ours. For example, Rodríguez-Brenes *et al*. [21] conclude that anti-tumour mechanisms should favour the evolution of shorter (few intermediate transient amplifying stages) differentiation cascades over longer ones. The delay-induced extinction mechanism described in this paper also favours the evolutionary prevalence of shorter cascades. The effect of self-renewal in the intermediate differentiation stages has not been explored in this paper and it is left as a subject for future investigation.

Related models have been proposed by Traulsen *et al*. [7,26,28,29], where they use a Moran process to study the competition between two populations as a function of the duplication probability to find that the probability of invasion is an increasing function of the duplication probability of the invader population. We observe similar results, only for small values of the duplication rate *p*. For larger values of *p*, the probability of invasion decreases, as the invader population becomes extinct due the extinction mechanisms that we have characterized in this paper. This difference is due to the fact that we are considering a model with regulated-SC compartment size. Similarly, the models presented in [16,22–25,27,40] also lack feedback, and, therefore do no exhibit any of phenomeology described in this work.

We also observe significant differences between deterministic and stochastic models. Whereas the mean-field behaviour model presents a critical point, characterized in terms of the net growth rate of the SC compartment, beyond which the resident population stops being robust, the stochastic model exhibits non-monotonic behaviour in which probability: the invasion probability is virtually zero below the deterministic critical point; however, beyond the critical point it reaches a maximum and starts decreasing (owing to the extinction mechanism studied in the symmetric division model).

## Funding statement

The authors gratefully acknowledge the Spanish Ministry for Science and Innovation (MICINN) for funding under grant MTM2011-29342 and Generalitat de Catalunya for funding under grant 2009SGR345.

## Appendix A. Stability analysis of the mean-field models

In this appendix, we analyse the stability of the solution of the mean-field behaviour of the different models we have proposed in this paper.

In §3, we have propose a stochastic model of a differentiation cascade with symmetric SC division. The corresponding mean-field behaviour is described by the following set of ODEs:
A1
A2
A3
A4which has two fixed points, namely, the trivial equilibrium (0, 0, … , 0), which corresponds to extinction of the population, and a positive equilibrium, *P _{S}*, given by
A5Let

*J*be the

*n*×

*n*Jacobian matrix of the vector field described by equations (A 1)–(A 5) at the point

*P*. If all eigenvalues of

_{S}*J*have strictly negative real part, then the solution is linearly stable [41]. In figure 15, we study how stability of

*P*depends on the SC death rate,

_{S}*λ*

_{0}, and self-renewal rates,

*p*and on the death rate of the MCs. We observe that the stability region of the fixed point

*P*expands with as

*λ*

_{0}and

*λ*increase. By contrast, the stability region shrinks when

_{n}*p*increases.

In §4, we have formulated a stochastic model of a differentiation cascade with mixed symmetric and asymmetric SC proliferation. The mean-field behaviour of this system is given by
A6
A7
A8
A9
A10which has two fixed points: the trivial equilibrium (0, 0, … , 0), a positive, *P _{A}*, given by
A11Here, we analyse the effect of the non-asymmetric division probability of the SCs,

*e*, on the stability region of

*P*in the

_{A}*p*−

*λ*

_{0}plane, which is shown in figure 16. In particular, we observe that asymmetric SC division (

*e*< 1) leads to a larger stability region of

*P*.

_{A}- Received March 12, 2014.
- Accepted March 20, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.