Phenotypic heterogeneity in modeling cancer evolution

The unwelcome evolution of malignancy during cancer progression emerges through a selection process in a complex heterogeneous population structure. In the present work, we investigate evolutionary dynamics in a phenotypically heterogeneous population of stem cells (SCs) and their associated progenitors. The fate of a malignant mutation is determined not only by overall stem cell and differentiated cell growth rates but also differentiation and dedifferentiation rates. We investigate the effect of such a complex population structure on the evolution of malignant mutations. We derive exact analytic results for the fixation probability of a mutant arising in each of the subpopulations. The analytic results are in almost perfect agreement with the numerical simulations. Moreover, a condition for evolutionary advantage of a mutant cell versus the wild type population is given in the present study. We also show that microenvironment-induced plasticity in invading mutants leads to more aggressive mutants with higher fixation probability. Our model predicts that decreasing polarity between stem and differentiated cells turnover would raise the survivability of non-plastic mutants; while it would suppress the development of malignancy for plastic mutants. We discuss our model in the context of colorectal/intestinal cancer (at the epithelium). This novel mathematical framework can be applied more generally to a variety of problems concerning selection in heterogeneous populations, in other contexts such as population genetics, and ecology.


Introduction
Cancer can be thought of as a complex ecosystem in which not only tumor cells but also other cell types (phenotypes) may influence the overall health of an organism. Experimental results have recently shown that cancer cells may mimic the functional features of normal cells (Kreso and Dick, 2014). The most important features are associated with a small subpopulation of cells, namely the stem cells. Stem cells (SCs) are defined to be cells with self-renewal capacity and pluripotency. For instance, they can replenish and regenerate the whole epithelial cell population in normal tissues. It has been proposed that cancer stem cells (CSCs) maintain invasive characteristics, such as (undesirable) multipotency and uncontrolled growth and tumor initiating capacity (Borovski et al, 2011;O'Brien et al, 2009;Reya et al, 2001;Sprouffske et al, 2013;Marjanovic et al, 2013). The differentiated progenies of SCs are the cells with specialized distinct functions, within the organism. They are produced via a hierarchical division scheme. As the differentiated cells (DCs) become more mature along the hierarchy, their replication potential decreases (Chaffer et al, 2011a;Legraverend and Jay, 2011;Magee et al, 2012;Ritsma et al, 2014).
CSCs reside in small niches and manifest characteristics similar to somatic SCs (Borovski et al, 2011). In solid cancers, CSCs are usually imputed as a result of the expression of similar biomarkers as those used to identify SCs (Ginestier et al, 2007;Natarajan and FitzGerald, 2007;Singh et al, 2003;Bao et al, 2006;Woodward and Sulman, 2008). In colon cancer, the over-expression of the polycomb ring finger oncogene BMI1 leads to the down-regulation of proteins p16INK4a and p14ARF. These proteins override cellular proliferation restriction and generate cancer SLCs (Dimri et al, 2002;Legraverend and Jay, 2011;Marjanovic et al, 2013). For mammary stem cells, CD44 + and CD24 − are reported as markers for stemness. In acute-myeloid leukemia (AML) CD34 + CD38 − cells are a leukemia-initiating subpopulation Roesch et al, 2010).
Despite the new established dogma that cancer cells originate from a small niche of cells Ritsma et al, 2014;Vermeulen et al, 2013), a range of experiments have now investigated and reported on the cancer initiating capacity of committed progenitor cells (O'Brien et al, 2009;Li and Laterra, 2012;Marjanovic et al, 2013). In other words, non-SCs can undergo a dedifferentiation process and regain stemness (these cells are also called stem like cells). In breast cancer, epithelial-mesenchymal transition (EMT) factors have been implicated in the production of stem like cells from non-stem cells (Mani et al, 2008;Chaffer et al, 2011bChaffer et al, , 2013. Gupta et al (Gupta et al, 2011) have also observed that the epithelial differentiated cells with basal markers can convert to cells with stem cell markers (see also (Chaffer et al, 2011b)). There are several other experimental observations supporting dedifferentiation of committed progenitor cells (Cabrera et al, 2015;Huels and Sansom, 2015;Li and Laterra, 2012;Marjanovic et al, 2013;Philpott and Winton, 2014). In addition, this dedifferentiation has been observed, under certain microenvironmental conditions, in normal SCs (Kreso and Dick, 2014;Tata et al, 2013;Dorantes-Acosta and Pelayo, 2012). In fact, it is becoming apparent that cellular trans-differentiation is activated in a number of organs to produce stem like cells in support of SCs in tissue regeneration (Kreso and Dick, 2014;Li and Laterra, 2012;Marjanovic et al, 2013;Vermeulen et al, 2013;Chaffer et al, 2011a;Fausto, 2000;Fillmore and Kuperwasser, 2008;Legraverend and Jay, 2011).
A variety of quantitative approaches have been utilized to investigate the effect of the stem cell hierarchy and phenotypic heterogeneity on tumor growth (Clayton et al, 2007;Dhawan et al, 2014;Turner et al, 2009;Zapperi and La Porta, 2012;Gao et al, 2013). In the context of cancer evolution, the deterministic population dynamics of the stem cell hierarchy in the absence of plasticity, is discussed by Werner et al (Werner et al, 2011(Werner et al, , 2013. Plasticity and dedifferentiation is explored under a diffusion approximation (Jilkine and Gutenkunst, 2014) and replicator equation (Kamran Kaveh, 2013).In (Shahriyari and Komarova, 2013), the authors discuss the rate of evolution in a simple hierarchical stem and non-stem cell population. They argue that stem cell symmetric division is preferred under natural selection for two-hit mutations.
The evolutionary dynamics of malignant and normal genotypes in the presence of phenotypic transformations (differentiation and dedifferentiation) is not well understood. In this study, we consider a general framework to study natural selection in heterogeneous populations. We analyze competition between resident and mutant populations which are genotypically different. Each of these types divide into phenotypically different subpopulations (stem cell and differentiated subtypes). Due to homeostasis, the size of SC and DC subpopulations are assumed to remain constant. Stem cells can self-renew and replenish their own population or contribute to differentiated cell population via differentiation events. Differentiated cells can also divide into differentiated cells or dedifferentiate into stem like cell states.
We investigate conditions for the successful selection of a malignant mutation in this complex population structure. Due to the plastic nature of the early malignant progenitors, there is a finite chance for an advantageous mutant to exit the differentiated group and become part of the SC niche. We derive analytic results that predict the fixation probability of a mutant (either in the SC or DC subpopulations) to establish a finite colony. We assume arbitrary population sizes and division rates and selection intensities as well as (de)differentiation rates. The analytic results are in excellent agreement with stochastic simulations in finite populations. We apply our findings to colorectal cancer and predict dedifferentiation rates that can confer a selection advantage for p53 mutants.
The paper is organized as follows: in the Material and Methods section, the generalized Moran process is defined and the generating function method to calculate long-term fixation or extinction probabilities is presented. The replicator dynamics for the model is derived in the absence of random drift. In the Results section, we discuss the fixation probability as a function of stem cell self-renewal rate, differentiation rate and dedifferentiation rate. The phase diagrams for advantageous mutants are also presented towards the end of this section. In the Discussion section, we summarize our findings and suggest some possible generalizations to characterize more heterogeneous environments.

Materials and Methods
Generalized Moran process with (a)symmetric division and plasticity Consider two populations of resident or wild type (type 1) and mutant or invader (type 2). Mutants are the result of an oncogenic mutation in the resident population. Each genotype is divided into phenotypically different subpopulations of stem cells (SC) and differentiated cells (DC). Stem cells can self-renew symmetrically where the offspring are stem cells. They can differentiate (symmetrically or asymmetrically) to produce differentiated progenies of the same genotype. We denote the probability of asymmetric differentiation (per division) byû 1 ,û 2 and symmetric differentiation byv 1 ,v 2 . The overall probability of differentiation is u 1,2 =û 1,2 + 2v 1,2 . This is due to the fact that symmetric differentiation produces two differentiated cells. Similarly the self-renewal probability is denoted by 1 − u 1,2 . The indexes 1 or 2 denote the corresponding probabilities for a wild type or mutant. The division rate of a normal (or mutant) stem cell is denoted by r 1 (r 2 ) respectively. Similarly, the division rates of progenitors/differentiated cells are denoted byr 1 andr 2 . For evolutionary Figure 1: Phenotypic-genotypic changes in individuals within a four-compartmental structure. We consider constant population sizes N S and N D for SCs and DCs respectively. With respect to the finite Markov chain, we consider a generalized model to take into account the competition between normal and malignant individuals in each of the SC and DC subpopulations. Differentiation and dedifferentiation events connect the selection dynamics between the two niches. In (a), all possible differentiation, dedifferentiation, and death events with their corresponding rates are represented. The SC-DC compartmental structure is depicted in (b) with the associated self-renewal and differentiation/plasticity possibilities. dynamics we consider a birth-death (BD) Moran process as follows: At each time step, an individual is chosen to reproduce proportional to its fitness within the SC or DC compartments. If a normal (mutant) cell in the SC compartment is chosen to reproduce, its offspring replaces a randomly chosen cell in the stem cell compartment with probability 1 − u 1 (1 − u 2 ). Otherwise, with probability u 1 (u 2 ), the (differentiated) offspring replaces a randomly chosen cell in the DC compartment. Similarly, if a differentiated cell is chosen to reproduce, its offspring replaces another cell in the differentiated cell compartment with the probability 1 − η. Alternatively, the offspring can dedifferentiate into a stemlike cell and replace a randomly chosen individual in the stem cell compartment with a rate η, where η = η 1 (η = η 2 ) denote the dedifferentiation probability for normal (mutant) DCs. For simplicity we assumed death rates of all types to be equal and set this to unity. Net reproduction rate of wild type and mutant stem cells r 1 ,r 2 Net reproduction rate of wild type and mutant differentiated cells u 1 , u 2 Asymmetric differentiation rate of normal and mutant stem cells η 1 , η 2 Dedifferentiation rate of normal and mutant differentiated cells The above dynamics models the differentiation mechanism with an effective asymmetric division with the probabilities u 1 , u 2 . Thus in the following we use the terms differentiation (of stem cells) and asymmetric division interchangeably.
The above Moran process can be written as a continuous time process (1/N is the duration of each time step for N = N S + N D ) 1 N ∂p(n S , n D ; t) ∂t = W + S (n S − 1, n D ) p(n S − 1, n D ; t) + W − S (n S + 1, n D ) p(n S + 1, n D ; t) + W + D (n S , n D − 1) p(n S , n D − 1; t) + W − D (n S , n D + 1) p(n S , n D + 1; t) − W + S (n S , n D ) + W + D (n S , n D ) + W − S (n S , n D ) + W − D (n S , n D ) p(n S , n D ; t). (1) where p(n S , n D ; t) denotes the probability of having n S mutant stem cells and n D mutant differentiated cells, at time t (given n 0 S and n 0 D at t = 0). The population of normal cells are given by N S − n S and N D − n D correspondingly. The probabilities W ± S and W ± D are the transition probabilities corresponding to an increase or decrease by one in the number of mutant SCs and DCs resp. They are given by W + D (n S , n D ) = Prob(n S , n D → n S , n D + 1) The denominator N r denotes the total fitness of SC and DC individuals: The above Markov process has two absorbing states corresponding to fixation or extinction of the mutant or WT. The competition between the two genotypes in the stem cell compartment is tied to the competition inside the differentiated compartment via differentiation and dedifferentiation mechanisms. In the absence of plasticity we have a hierarchical population structure where only mutations in the stem cell compartments can give rise to fixation in the whole population.
In the next section, we present analytic solutions for the probability of fixation, of an invading mutant as a function of division rates and (de)differentiation rates. Our calculations match with simulation results for a wide range of parameters and population sizes.

Fixation probability in a heterogeneous Moran process
One of the most important questions to address within a heterogeneous population is the chance of success for a mutation in different subtypes.
The fixation probability of a mutant originating in the stem cell compartment, ρ S or the differentiated cell compartment, ρ D is a measure of the tumor initiating capacity of each subpopulation. For a completely hierarchical population, only mutants that arise in the stem cell niche have a chance of fixating in the whole population thus ρ S = ρ. If the progenitors can dedifferentiate into stem-like cells, the comparison between the two fixation probabilities, (ρ S and ρ D ), is a good measure of how the tumor initiating capacity correlates with the notion of stemness.
The use of the probability generating function (PGF) method to study a constant population Moran process is discussed in Vallade, 2010, 2011;Kaveh et al, 2015). It is used to present an alternative derivation of the (well-mixed) Moran fixation probability, by identifying a martingale for the process. Here we generalize this technique for a heterogeneous population under selective pressure in the presence of phenotypic plasticity.
A martingale for the above four population model, Eqs.
(1), (3), can be written as where · denotes the stochastic average and the (auxiliary) variables z S and z D satisfy the following system of algebraic equations By matching the initial conditions for t = 0 and the steady state solutions of the PGF, we can obtain analytic expressions for the fixation probability of mutants of each subtype (stem or differentiated). In general, the fixation probability beginning with i mutant SCs and j mutant progenitors is (See SI for details) For i = 1 and j = 0 (starting with one initial SC mutant) the fixation probability is Similarly, the fixation probability of a newborn mutant in the DC compartment Moreover, assuming random mutations, i.e. uniform mutation rates in both compartments, the average fixation probability is given by with N tot = N S + N D . The probability of a successful emergent mutant before time t (from a background of N tot normal cells) is given by where µ denotes the mutation rate.
Stochastic simulation Using the model described above, we performed numerical simulations using such updates until each of the runs tends to saturation in the fraction of SCs and DCs, or until we reach the maximum updating time of T=15,000 for each realization. Then running the whole procedure for 20,000 realizations, we calculated the fraction of results for the fixation probability of SCs and DCs in those runs. Then repeating each calculation for a set of five iterations, we calculated the mean and error bars. Errors are calculated as the standard deviation of the mean.

Replicator equation
The Markov process considered in our four-compartment model exhibits deterministic dynamics in the absence of stochastic fluctuations, i.e. infinite population limit. This replicator equation captures the average frequency of various phenotypes which provides insight into the evolutionary dynamics of the system. Firstly, we detect the frequency of each phenotype at the fixed points and secondly, analyze the phase diagram by varying different parameter values (see the Results Sec.) at equilibrium. Starting from the master equation (1), one obtains the following system of deterministic equations for malignant SCs and DCs (see SI for derivation) where x S = n S (t)/N S and x D = n D /N D are the average frequencies of mutant SCs and DCs resp. Stationary state frequencies, x S and x D satisfy the following coupled system of equations In the next section, after obtaining the exact analytic form of the fixation probability for this generalized Moran process, we will use Eqs. (41) to derive the condition for evolutionary advantage of the mutant genotype and the phase diagram for the model.

Exact analytic results for the fixation probability
We first begin with the analytic expressions for the survival probabilities Eqs. (6)-(9) where z S and z D are the solutions of Eqs. (5). We begin with some simpler limiting cases where compact algebraic results can be obtained, and then proceed to general solutions of Eqs. (6)-(9).
Standard Moran process.. Let us consider the simple case where the differentiation and dedifferentiation rates are set to zero. In this case our model should reduce to two disjoint Moran processes, one for the stem cell and one for the differentiated cell compartments.
For this case, we obtain z S = r 1 /r 2 and z D =r 1 /r 2 and the fixation probability for a mutant to dominate a SC (DC) niche is Another interesting limit that resembles a well-mixed Moran process in the stem cell compartment occurs whenr 1,2 0 but u 1,2 = 0. This is an exaggerated case showcasing the limited proliferation capacity of differentiated progeny. The corresponding average fixation probability for an emerged mutant in the SC or DC compartment is since the fixation probability of a newborn mutant in the DC class will be zero in this case.
Invasion in hierarchical model (zero plasticity).. We now consider a more general case with no plastic potential, η 1,2 = 0. For the moment, we assume r 1 =r 1 = 1, r 2 =r 2 = r. Fig. 2 shows how the survival chance of an initial mutant SC varies in terms of the relative fitness of mutant SCs/DCs as well as the probability of asymmetric division in the SC compartment. The population size is set to N S = N D = 10. The results are plotted for three sets of differentiation rates (u 1 = u 2 = 0.5), (u 1 = 0.1, u 2 = 0.5) and (u 1 = 0.5, u 2 = 0.1) as r varies ( Fig. 2(a)). If the normal cell differentiation rate u 1 is decreased from the balanced limit of u 1 = u 2 = 1/2 the fixation probability decreases as well. Conversely, if the mutant cell differentiation rate, u 2 is decreased, away from u 2 = 1/2, ρ S would increase. We rewrite differentiation rates as u 1 = u and u 2 = u where u stands for an overall measure of differentiation, in both genotypes, and measures the asymmetry in them. If < 1 the normal type differentiates more often per division and > 1 indicates that invader cells differentiate more often. Fig. 2(b) shows ρ S for = 1/2 as a function of u for several division rates, r = 1, 1.2, 1.5, 2. For high differentiation rates, u, the fixation probability converges to unity. This is consistent with Eq. (15). The relative fitness of the mutant versus normal types are given by r(1 − u 2 )/(1 − u 1 ) = r(2 − u)/(2(1 − u)). As u → 1 the relative fitness approaches infinity and thus the value of fixation probability tend to unity.  (c) Figure 2: The fixation probability of mutants in the absence of plasticity. We assume N S = N D = 10, r 1 =r 1 = 1, and η 1 = η 2 = 0 in simulations (points) and analytic calculations (solid lines). Each error bar shown at each point is the standard deviation of mean. In (a) changing parameters u 1 , u 2 , that are the differentiation rates of normal and tumor SCs resp., the trends for the fixation probability of mutants is given as a function of relative fitness of mutants, referred to as r 2 =r 2 = r. In (b) and (c) the fixation probability variation is given in terms of asymmetric differentiation rates u 1 = u 2 = u and various values of r and the ratio of the differentiation rates of normal SCs = u 2 /u 1 . In (b) = 0.5 and in (c) = 1.5.
Invasion in the presence of plasticity (dedifferentiation).. Now we consider a more general case with non-zero differentiation and dedifferentiation rates. As before we set r 1 =r 1 = 1. We parameterize the fitness of mutant subtypes as r 2 = αr andr 2 = βr. (α is different from the one introduced in the previous figure.) The parameter r denotes the overall proliferation advantage of mutants over normal cells. Fig. 3 shows ρ as a function of r for various values of α and β. Fig. 3(a) assumes plastic mutants only (η 1 = 0, η = 0.5) whereas Fig. 3(b) assumes equally plastic genotypes (η 1 = η 2 ). The differentiation rates are equal and set to 1/2 for simplicity. As can be seen, the overall increase in r increases the fixation probability monotonically in both cases. However, varying values of α  Relative Fitness, r Fixation Probability, ρ S Analytic, α=0.5, β=1 Analytic, α=1, β=1 Analytic, α=2, β=1 Moran, α=0.5, β=1 Moran, α=1, β=1 Moran, α=2, β=1 (b) Figure 3: The fixation probability of mutants in the presence of phenotypic plasticity. We suppose that N S = N D = 10, r 1 =r 1 = 1, r 2 = αr, andr 2 = βr in the analytic results (solid curves) and stochastic simulation (points with bars as the standard deviation of mean). Changing α from 0.5 to 1 and then to 2 while β remains fixed, the behavior of the system is shown in terms of the fixation probability with respect to the relative fitness r. In subfigure (a) plastic potential has only been considered for malignant individuals: η 1 = 0, η 2 = 0.5 while in (b) both WT and mutant cells can dedifferentiate to the stemness state (η 1 = η 2 = 0.5). Straightforward calculations reveal that for the given parameter values in (b), ρ = ρ S = ρ D . and β leads to the neutral value for the proliferation potential r. For example, in the case of Fig. (a) for α = β = 1, i.e. r 2 =r 2 = r, the mutant is advantaged for r > 0.75. For α = 0.5, β = 1, however, r > 1.25 has ρ > 1/N . Interestingly if stem cells divide faster than progenitors, i.e. α > β, even for r ≈ 1/2 the fixation probability is larger than the neutral value. We used the neutral fixation probability as 1/N S = 1/N D = 1/10. Which is the neutral ρ in the absence of (de)differentiation. Similar observations can be made in the case of η 1 = η 2 in part (b). As can be seen our results are in very good agreement with exact stochastic simulations.
In Fig. 3, ρ S is plotted as a function of r for several values of u 1 and u 2 (part (a)) as well as η (part (b)). The effect of dedifferentiation on increasing the value of the fixation probability is most significant for r values corresponding to the neutral limit ( Fig. 3(b)) and Fig. 2(a). For example in Fig. 3(b) value of ρ near r = 1 increases from approximately 0.1(≈ 1/N ) to 0.3 as η increases to 0.5 (from 0.01). However for r = 4 the difference between the two values of ρ for η = 0.01, 0.5 is negligibly small.
In the presence of dedifferentiation, mutations arising in the differentiated compartment can now undergo clonal expansion. In Fig. 5 we plotted values of the fixation probabilities ρ S , ρ D and the average fixation ρ as a function of η. As expected for η = 0, ρ D tends to zero for various values of u 1 and u 2 . However ρ S approaches the results obtained before (Eq. (15)). In Fig.5(a) for u 1 = u 2 = 0.5 and r = 1.1 ρ S approaches 0.15 as η = 0. This is in agreement with the extended Moran result, Eq. (15), with effective fitness r(1 − u 2 )/(1 − u 1 ) which for N S = 10 gives ρ S = 0.15. As η increases ρ S increases further and approaches 0.4 for large values of η. We can see that finite values of η have now conferred  Figure 4: Effect of change in asymmetric differentiation and plasticity rates on survivability of mutants. We assume that N S = N D = 10, r 1 =r 1 = 1, and r 2 =r 2 = r. In subfigure (a), the fixation probability of SCs as a function of η is given, where η = 0.01, 0.1, 0.5 while u 1 = u 2 = 0.5, and η 1 = 0. In (b) η 1 = 0 and η = 0.1. Changing parameters u 1 and u 2 , which are the asymmetric division rates of normal and tumor SCs resp., the fixation probability as a function of u 1 , u 2 is shown. Solid lines represent the analytic calculation and points correspond to simulation results (error bars are based on the standard deviation of mean). a selection advantage on the previously deleterious mutants. For example for u 1 = 0.5, u 2 = 0.1 and r = 1.1, the effective fitness is less than unity with ρ S ≈ 0. For finite values of the dedifferentiation rate, ρ S will exceed the neutral value 1/N . For example for η = 0.5 we can read ρ S ≈ 0.25 in this case.

Phase diagram
To derive the phase diagram that determines the evolutionary advantage of a malignant genotype, we begin with the system of equations (40). The fixed points of the replicator equation determine the stationary frequencies of mutant subtypes: x S , x D . There are four fixed-points for the system of equations (40). The fixed point (x S = 0, x D = 0) corresponds to the extinction of both phenotypes and other fixed points (x S = 1, x D = 1) or any (x S = 0, x D = 0) corresponds to successful invasion by the mutant. We do not distinguish between these fixed points as the stochastic fluctuations prohibit coexistence of mutant and WT phenotypes. We determine the phase boundary when a fixed point (x S , x D ) merges with (0, 0). This indicates that the invasion has become unstable to parameter changes, and determines the phase boundary in the phase space of parameters.
Similar results can be obtained from the Eqs. (6)-(8) requiring that ρ S = 1/N S (corresponding to the neutral selection). In a large population, this approach would be identical to the results from the replicator equation.
Next, we consider the case with r 1 =r 1 = 1, r 2 =r 2 = r, u 1 = u 2 = u, η 1 = 0, and η 2 = η. We have plotted the phase diagram for the model parameter , and (c) we have resp. considered the fixation probability of mutants while the initial malignant mutation resp. occurs in SC, DC, and SC+DC compartments at random. In these plots, we assumed that N S = N D = 10, r 1 =r 1 = 1, r 2 =r 2 = 1.1, u 1 = 0.5 and η 1 = 0. Different asymmetric division rates of normal and malignant individuals have also taken into account when the fixation probability is given as a function of η. Solid curves and points are respectively the results of analytic and simulation analysis. At each given discrete point, the error bar depicts the standard deviation of mean.
space of η − u − r. The phase boundary that determines the condition for neutrality is given by the algebraic relation (see SI for details) (u + η − 1)r 2 + −u 2 + (η − 1)u + 2 − η r − 1 + u 2 = 0, The phase diagram is plotted in the space of u − r ( Fig. 6(a)) and η − r ( Fig.6(b)). The observation that an increase in η can turn a previously deleterious mutant into a beneficial one, can be seen here as well. For example in Fig.6(b) for u = 0.7 and r = 0.9 < 1 where we expect a deleterious mutant for large enough values of η, we can cross the phase boundary into an advantaged region. Interestingly we observe a similar trend as the overall differentiation rate increases.

Application to colorectal cancer
Clonal expansion in colorectal cancer is known to be initiated as a result of mutations occurring at the bottom of the crypt (Vogelstein et al, 2013). More Phase diagram of plastic mutant SCs. The phase boundary for advantageous and disadvantageous mutant populations are given as differentiation and plasticity rates change. We assume that r 1 =r 1 = 1, r 2 =r 2 = r, u 1 = u 2 = u, and η 1 = 0. Different regions for advantageous and disadvantageous mutant SCs are given in (a) as u changes. A similar analysis has been carried out in (b) as η varies. In (a) η = 0.1, 0.3, 0.7, here the alteration in the plasticity rate of DCs results in a tendency to approach various regions of fixation for mutant SCs, while the extinction domain shrinks with increasing η. In (b) u = 0.1, 0.3, 0.7. Increasing the asymmetric division rate u, the region for advantageous mutants expands to provide a higher survival chance for mutant SCs. In both cases, advantageous criteria relate to either fixation of mutants or coexistence of mutants and WT individuals.
recently, Vermeulen et al (Vermeulen et al, 2013) considered the dynamics of cells at the bottom of a normal colonic/intestinal crypt. Such controlled cells in their circular model of size 5 undergo a selection process. The selection mechanism is investigated for several oncogenes and tumor suppressor genes imposed at the crypt base. Fixation probability of a single mutation and the relative fitness r of the mutant cells, compared with the normal host, are reported. Their estimates reveal that the relative fitness of the original cell containing APC −/− is r = 1.58, while r = 1.56 for Kras G12D , and r = 0.96 for P53 R172H (compared with normal control cells mice). P53 mutation seems not to confer a fitness advantage and is weakly deleterious. However for P53 R172H , it has been observed that the fitness of a mutant elevates from 0.96 to 1.16 in comparison with the DSS-treated cells (colitis) which is in the presence of inflammatory injury. Thus, under inflammatory signaling effects, the P53 R172H mutants appear to gain a selection advantage and thus a higher fitness (Vermeulen et al, 2013(Vermeulen et al, , 2010. According to the recent in vivo study by Schwitalla et al (Schwitalla et al, 2013), inflammatory signaling plays a role in elevating the rate of dedifferentiation. It has been also shown that inflammatory disease activates the transcription factor NF-κB. NF-κB can, in turn, elevate Wnt-signaling which leads to the phenotypic plasticity of non-SCs (Karin and Greten, 2005).
Thus we suggest that the higher survival chance of the P53 mutated SCs along with the DSS-treated cells, may be the result of dedifferentiation in the presence of inflammatory stroma (Schwitalla et al, 2013). The survival probability of mutants in colon/intestine, in the presence of inflammatory signaling, is presumably correlated with both the fitness and plastic nature of the epithelial cells. Thus the fixation probability reported in (Vermeulen et al, 2013) of cells with P53 DSS colitis mutation, can be derived by the same fitness r = 0.96 for P53 when dedifferentiation occurs.
For instance, when r = 0.96 and u 1 = 0.5, u 2 = 0.25, having the plasticity rate at η = 0.12, one obtains the same fixation probability as (Vermeulen et al, 2013) (see Table 2). This finding strongly suggests that there is an elevation in the fixation of a deleterious mutant into an advantageous trait due to the plastic properties of the mutant (see also the next section). In Fig. 7, we considered a cylindrical model of the crypt-base, and immediate adjoining transit amplifier layers, cells in the colon/intestine. The bottom layer consists of central stem cells and border stem cells, while transit ampflifying and non-SC cells reside on top of this layer and at higher layers. A mouse crypt is comprised of 5-7 functional stem cells (Potten et al, 1992;Vermeulen et al, 2013). To compare the results of our model with experimental observations of (Vermeulen et al, 2013) for the fixation probability of P53 DSS colitis, we consider the first two circular layers of SCs (see Fig. 7) in which we roughly assume N S ≈ 10 (Ritsma et al, 2014) and note that the fitness reported in the experimental data (Vermeulen et al, 2013) is for functional SCs and may be considered unchanged for the two layers (N S = 10) of central (functional) and border SCs as well.

Discussion
The evolutionary implications of epigenetic heterogeneity are not very well understood in cancer biology. A known picture for phenotypic heterogeneity (when the genotypes are assumed to be identical) relates to the cancer stem cell hierarchy. In this picture, pluripotent cells with tumor initiating capacity can undergo mitotic events and either replenish their own population or produce a lineage of partially differentiated cells (including precursor and/or transit amplifying cells).
In the current study, we present a general model of four distinct subpopulations to investigate Darwinian evolution in such a hierarchical structure. We consider two genotypically different populations (mutant and wild-type). Mutations are results of unwanted oncogenic or TSG mutations. Each genotype has phenotypically different subtypes of stem cells (SCs) and differentiated cells (DCs). SCs ultimately generate their associated progenies, DCs, through proliferation and asymmetric differentiation. DCs have restricted proliferation capacity. Due to the tissue structure of the crypt, the population of different subtypes remains approximately unchanged. Genetic mutations can occur among SCs or DCs which we assume to be occurring through a uniform probability. Mutations can confer not only higher division rates but also different rates of differentiation and plasticity among mutant subtypes. These changes can be triggered, for example, by microenvironmental conditions.
Our model predicts the fixation probability of a newborn mutant as a function of the division rates of mutant and resident SCs/DCs, differentiation and dedifferentiation rates. Exact analytic calculation and numerical simulationswhich are in almost perfect agreement -suggest that the asymmetric differentiation in the SC group has a major effect on the fate of mutants compared with dedifferentiation. More specifically, a greater impact on the fixation probability of SCs can be observed by the change in asymmetric differentiation of normal cells compared with that of malignant cells. Furthermore, we observe that the more plastic trait has an evolutionary advantage. This is most notable close to the neutral limit, i.e. when the proliferation rates of the mutants and residents are very close in value. Most interestingly, a sufficient increase in the rate of plasticity can turn a previously deleterious mutant into a beneficial one.
As an important application of this model, we consider the intestinal/colonic crypt with two groups of SCs and their neighboring, partially differentiated cells. The competition between malignant mutations and normal cells in the base of the crypt has attracted much research interest, and is one the most studied scenarios in cancer evolution. It has recently been shown that the Moran type process, initially suggested by (Kaveh et al, 2016(Kaveh et al, , 2015Komarova, 2006Komarova, , 2007Nowak, 2006;Nowak et al, 2002) are in perfect agreement with the experimental observation (Vermeulen et al, 2010). In this in vivo experimental analysis, KRAS, APC +/− , APC −/− , and P53 DSS colitis mutations separately induced in the crypt base. Assessing clonal lineage tracing (Vermeulen et al, 2010), the authors were able to observe the growth of mutants in the populations and thus measured the fixation probability of mutants. However, the experimental observations sometimes ignore the microenvironmental interaction between the neighboring transit amplifying (TA) cells and SCs in the crypt (See (Ritsma et al, 2014;Vermeulen et al, 2010) for the models and estimated parameter values of population dynamics of the crypt, as well as numbers of central and border SCs at the base of the crypt).
In this investigation, we have suggested a plausible estimation method for the dedifferentiation rate which can give the same value of the fixation probabilities of P53 inactivation in the crypt base as observed by (Vermeulen et al, 2013). As the authors in (Vermeulen et al, 2013) noticed, P53 mutations can maintain advantageous features to invade the crypt base in the presence of inflammatory signaling. We propose that this could potentially be a result of dedifferentiation caused by an inflammatory stroma (Schwitalla et al, 2013). This finding supports our idea about the elevation of a deleterious mutant into one with an advantageous trait due to the plastic properties of the mutant.
This research is one of many mathematical approaches recently used to investigate various aspects of tumorigenesis. In our approach, we attempt to dissect a part of the complex machinery in multi-compartmental models, and understand this in greater detail. This study provides various insights concerning some features of initiation and progression of cancer, and suggests possible experimental investigations to confirm some of the theoretical results. The approach taken in this paper may lead to a better understanding of the natural pathological mechanisms in the colonic/intestinal epithelium or any similar four-compartmental structures in ecology, population genetics, and social networks. The particular application of this approach to carcinogenesis may have some results at prognosis.

A. Characteristic equation
Starting from the master equation of the given model (where for N = N S + N D , we assume that 1/N is the duration of each updating time) the generating equation can be derived by assuming the probability generating function (PGF) in which coefficients define the probabilities of being at different possible states after a given time.
Starting from the master Eq. (17) discussed in the previous section, we can analyze the probability of absorption (fixation or extinction) using generating function techniques recently developed for the constant population birth-death processes (Houchmandzadeh and Vallade, 2011). Defining the probability generating function (PGF) of the probability density function p(n S , n D ; t) of having n S mutant SCs and n D mutant DCs at time t. Then the PGF is Using the PGF one can rewrite the master equation (the latter equation) for PGF to derive the characteristic equations as follows where the operators W ± S and W ± D are W + S (n S , n D ) = Prob(n S , n D → n S + 1, n D ) W ± D (n S , n D ) and W ± D (n S , n D ) define the transition probabilities for increase/decrease in the number of mutant SCs and DCs resp. Acquiring these probabilities with those for no change in the number of mutant cells in either compartments, one can define the transition matrix which incorporates random walks among various states prior to fixation.
Substituting Eqs. (20) into generating equation (19), we obtain: For large -N S (N D ), we can simplify the latter equation for the probability generating function by keeping the linear derivative terms to leading order in N S (N D ). This tends to where the operator can be considered constant when N S and N D are set to be large, N S , N D 1. Setting the coefficients of the derivatives ∂F ∂zS and ∂F ∂zD to zero, we obtain approximate quasi-stationary points of this constant population dynamics at the large-t limit. This relates to the corresponding martingales for N S , N D → ∞ branching process limit. Denoting the solutions with z S and z P , one attains

B. Fixation probability
Taking advantage of the generating function, an exact analytic approach for the fixation probability can be derived even in the presence of plasticity (when mutation and mutation-back do not occur). The results are obtained for a BD Moran process; however, a similar calculation can be performed for a Voter (DB) Model with presumably different fixed points for the same initial conditions. The boundary and initial conditions for the corresponding generating function are resp. as follows where i and j are the initial number of cancer SCs and DCs resp. For the following special cases, the system has two absorbing states at equilibrium which signify fixation and extinction for mutant cells in either compartments. We denote the probability of reaching extinction and fixation states by B 0 and B 1 resp. Thus we obtain the following result from the PGF at steady state based on the boundary condition (26), B 0 + B 1 = 1. Biologically, since the mutant DCs are produced by cancer SCs, extinction of cancer SCs will result in replenishment of mutant DCs. Moreover, the co-operation between SCs and DCs suggests that it is impossible to have mutant SCs fixate while DCs become completely extinct. We conclude that Finally, applying the initial condition, F (z S , z D , t = 0) = z i S z j D and the boundary condition, F (z S = 1, z D = 1, t) = 1, the fixation probability for an initial population of i malignant stem cells and j progenitors is derived as where (z S , z D ) is the nontrivial fixed point of the generating Eq. (25).
C. Finding the quasi-fixed points and the associated survival probabilities An interesting scenario occurs when the normal component is not plastic, i.e. η 1 = 0, in which case the above equations can be simplified to the following expression of z S , z D = 1 z D = r 1 (1 − u 2 ) (r 1 + r 1 u 1 ) z S + r 1 (1 − u 1 ) (r 1 + r 1 u 1 ) Eqs. (5) can be reduced to the following closed form where The latter equations together with the original system suggest the following solution for z S and thus represent another fixed point of the problem: Now let us consider the following particular cases which give rise to some interesting consequences of the model.

C1. Standard Moran process
Assume that there is no transition (migration) between SC and DC compartments (two islands), that is, u i , η j 1 for i, j = 1, 2. Then each compartment will follow the mass-action BD Moran model with the fixed points for z S and z D where the overall behavior is akin to two disjoint Moran processes. The solutions to Eqs. (25) are (1) z S = 1, z D = 1, (2) z S = 1, z D =r 1 r2 , (3) z S = r1 r2 , z D = 1, (4) z S = r1 r2 , z D =r 1 r2 . where the solution (4) accounts for the fixed points of two separated Moran models in SC and Dc groups. Therefore, the fixation probabilities of starting from one imposed mutant in SC and DC groups resp. are Another interesting limit that results in a well-mixed Moran model in the SC class occurs whenr 1,2 0 where SCs are also committed to reach a stationary state. Then the solutions to Eqs. (25) are given by (1) z S = 1, z D = 1, (2) z S = r1 r2 , z D = 1.
Compared with the result from the previous case, the second solution relates to the solution of a well-mixed model occurring in the SC compartment where the only absorbing state for mutant DCs is to take over the whole population where the evolutionary scenario accounts for cancer DCs domination. In this scenario, the average fixation probability is C2. Invasion in hierarchical model (no plasticity) Now we consider the case where no plastic potential is taken into account during the whole process, we assume that r 1 =r 1 = 1, r 2 =r 2 = r, u 1 = u 2 = u, and η 1 = η 2 = 0. in this situation, one obtains the following solutions (1) z S = 1, z D = 1, where u eff = r(1−u2) 1−u1 is the effective asymmetric division rate and Γ 2 = 1 − 6ru 2 u 1 + u 2 1 + 2ru 2 − 2r − 2u 2 − 2u 1 + 2ru 1 + r 2 + 2u 2 u 2 1 + 2u 2 2 u 1 + u 2 2 + u 2 2 u 2 1 .
The third solution cannot be maintained in reality since there is no cooperation between the compartments, which could result in support for the mutant SC by the mutant DC group. In such a case, the only possible outcome would be the fixation of the mutant DCs. Collectively, we find C3. Invasion with phenotypic plasticity (dedifferentiation) Suppose that dedifferentiation only occurs for the mutant DCs at a rate η 2 = η for η η 1 ≈ 0. Also let assume that r 1 =r 1 = 1, r 2 =r 2 = r, u 1 = u 2 = u and η 2 = η, then the solutions for z S and z D are (1) z S = 1, z D = 1, (2) z S satisfies in the equation The values for z S and z D , in this case, satisfy the following relations (for more details refer to Supplementary Information): where A = r(1 − u), B = 1 − u, C = ru, E = r(1 − η), F = u + 1, G = rη. The derived solutions introduce the possible fixed points of the characteristic equations. We investigate this case in more detail later and through analyzing the phase diagram of the generalized model for non-zero plasticity in the replicator dynamics.

D. Phase diagram of the system
Now, acquiring the replicator dynamics of the four compartment model which depict the alterations in the the average frequencies of mutant SCs and DCs, resp. x S = n S (t)/N S and x D = n D /N D , one obtains the following system of equations: At equilibrium, the fraction of mutant stem and non-stem cell groups would lead to a specific states which are the pseudo-fixed points of the problem. These type of fixed points can be attractive or repulsive, depending on the initial conditions of the system. As we described in the paper, according to the cooperation between mutant SCs and DCs (and similarly between normal SCs and DCs), the malignant individuals may become extinct or survive together. Now defining the fixed point as (x S , x D ), it may tend either to (0, 0) or to (1, 1). Such a criteria suggests two distinct phases for the fate of the malignant cells which are separated by the phase boundary. At steady state, the following system can be derived for the fixed points x S and x D : To analyze the latter system, we may consider two important cases: Case 1.. At first, we assume u 1 = u 2 = η 1 = η 2 = 0, then the solutions to the system (41) results in the extinction or fixation of mutant SCs and DCs. In this case, the phase diagram is simply divided to two advantageous (r > 1) and disadvantageous (r < 1) cases.
Among the given solutions (1)-(4) of the present case, the only acceptable non-trivial solution is (4). The solution (3) does not satisfy the condition 0 ≤ x D ≤ 1 for all possible values of r, u, and η. Moreover, as we described above, the cooperation among mutant cells results in the fixation of both mutant SC and DC groups to the state (1, 1).
The last solution, then, implies that having (x S , x D ) → (0, 0) leads to the following solution A = (u + η − 1)r 2 + −u 2 + (η − 1)u + 2 − η r − 1 + u 2 = 0, which characterizes the advantageous and disadvantageous regions for mutants (see Fig. (6) for more details). Another limit relates to the case where (x S , x D ) approaches (1, 1). One obtains r η = (r − 1)(u − 1) This condition, which is not defined for r > 1, does not change the phase diagram and would not have any effect on the selection pressure of the system on mutants.