Effect of Dedifferentiation on Time to Mutation Acquisition in Stem Cell-Driven Cancers

Accumulating evidence suggests that many tumors have a hierarchical organization, with the bulk of the tumor composed of relatively differentiated short-lived progenitor cells that are maintained by a small population of undifferentiated long-lived cancer stem cells. It is unclear, however, whether cancer stem cells originate from normal stem cells or from dedifferentiated progenitor cells. To address this, we mathematically modeled the effect of dedifferentiation on carcinogenesis. We considered a hybrid stochastic-deterministic model of mutation accumulation in both stem cells and progenitors, including dedifferentiation of progenitor cells to a stem cell-like state. We performed exact computer simulations of the emergence of tumor subpopulations with two mutations, and we derived semi-analytical estimates for the waiting time distribution to fixation. Our results suggest that dedifferentiation may play an important role in carcinogenesis, depending on how stem cell homeostasis is maintained. If the stem cell population size is held strictly constant (due to all divisions being asymmetric), we found that dedifferentiation acts like a positive selective force in the stem cell population and thus speeds carcinogenesis. If the stem cell population size is allowed to vary stochastically with density-dependent reproduction rates (allowing both symmetric and asymmetric divisions), we found that dedifferentiation beyond a critical threshold leads to exponential growth of the stem cell population. Thus, dedifferentiation may play a crucial role, the common modeling assumption of constant stem cell population size may not be adequate, and further progress in understanding carcinogenesis demands a more detailed mechanistic understanding of stem cell homeostasis.


Introduction
Most tissues consist of three classes of cells: stem cells, transitamplifying progenitor cells, and differentiated cells. Multicellular organisms require a tight control of cell division to ensure a proper balance between these different cell populations. The cancer stem cell (CSC) hypothesis states that tumors are also hierarchically organized, with a small sub-population of cancer cells driving cancer growth [1]. Individual cell tracing studies of tumor development strongly support the cancer stem cell hypothesis in many (but not all) types of cancer [2,3], and identifying these cells in tissues is an ongoing goal in cancer research. Lineage studies find that malignant tumors contain more cancer stem cells compared to benign tumors and that cancers gradually lose their tissue-like hierarchical organization as they evolve from benign to malignant [2].
Cells escape proliferation control after acquiring a series of mutations in a multi-step process [4]. While some cancers may require only a few mutations [5], the number of required (driver) mutations in solid cancers is larger, with up to twenty driver mutations being required [6]. In order to accumulate this critical number of mutations during a lifetime, cells either have to be longlived or the mutation rate has to be large [7]. Stem cells have been proposed to be likely candidates for the initial cell of mutation due to their long lifetime and sustained self-renewal capacity [1]. In addition to their long life span, stem cells are able to generate full lineages of differentiated cells, thereby perpetuating mutations through clonal expansion. Given known division and mutation rates, theoretical studies have argued that the necessary number of mutations for carcinogenesis cannot be obtained in the stem cell population on a reasonable time scale without assuming either significant selective advantage or elevated mutation rates [4,7]. However, there is conflicting evidence as to how early in tumor development cancers acquire an elevated mutation rate [8,9] and several cancer genome sequencing studies have estimated mutation rates during cancer initiation to be normal for some types of cancer [10][11][12].
Although a stem cell may sustain the first oncogenic hit, subsequent alterations required for development of CSCs can occur in descendent progenitor cells [13]. Dysregulation of pathways involved in stem cell self-renewal may lead to progenitor cells acquiring a stem cell-like phenotype. It remains an open question whether cancer stem cells originate from stem cells that escape homeostasis or from dedifferentiated progenitor cells that acquire infinite proliferating potential [14]. There is significant evidence that dedifferentiation can play a role in establishment of certain cancers [15][16][17]. For example, cell sorting has demonstrated that stem-like cells can arise de novo from non-stem-like cancer cells in in vitro breast cancer cell lines [18,19]. In the hematopoietic system, it has been shown that leukemic stem cells can be generated from committed progenitor cells that acquire stem cell-like behavior [20]. It has been suggested that acute myeloid leukemia (AML) is a progenitor disease, where a progenitor acquires abnormal self-renewal potential and ''dedifferentiates'' to a stem cell-like state [21,22]. Other myeloid leukemias such as CML (chronic myeloid leukemia) are thought of as stem-cell diseases [23]. However, although a hematopoietic stem cell is thought to be the cell of origin in the early phases of CML, in patients with CML blast crisis, granulocyte-macrophage progenitors are thought to acquire self-renewal capacity through a b-catenin mutation and emerge as the probable CSCs [24]. Using mathematical modeling to investigate the likelihood of mutation occurring in a progenitor versus a stem cell is a continuing line of investigation [25]. We treat the probability of a mutant progenitor cell acquiring stem cell-like state as a ''dedifferentiation'' rate, and we study how this parameter influences the time to carcinogenesis. We are primarily interested in whether dedifferentiation can speed up the time to tumor development in hierarchically organized cancers and in what rates of dedifferentiation are necessary for a noticeable effect.

Prior Related Mathematical Modeling
Certain aspects of the cancer stem cell hypothesis have previously been addressed by mathematical models. It has been shown that having a hierarchical tissue design, where a small population of stem cells maintains a transient population of differentiating cells, may slow the accumulation of mutations and protect against cancer [26][27][28]. The question of whether genetic instability (resulting in hyperactive mutation rate) is an early or later event in mutation acquisition leading to cancer has been addressed by several groups (see [4] for review). Most mathematical models find that the onset of genetic instability should be an early event, if at least some of the mutations are neutral. However, sequencing suggests that the mutator phenotype is expressed relatively late in cancer progression [9].
Stem cell populations are typically small. Hence, the dynamics of mutant cells in the stem cell population are highly sensitive to stochastic fluctuations. A tumor begins with a single mutated cell, so there is a substantial chance of mutant extinction due to random events. Genetic drift and stochastic clonal extinction in stem cell lineages have been experimentally demonstrated for both normal tissue stem cells [29][30][31] and cancer stem cells [2] in several tissue types. Consequently, a deterministic model of mutation acquisition in stem cells will significantly underestimate the time to cancer establishment [32]. Many models of mutation acquisition use a stochastic approach and are concerned with calculating time to emergence or fixation (or when the number of mutant cells reaches some threshold value used in diagnosis) of a mutant cell with fitness r~1zs in a population of size N sc .
The waiting time for cancer is often defined as the time until a particular number of mutation events have occurred in at least one cell. Iwasa et al. [33] considered a two-stage Moran model and described conditions under which ''stochastic tunneling'' can occur. (In this phenomenon, cells with two mutations reach fixation before cells with one mutation reach fixation.) Durrett et al. [34] obtained asymptotic estimates of waiting times until a cell with i mutations first appears under the assumption of neutrality (s~0). These models typically consider a fixed population size [5,23,25,[35][36][37][38][39]. The fixed population assumption is supposed to reflect homeostasis in the stem cell population, though how homeostasis is achieved is typically not addressed. Although the Moran model captures the stochastic nature of mutation acquisition, this type of model is not capable of describing mutations that change the stem cell division pattern and result in possible expansion of the stem cell pool, which in turn leads to tumor growth. Some recent models also consider mutation accumulation in exponentially growing cell populations [40][41][42][43]. Beerenwinkel et al. [6] used the Wright-Fisher model with exponentially growing population size to look at the effect of selection on the waiting time to cancer, and they predicted that the observed genetic diversity of colorectal cancer genomes can arise under a normal mutation rate (taken to be u~10 {7 per cell division) if the average selective advantage per mutation is on the order of 1%. Similar calculations using a discrete branching process found s~0:4% given u~10 {5 [40]. Note that increased mutation rates due to genetic instability would allow even smaller selective advantages during tumorigenesis, but neutral mutants (s~0) result in waiting times that are too long compared with disease incidence. Other groups have also concluded that for normal mutation rates and neutral mutants, mutations in multiple genes in acquired hematopoietic disorders are most likely very rare events, as acquisition of multiple mutations typically requires development times that are too long compared to disease incidence [36].
Spencer et al. [44] and Ashkenazi et al. [45] have focused on the sequential order of mutations associated with increased rate of proliferation, decreased rate of death, increased mutation rate, and other hallmarks of cancer that must accumulate before emergence of cancer. The sequence of mutations with the shortest waiting time to getting all the necessary mutations is considered the most likely mutational pathway [25,44]. However, these models do not consider the possibility that dedifferentiation of progenitor cells can affect the time to multiple mutation acquisition.
The dividing progenitor cell population has previously been described by multi-compartment ODE models, with cells moving between compartments as they age [45][46][47]. Note that in these models the exact number of different stages of differentiation is ambiguous and does not exactly correspond to mitotic events, as cells may undergo more than one division in each compartment stage [46]. Most of these models of age-structured cell populations assume a stem cell proliferation rate that is dependent on the total number of cells and thus incorporate negative feedback as a means of achieving homeostasis [48,49]. These deterministic models have

Author Summary
Recent evidence suggests that, like many normal tissues, many cancers are maintained by a small population of immortal stem cells that divide indefinitely to produce many differentiated cells. Cancer stem cells may come directly from mutation of normal stem cells, but this route demands high mutation rates, because there are few normal stem cells. There are, however, many differentiated cells, and mutations can cause such cells to ''dedifferentiate'' into a stem-like state. We used mathematical modeling to study the effects of dedifferentiation on the time to cancer onset. We found that the effect of dedifferentiation depends critically on how stem cell numbers are controlled by the body. If homeostasis is very tight (due to all divisions being asymmetric), then dedifferentiation has little effect, but if homeostatic control is looser (allowing both symmetric and asymmetric divisions), then dedifferentiation can dramatically hasten cancer onset and lead to exponential growth of the cancer stem cell population. Our results suggest that dedifferentiation may be a very important factor in cancer and that more study of dedifferentiation and stem cell control is necessary to understand and prevent cancer onset.
focused on mechanisms that could regulate cell numbers that are necessary for homeostasis and efficient repopulation. We use a similar mathematical approach to model the progenitor population as [49], but we couple it to stochastic dynamics in the stem cell compartment.
Upon division a stem cell can produce zero, one, or two stem cells with probabilities a D , a A , and a S , respectively (Fig. 1A). The mean number of stem cell offspring is given by a A z2a S . If symmetric divisions are permitted, the stem cell population can be described by a branching process with the expected number of cells at time t given by (a A z2a S ) t . However, a branching process either goes extinct or undergoes exponential growth, and thus it cannot capture stem cell dynamics at equilibrium. One solution is to use a conditional branching process [50], where the probabilities for a branching process are conditioned to the total population size remaining constant by an unspecified sampling mechanism (i.e., assuming that the stem cell population remains in homeostasis). Some theoretical studies have previously considered the impact of the asymmetry of cell division on stem cell dynamics. However, these stochastic models all assumed a fixed stem cell population size, either through a variant of the Moran process [35,51] or conditional branching process [39]. We utilize a different approach to get a time-varying but bounded stem cell population size in our models.

Our Modeling Approach
We use mathematical modeling to study how the possibility of ''dedifferentiation'' of mutant progenitor cells into a stem cell-like state affects the waiting time to carcinogenesis. Dividing progenitor cells have large growing populations, so we use a deterministic model to describe their evolutionary dynamics. For stem cell populations, stochastic effects are important, because the proliferating stem cell population is typically small. We use a stochastic model for stem cell dynamics as a boundary condition to the PDE governing differentiated cell expansion ( Fig. 1B and C.) There is also feedback from the deterministic progenitor population to the stochastic stem cell population as a rate of ''dedifferentiation''.
To assess the effect of dedifferentiation on time to carcinogenesis, we consider models for stem cell dynamics with both fixed and variable stem cell numbers (Fig. 1D)

Models
We make the following assumptions: 1. The population of progenitor cells is large enough that maturation can be treated as a continuous variable. Discretizing the progenitor cell population based on the number of divisions a cell has completed, we obtain an age-structured partial-differential-equation model for the number of differentiated cells of age a at time t. We assume that when progenitor transit-amplifying cells carrying i mutations divide, they produce progenitor cells of the same maturity stage, obtaining the linear PDE where p i (a,t) is the progenitor cell density at age a and time t, s(a) is the age-dependent proliferation rate, and m(a) is the age-dependent mortality rate. We assume that the rate of maturation da dt does not depend on age a and, without loss of generality, set it equal to 1. Similar age-structured population equations have been previously studied, with focus on the regulatory feedback mechanisms that are necessary for homeostasis and structural stability of the steady state solution [47,49]. 2. We assume that there is a separate stem cell population that gives rise to newly born differentiated cells p i (0,t) that serves as a boundary conditions to the PDE system in Eq. (1). We treat this population stochastically. We consider two different models for the stem cell population: one where the total stem cell population is fixed and a variant that allows a time-varying but bounded stem cell population size (Fig. 1D). 3. We assume neutral fitness of mutant stem cells, with the proliferation advantage of the mutant phenotype appearing only in the progenitor stage, in line with what is known for some cancers of the hematopoietic system [52,53]. 4. We require M~2 mutations to appear in the progenitor population before dedifferentiation to a stem-cell like state is possible, because sequencing of acute myeloid leukemia genomes suggest that there are two driver mutations present [54]. Additional justification for requiring M~2 mutations is considered in the Discussion. 5. Due to lack of data on dedifferentiation capacities, we assume progenitor cells of all maturity stages have an equal probability of dedifferentiating.

Progenitor Cells
Extending Eq. (1) to account for mutations between multiple subpopulations of progenitor cells (Fig. 1C) we obtain Here u Ã is the mutation rate per cell per unit time and p i (a,t) is the number of progenitor cells of ''age'' a from the subpopulation with i mutations. We assume 0ƒiƒM, and no back mutation is allowed. Let n i (t) be the number of stem cells with i mutations at time t. Let a D,i be the probability of a symmetric division that gives rise to two differentiated cells, a A,i be the probability of an asymmetric division that gives rise to one stem cell and one differentiated cell, and a S,i be the probability of a symmetric division that gives rise to two stem cells. Then If we neglect mutation, the steady wave-form solutions of Eq. 2 have the form where a~2a D,i za A,i is the average number of stem cells of type i produced per division and r i (a)~Ð a 0 (s i (s){m i (s))ds is the agedependent growth rate of the differentiated cell population (Text S1). Hence, the long-term age distribution is largely determined by the functional forms of the differentiated cell birth and death rates ( Fig. S1 and S2.) Altered birth and death rates due to mutations can result in mutant subpopulations growing to higher plateaus in size, but the final population size will be bounded. Our PDE system can be easily modified to have a maximal carrying capacity K i for each sub-population. This does not qualitatively change the age distribution of progenitor cells (Fig. S1) and does not significantly affect the fraction of i-mutation cells in the total progenitor population (Fig. S3), so we do not consider it further.
To mimic a maturity switch for cellular proliferation and death, we took the proliferation and death rates of differentiated cells per unit time to be Here b and d are the maximal proliferation and death/removal rates of progenitor cells. The age at which the proliferation switch occurs (i.e., half the progenitor cells stop dividing) is given by v b , and the steepness of the proliferation switch is determined by r b . Similarly, the age at which half the cells begin to undergo apoptosis is given by v b , and the steepness of the death switch is controlled by r d . If v b vv d , then differentiated cells between the ages of v b and v d are not replicating (senescent). Note that setting either of these values to zero results in a uniform rate of birth/ death. Effects of varying proliferation/death parameters are shown in Fig. S2. The parameters governing proliferation, in particular b and v b , have much larger influence on the final differentiated cell population size than parameters governing death/removal. The steepness of the switch does not substantially change the age distribution.

Stem Cells
Constant stem cell population size. To model the evolutionary dynamics of a stem cell population under strict homeostasis (resulting in a fixed stem cell population size), we used the Moran stochastic process for Mz1 different types, with mutations between types and neutral fitness [50]. Let the number of individuals carrying each possible number of mutations be given by n~(n 0 ,n 1 ,:::,n M ), where X M i~0 n i~Nsc . We considered two versions of this model, with and without dedifferentiation. In both cases, we assumed that each stem cell divides, on average, every T gen chronological time units. Thus, in a population of size N sc , the average time between divisions was T gen =N sc .
In the first model, no dedifferentiation of progenitor cells was possible. Every T gen =N sc time units, a single randomly chosen stem cell j was removed and one cell i was born with probability given by where m i,h is the probability of changing to type i from type h per replication event, and e j is a unit vector with 1 in the jth column. We considered a linear cascade of mutations in which the mutation matrix ½m i,h is given by We also considered a version of the model in which dedifferentiation of two-mutation differentiated cells was allowed, but the total stem cell population size remained fixed. In this model, the probability of death of a j-mutation stem cell and birth of an imutation cell was given by where e is the proportion of cells in the stem cell pool that come from dedifferentiated cells at each replication event, and d i,2 is the Kronecker delta function signifying that that only two-mutation progenitor cells can dedifferentiate. Here is the proportion of two-mutation cells of all ages in the progenitor population, given that p i (a) is the density of differentiated cells of age a carrying i mutations. We also considered a version of the model in which all progenitor cells, regardless of the number of mutations, could dedifferentiate (Text S1). Because the Moran model has been studied extensively, we were able to use several existing results on the time to emergence and fixation of mutants. Let t M be the first time at which an individual carrying M mutations emerges who will go on to fix in the population. We focus on the case M~2 because sequencing of acute myeloid leukemia genomes suggests that there are 2 driver mutations present [54]. (See Discussion for more details.) Using branching process approximations, Durrett et al. [34] calculated the waiting time for the Moran model under neutral drift of prior mutants. For M~2, the probability density function for t 2 is given by This simplifies to where l~N sc u 1 , and r fix is the probability that a single mutant individual will fix in a population of size N sc . For neutral drift, r fix~1 =N sc , and for weak selection where M-mutation cells have advantage s%1 [50]. The time to fixation of the subpopulation with M~2 mutations is a sum of two random variables: the time t 2 until appearance of a successful two-mutation cell (Eq. (9)) and the waiting time t fix from the time that mutant first appears until that mutant fixes [55]. Note that this time is given in units of stem cell generation times T gen . The probability density function of the total fixation time, T fix , is given by the convolution of the probability density functions w for time to first appearance of successful mutant andw w for time it takes that mutant to fix.w w can be obtained from the backward Kolmogorov equation for the probability of fixation f of a gene with initial frequency p 0 before time t: subject to boundary conditions f (1,t)~1 and f (0,t)~0 and initial condition f (p 0 ,0)~d(p 0 ). Dividing f by the ultimate probability of fixation r fix and differentiating with respect to t, we obtain the probability density function forw w as a function of initial allele frequency p 0 [56].
Variable stem cell population size. Our previous stem cell models couple birth and death events to keep the population size fixed, but we next decoupled these events to allow for a stochastically varying population size. For clarity, we refer to the total stem cell population size in this model by S(t). Again assuming that the average replication time of a stem cell is T gen , the interval between birth/death events in this stochastic stem cell model corresponds to T gen =S(t) time units in the progenitor cell model. We assume that homeostasis in the stem cell pool is maintained by control of cell fate upon division, and that each stem cell can produce zero, one, or two stem cell offspring. For example, the possible offspring from a zero-mutation stem cell are: two differentiated zero-mutation cells with probability a D , one zero-mutation differentiated cell and one zero-mutation stem cell with probability a A (1{u), one zero-mutation differentiated cell and one one-mutation stem cell with probability a A u, two zeromutation stem cells with probability a S (1{2u), and a one zeromutation stem cell and one one-mutation stem cell with probability 2ua S . (For simplicity we assume the probability of both offspring carrying new mutations to be negligible.) In general, a stem cell carrying i mutations can produce k stem cell offspring carrying j mutations with probability P j i (k) given by With constant division probabilities a D ,a S and a A , this model is a rescaled Galton-Watson branching process, and the stem cell population either goes extinct in finite time or undergoes exponential growth when the mean number of stem cell progeny per cell, m~2a S za A , is greater than one [57]. To describe a stem cell population under homeostasis, the probabilities a D , a A , and a S must depend on the total stem cell population size S(t). We model a carrying capacity K i for stem cells carrying i mutations, such that As g approaches one, a A approaches one, so that most of the divisions that occur do not change the stem cell population size, and we recover the Moran process for all carrying capacity values. The parameter j controls strength of fluctuations about the carrying capacities; as j increases, the fluctuations become smaller. In this model, newly emerging mutants can still go extinct due to stochasticity, but the total population reaches a quasi-stationary regime at population size K i . Note that, because this is a quasistationary regime, eventually the stem cell population will go extinct, but the expected time to extinction is exponentially proportional to K i for j~1 and g~0 [58]. We chose the carrying capacity to be large enough that extinction of the stem cell population does not occur on a physiological timescale, and we initialized the stem cell population to be at carrying capacity with zero-mutation cells.
We considered two versions of the variable stem cell population size model. In the first case, no dedifferentiation was possible. In the second case, differentiated cells with K mutations were allowed to dedifferentiate and re-enter the stem cell population. Let d be the dedifferentiation rate per two-mutation progenitor cell per unit time. Then the mean number of cells dedifferentiating in the interval T gen =S(t) between two stem cell replications is To introduce dedifferentiated cells into the stem cell population, at each replication event we calculate the mean number of dedifferentiated progenitor cells l and update the stem cell population: ½n 0 ,n 1 ,n 2 ?½n 0 ,n 1 ,n 2 zPoisson(l), ð17Þ After the total stem cell number is updated, the probabilities of reproduction are re-calculated using Eq. (14), and reproduction is carried out.

Model Parameters
Parameters used are summarized in Table 1. We used parameter estimates from the human hematopoietic system because parameters for other cancers are less well known. We used M~2 as the number of necessary mutations to develop a cancerous phenotype. Although it has been estimated that for the human hematopoietic system there are 11,000-22,000 stem cells [59], which give rise to all blood and immune system cells, most of these cells are quiescent and only divide when body sustains an injury and needs to repopulate the hematopoietic system. Our model only considers actively dividing stem cells, which have been estimated by various methods to number around 100 [32,60]. The entire actively dividing stem cell population has previously been modeled as turning over once per year [32], but most recent estimates have an individual stem cell dividing every 25-50 weeks [61]. However, this is likely an overestimate, as it is difficult to distinguish between actively dividing and quiescent stem cell populations. We assume that an active stem cell divides every 20 weeks, which when multiplied by N sc results in an active stem cell population turnover time of T gen~5 weeks. (The entire stem cell population including quiescent cells turns over on a much longer timescale.) Whereas the size of the active hematopoietic stem cell pool is small, the number of progenitor cells such as granulocyte, erythroid, monocyte, and megakaryocyte colony-forming units (CFU-GEMM) and granulocyte and monocyte colony-forming units (CFU-GM) is much larger. There are approximately 10 5 CFU-GEMM cells and 10 8 CFU-GM cells [62]. There are estimates that each CFU-GEMM may contribute to hematopoiesis for an average of 60 days (range of 40-340 days) and that it replicates at an average rate of once every 50 days (range of 35-285 days) [62]. We track the progenitor populations for L~20 weeks, and assume that their proliferative potential rapidly drops off after 10 weeks. The maximal proliferation and death rates, b i and d i were chosen so that 100 stem cells results in 10 5 {10 6 progenitor cells of all ages.
Not much is known about the selective advantage s provided by driver mutations for different cancer types, except that it is small (r~1zs&1). Unless stated otherwise, we assume neutral fitness in the stem cell pool (s~0) in our stochastic models throughout the paper, to focus on the effect of dedifferentiation. We use a range of s~0 Á Á Á 0:4 for the progenitor cells in the deterministic model.
Mutation estimates per cell division per gene range from about 10 {7 in normal cells to 10 {2 in the case of chromosomal instability [63]. (Note that the rate of epigenetic change has been estimated to be orders of magnitude higher than that of genetic change and could also play a role in cancer initiation [10].) A common value used in many mathematical models is a driver mutation rate of u~10 {5 per division, obtained by assuming a somatic mutation rate of 10 {7 per gene, and about 100 genes that could be mutated to give same phenotype [40,45]. In normal hematopoietic cells the mutation rate has been measured as u~10 {6 per division [64].
Note that in the stochastic model, which considers every cell division, the mutation rate u can be used as is, but using chronological time (i.e., weeks or months) means that this value should be multiplied by the average number of divisions per unit time to obtain u Ã . (Mutations that speed up the cell cycle will then speed up the apparent mutation rate per unit of chronological time in our progenitor model.) The expected number of doublings from n i stem to p i progenitor cells is log 2 p i =n i ð Þz2 [46], and the total number of progenitors cells of type i is p i &an i e Ð r(a), da . Using values from Table 1, this results in 8{10 cell divisions that take place over 10 weeks, so u Ã &u in equations (2).

Results
The coupled system of stem cells and progenitor cells undergoing mutation and dedifferentiation we modeled is complex. To disentangle the effects of different phenomena, we systematically built up the model. We first considered the progenitor population alone. We then considered the stem cell population alone, in models with strict and variable stem cell homeostasis. Finally, we coupled the stem and progenitor populations through dedifferentiation.

Progenitor Population Alone
We first considered whether mutation and reproduction in the progenitor population could by itself generate a sustained population of two-mutation cancerous cells. We thus modeled a scenario in which no stem cell mutations occur, so the boundary condition to the progenitor population system in equations (2) is simply (N sc ,0, Á Á Á ,0). Because selection in the progenitor population might favor mutants, we also assumed that progenitor cells with i mutations have a proliferation rate b i~( 1zs) i b 0 (Eq. (5)). This yields a steady-state age distribution of normal and mutant progenitor cells (Fig. S2). Fig. 2 summarizes results for typical parameter values, showing that for M~2 mutant cells to be an appreciable fraction of the population, the mutation rate u and proliferative advantage s must both be unreasonably high. This is true both if the total progenitor population can grow without bound (Fig. 2a) and if its growth is restricted (Fig. 2b). Similar findings are obtained if competition between progenitor subpopulations is included in the model (Fig.  S3). Consistent with previous work [26,27,36], these results show stem cell dynamics cannot be ignored in considering time to carcinogenesis, so we next considered stochastic models of the stem cell population.

Stem Cell Population Alone
Our first models for stem cells did not incorporate dedifferentiation, so the dynamics were entirely governed by the stem cells.
In modeling cancer, the time to carcinogenesis can be defined as the time for a single M-mutation cell to emerge, the time for Mmutation cells to pass some threshold number or fraction, or the time for M-mutation cells to fix in the population. If the mutation rate is low (such that N sc u%1), then all three definitions are similar, because the time to emergence of a successful M-mutation cell is long compared to the time from emergence to fixation. However, there is large uncertainty regarding effective mutation rates in carcinogenesis (Table 1), so the assumption of low mutation rate may not always be valid, and we thus calculated times to fixation.
We began our stem cell modeling by considering fixed population size, corresponding to strict homeostasis. In this constant N sc case, we could leverage several analytic results, with which our simulations agreed well. Fig. 3A shows a typical simulation. The full probability density distribution of time to fixation is given by Eq. (12) and agrees well with our simulations for high mutation rates (Fig. 3B). The time to emergence of a successful mutant is of order 1/( ffiffiffiffiffiffiffi N sc p u 3=2 ) stem-cell generations (Eq. (9)). For normal mutation rates of u~10 {6 {10 {8 per cell division, the mean time until emergence of a two-mutation cell is 10 8 {10 11 stem cell generations, which is very long even with a short stem cell generation time.
Because homeostasis is likely imperfect, we also considered a stochastically fluctuating stem cell population size. We found that, without dedifferentiation, the distributions of times until fixation are very similar for models with and without fluctuations in the stem cell population size, as long as we condition on nonextinction of the stem cell population (Fig. 3B). This is true for a wide range of probabilities of asymmetric division g and strengths of mean reversion j (Eq. (15)). This agrees with previous findings that demographic stochasticity does not alter fixation times of neutral mutants in a large population [65], provided that the carrying capacities of the mutants are the same.
Our results suggest that dynamics within either the progenitor or stem cell compartments considered separately do not result in carcinogenesis in the hematopoietic system on a realistic timescale, provided that cancer-causing mutations occur at normal mutation rates, selection advantages relative to wild-type stem cells do not appear until M~2 mutations, and the stem cell population size is constant or varies stochastically around a carrying capacity. We thus turned our attention to coupled model systems in which progenitor cells can dedifferentiate into stem cells.

Dedifferentiation with Constant Stem Cell Population Size
For the coupled system, we first considered stem cell homeostasis caused by strict asymmetric division in the stem cell population, so the stem cell population size remains fixed. To model dedifferentiation in this case, we built off the Moran model and assumed that when a stem cell dies and another enters the population, the new entrant comes from the two-mutation progenitor population with probability equal to e times the proportion of two-mutation cells in the progenitor population. Otherwise the new stem cell comes from replication of another stem cell. Roughly speaking, in this model the death of a stem cell leaves a opening in the niche, which can potentially be filled by a dedifferentiated progenitor cell. The number of progenitor cells which can successfully dedifferentiate is controlled by the number of niche openings (stem cell deaths), not by the absolute number of progenitor cells.
Typical simulation results are shown in Fig. 4A. We found that dedifferentiation dramatically shortens the time to fixation of twomutation cells (Fig. 4B). For small dedifferentiation rates e 0:05, we also saw good agreement between our simulations and a semianalytical approximation for the time to fixation of two-mutation cells with selective advantage e (Eq. (12)). This agreement suggests that under strict stem cell homeostasis, dedifferentiation is effectively equivalent to a growth advantage for mutant stem cells.
Distributions of times to fixation of two-mutation stem cells are plotted as a function of both dedifferentiation rate e and mutation rate u in Fig. 4C. Dedifferentiation had two major effects in this model: increasing the probability that an emergent two-mutation stem cell would fix and reducing the time between emergence and fixation. Both of these effects act only after a two-mutation cell has been generated in the stem cell population. (Recall that, as shown in Fig. 2, the mutation rate and selective advantage must be unrealistically high for a nontrivial fraction of two-mutation progenitor cells to exist in the absence of underlying two-mutation stem cells.) For all mutation rates u, the distribution of times to fixation was roughly constant for dedifferentiation rates e 1=N sc , consistent with population genetics theory that selection is only effective when the selection coefficient is greater than the reciprocal of the effective population size. For small mutation rates u, increasing e beyond this threshold only marginally shortened the total time to fixation. This is because in this case the total time to fixation is dominated by the time for a successful two-mutation cell to emerge, and dedifferentiation only reduces this time by a factor of 1= ffiffiffiffiffiffiffi r fix p (Eq. (9)), where r fix is the probability of a emergent two-mutation stem cell fixing. Under and sigmoidal birth rate with maximal growth rate b 0~2 , b i~( 1zs)b i{1 for i~1,2. In (B) the carrying capacity used is N 1~2 00N sc ,N 2~2 50N sc ,N 3~3 00N sc : Other parameters are as in Table 1. For two-mutation cells to reach appreciable levels in this scenario, both the mutation rate and the proliferative advantage must be unreasonably large. doi:10.1371/journal.pcbi.1003481.g002 neutrality r fix~1 =N sc , so for our model with N sc~1 00, dedifferentiation can shorten the time to emergence by at most a factor of 10. The dedifferentiation rate needed to significantly change this waiting time scales linearly with N sc (Fig. S4C). Hence, for larger stem cell population sizes, a small dedifferentiation rate would have a larger effect. For high mutation rates u, the effect of dedifferentiation is more dramatic, because the time from emergence to fixation of two-mutation cells, which dedifferentiation also shortens, is comparable to the time to emergence (Fig. 4D).
The model considered in Fig. 4 assumes that only two-mutation progenitor cells can dedifferentiate. We also considered an alternate model in which any progenitor cell can dedifferentiate (Text S1). In this alternate model, dedifferentiation again had little effect for e v * 1=N sc . Past that threshold the effect was substantial, because in this model dedifferentiation speeds up the time to emergence of two-mutation cells, because one-mutation cells fix much more quickly when they too can dedifferentiate (Fig. S4D). In addition, we considered the case in which the dedifferentiation rate is additionally weighted by the progenitor proliferation rate, and our results did not change qualitatively (Text S1, Fig. S4B).
Our analytical and numerical results suggest that, with intact homeostasis in the stem cell population and normal mutation rates, dedifferentiation plays a fairly minor role in speeding up the time to cancer initiation. We thus turned to consider the case in which homeostasis is not strict.

Dedifferentiation with Variable Stem Cell Population Size
In the previous section, we assumed that the stem cell population size was constant because homeostasis was maintained by all divisions being strictly asymmetric. Consequently, dedifferentiated progenitor cells could only occupy newly created openings in the stem-cell niche created by a death event in the stem cell population. Because homeostasis is likely maintained at the population level [66], with each stem cell division producing not strictly one stem cell but rather on average one stem cell, we next considered a model in which the stem cell population could stochastically fluctuate around a carrying capacity. In this model, stem cell homeostasis was maintained by dynamically altering the probabilities of the three possible outcomes of a stem cell division: two stem cells, one stem and one progenitor cell, or two progenitor cells (Eq. (15)). Two-mutation progenitor cells each had a probability per unit time of dedifferentiating, and dedifferentiated cells were simply added to the stem cell pool. Thus in this model the total influx of dedifferentiated cells depended on the total number of two-mutation progenitor cells, not on the creation of openings in the stem cell niche. (Note that, in our previous model with constant stem cell population size, the rate of dedifferentiation per reproduction event was denoted e. To distinguish the present model, we denoted the progenitor dedifferentiation rate per cell per unit time as d.) Again, we asked whether dedifferentiation substantially speeds the time to carcinogenesis. Fig. 5A and 5B show typical results from this model for a moderate dedifferentiation rate d. After a waiting time, the population of stems cells began to grow exponentially, because the influx of dedifferentiated two-mutation progenitor cells exceeded the capacity of stem-cell division homeostasis. For larger dedifferentiation rates, the exponential growth rate is larger  Fig. 5C and 5D), and the distribution of progenitor ages can be distorted, with many young cells, as seen in Fig. 5E and 5F. Exponential growth eventually occurs whenever the dedifferentiation rate exceeds a threshold d crit . Solving self-consistently for the influx of dedifferentiated cells and the growth rates of the stem and progenitor cell populations, we obtained an integral equation for the growth k k~ad which provides an excellent fit to the numerical simulations ( Fig. 5  and 6A,B). (For derivation details, see Text S1.) Setting this growth rate k to zero, we found Here g is probability of asymmetric stem cell division (producing one stem and one progenitor cell), and T gen is the mean time between stem cell divisions. (Note that if g~1, this model reduces to the Moran model with the population size monotonically increasing due to dedifferentiation.) Lastly, a~2a D,2 za A,2 in Eq. (19) is the average number of progenitor offspring produced by a two-mutation stem cell. Because a D,2 changes as the system attempts to maintain stem-cell homeostasis, a is actually a stochastic variable that depends on the stem cell population size. During exponential growth a&2{g, because the probability of symmetric divisions that give rise to two stem cells goes to zero, and all new stem cell growth comes from dedifferentiated progenitor cells. In Eq. (19), r(a) is the growth rate of twomutation progenitor cells as a function of age a, so Ð ? 0 e r(a) da is the number of progenitors produced by one two-mutation stem cell. Increasing the amplification of mutant stem cells into progenitors increases the net dedifferentiation rate, lowering the threshold d crit . Because the threshold d crit depends on the age distribution of the two-mutation cells, for a given (small) rate of dedifferentiation d, evolving a mutant that proliferates faster (increasing e r(a) ) can destabilize a system in which the number of cancerous cells is stable and take it into exponential growth regime.
The dependence of the critical dedifferentiation rate d crit on the growth-rate advantage s of two-mutation progenitor cells and probability g of asymmetric cell division is shown in Fig. 6B. The critical d decreases rapidly as the selective advantage of twomutation cells increases. Increasing g or T gen also lowers the critical dedifferentiation rate, because homeostasis is less effective when asymmetric stem cell divisions are less frequent. Note that the exponential growth rate k does not depend on the mutation rate (Fig. S5A), and although the critical d given by (19) needed for exponential growth is a function of the probability of asymmetric division g, the actual growth rate k and the time to exponential growth are not significantly affected by changing g (see Fig. S5B).
For dedifferentiation rates d below d crit , two-mutation stem cells eventually fix in the population, but for dwd crit , the stem cell population is likely to begin exponential growth before fixation of two-mutation stem cells. Thus in Fig. 6C and 6D we report the time to carcinogenesis as the time for the two-mutation stem cell population to exceed N sc , the nominal carrying capacity of the stem cell compartment. In this case of stochastic stem cell homeostasis, dedifferentiation can dramatically shorten the time to carcinogenesis, even for low mutation rates u. This is because the first two-mutation stem cell often arises not from direct mutation of a stem cell, but rather from dedifferentiation of a progenitor cell generated by mutations within the progenitor compartment (Fig. 6E). Although mutations in the progenitor compartment do not affect a large fraction of progenitors, because the number of progenitor cells is so large, the absolute number of two-mutation progenitor cells is non-negligible. Thus even small rates of dedifferentiation can have dramatic effects. This is in contrast to the case of strict stem cell homeostasis, in which the absolute number of two-mutation progenitor cells was unimportant, because they needed an opening in the stem cell niche to successfully dedifferentiate.
Our results show that the case of stochastically controlled stem cell homeostasis is qualitatively different from the case of strict homeostasis. If homeostasis is controlled at the population level (where stem cell decisions between symmetric and asymmetric division are stochastic), dedifferentiation can overwhelm it, leading to exponential growth of the stem cell population. Moreover, if dedifferentiated cells do not depend on openings to colonize the stem cell niche, dedifferentiation can dramatically hasten the time to carcinogenesis, even for low mutation rates.

Discussion
Progression to cancer is associated with expansion of the cancer stem cell (CSC) population, but the origin of these CSCs remains unclear. Although CSCs may arise directly from adult stem cells, they may also arise from somewhat differentiated cells that have dedifferentiated and acquired stem cell-like characteristics [13,14,18,19,67]. Stems cells replicate indefinitely, giving them a long time to accumulate the mutations that drive carcinogenesis, but the population of actively dividing stems cells (N sc ) is small. Progenitor cells replicate only a small number of times, but the population of progenitor cells is typically several orders of magnitude larger than the stem cell population. Thus, as a population, progenitors undergo many more divisions, potentially letting some of these cells acquire mutations that enable them to dedifferentiate and drive carcinogenesis. Here, using mathematical modeling, we have shown that even a small rate of dedifferentiation may drastically shorten the time to cancer emergence, even for low mutation rates.
Recent studies suggest stem cell dynamics during homeostasis are governed by neutral competition and genetic drift [10,29,30]. Traditionally, stem cells were thought to always undergo asymmetric division, always yielding a stem cell and a progenitor cell, resulting in a fixed stem cell population size. This scenario is represented by our first model for stem cell dynamics, based on the popular Moran model. It has been recently shown, however, that symmetric divisions also occur in adult stem cells and may be the predominant form of division [68,69]. Moreover, cancer stem cells have been shown to undergo more symmetric divisions than normal stem cells [70]. Little is known, however, about how the stem cell population size is regulated [29]. Hence, in our second model for stem cell dynamics, we made the simplifying assumption of an a priori carrying capacity K i . We considered a densitydependent stochastic process, in which the degree of mean reversion is controlled through the probabilities of producing zero, one, or two stem cell offspring. In this model, the non-constant stem cell population size S(t) tends to return to the carrying capacity K i , because the mean number of stem cells produced per division is greater than one when S(t)vK i and less than one when S(t)wK i . (Although the stem cell population size could, in principle, be maintained by regulating apoptosis rather than biasing division, previous modeling suggests that regulating division probabilities rather than cell cycle time or removal is more important for maintaining homeostasis [46,71].) If the stem cell population size varies and is regulated by biasing division, we found two distinct regimes. If the dedifferentiation rate is much less than a critical value, then the initial two-mutation stem cell often arises from a normal stem cell, so the time to fixation of such a cell is similar to the case with constant population size. If the dedifferentiation rate exceeds the critical value, however, then the initial two-mutation stem cell often arises from a dedifferentiated progenitor cell, so the time to fixation is dramatically shorter than the case with constant population size. Moreover, in this regime the stem cell population eventually grows exponentially, as dedifferentiating progenitor cells overwhelm stem cell homeostasis. Note that the threshold between these two regimes is independent of the overall mutation rate, if stem and progenitor cell mutation rates are proportional.
When the stem cell population size is constant, dedifferentiation simply acts like a selective advantage for mutant stem cells. When the stem cell population size is allowed to vary, however, dedifferentiation can additionally drive exponential growth of the stem cell population. If the stem cell population size N sc is constant, our results imply that stem cell dynamics in the coupled stem cell-progenitor system can be approximated by a population genetics model of the stem cells alone, as long as that model includes positive selection. In this case, we found that the dedifferentiation rate e must exceed 1=N sc to substantially shorten the time to cancer acquisition, similar to classical population genetics results that the selection coefficient must exceed the inverse population size to be effective. Hence, our model predicts that in tissues where the niche contains fewer cells, smaller rates of dedifferentiation are sufficient to influence the time to cancer. It is interesting to note that cancers where dedifferentiation has been shown to occur have a small niche size (i.e. intestinal crypts) [15,16]. For the hematopoietic system, based on the available literature, we assumed that the number of actively dividing stem cells is N sc~1 00, so the dedifferentiation rate must be &0:01 or higher to significantly shorten the time to cancer.
Here we focus on the hematopoietic system, in which the stem cell compartment consists of N sc &100 active cells, and two mutations are necessary for carcinogenesis. For some other cancers, such as colon cancer, the number of stem cells per compartment is much smaller, there are many compartments, and the number of necessary mutations is larger. For high mutation rate, the mean time to fixation scales linearly with N sc (see Fig.  S4C). So in cancers with small N sc two-mutation stem cells will fix much faster. However, the need to accumulate more mutations will slow carcinogenesis. We expect, however, that the qualitative effects of dedifferentiation will be similar to the hematopoietic system we analyzed.
The fact that mutants take a long time to reach an appreciable fraction of the stem cell population is not typically considered in the cancer modeling literature, which often makes an implicit assumption that a newly emerged mutant cell will not go extinct and will fix quickly. Our results show that, for high mutation rate, the time for a mutation to fix in the population is comparable to time for a successful mutant to first emerge, in accordance with classical results of Kimura and Ohta [56]. This is especially important if division events are rare and the population size is large. Considering the time to some predetermined diagnosis threshold is similar to considering the time to fixation, because the time between a selected mutation becoming common and fixing is typically short [72]. Hence, elevated mutation rate (genetic instability) may not speed up time to carcinogenesis as much as is typically assumed, suggesting that some form of selection (potentially through dedifferentiation) is necessary. Most tumors accumulate hundreds of mutations, but the number of necessary ''driver'' mutations depends on the type of cancer. We considered M~2 mutations, because sequencing of acute myeloid leukemia Figure 6. Fixation and exponential growth of two-mutation cells with dedifferentiation for variable stem cell population size. A: Observed growth rate k of the stem cell population (black curve) and the semi-analytic approximation Eq. (18) (green) for g~0, j~1, and u~0:01. The vertical line denotes d crit . B: Analytically predicted critical dedifferentiation rate d crit as a function of asymmetric division probability g and the growth advantage s of the two-mutation progenitor population. Exponential growth occurs for dwd crit . C: Normalized histogram (stars) of waiting times for exponential growth of the stem cell population with stochastic homeostasis and dedifferentiation for u~0:01, d~0:01. For comparison the histogram (red and black dots) as well as the analytical distributions of times to fixation given strict homeostasis for e~0:01 and e~0 are also shown. D: The median and inter-quantile range of times to first occurrence of N sc~1 00 two-mutation stem cells, given stochastic homeostasis and a range of dedifferentiation rates d. For comparison, the waiting times to fixation for N sc~1 00 given strict homeostasis (shaded areas) for the equivalent value of e are also shown. E: The probability that the first two-mutation stem cell arose from mutation in the stem cell compartment, rather than dedifferentiation. Vertical line denotes d crit . Parameters for all simulations given in Table 1 genomes suggest that there are two driver mutations present [54]. Moreover, recent findings on induced pluripotent stem cells also suggest M~2, as loss of both copies of the tumor suppressor protein p53 [73] or the activation of two oncogenes [67] may be necessary for dedifferentiation. Disabling both copies of p53 improves the efficiency of reprogramming to a stem-like state and greatly enhances the production of induced pluripotent stem cells [74,75]. The loss of p53 also leads to the emergence of tumor cells bearing functional and molecular similarities to stem cells [22,73]. Finally, inactivation of p53 changes the ratio of symmetric to asymmetric division in mammary stem cells, allowing the total stem cell population to escape homeostasis [70].
Our model only considers actively dividing stem cells, which in the human hematopoietic system have been estimated to be roughly 100 [32] out of 11,000-22,000 total stem cells [59]. A more complete model would consider both the active and quiescent stem cell populations. Transitions between these states may be influenced by the progenitor population size, potentially acting as a negative feedback and regulating the proliferation of cancer stem cells. In our models, cancerous cells take over the stem cell population, but the ratio of cancer progenitor cells to cancer stem cells is fixed by the progenitor growth process. Even when dedifferentiation drives exponential growth of the stem cells, it is their absolute number that increases, not their proportion in the population. This is in concordance with some in vitro studies, which suggest a fixed proportion of CSCs in a tumor [18].
Many theoretical models find that in order to accumulate multiple mutations on a reasonable time scale, the onset of elevated mutation rate (i.e., genetic instability) should be an early event in tumorigenesis (reviewed in [4,7]). The importance of genetic instability, however, depends on assumptions about symmetric self-renewal and differentiation of stem and progenitor cells. In particular, mutations that alter stem cell division or make committed progenitors somewhat immortal may also lead to an early onset of cancer, diminishing the impact of genetic instability [45]. Similarly, our results show that different assumptions about how dedifferentiation occurs (frequency-dependent reproduction versus absolute numbers of dedifferentiating cells) dramatically alter time to carcinogenesis.
A large body of modeling work in this area (reviewed in section Prior Mathematical Models) has focused on calculating the time to carcinogenesis under the assumption of constant population size (not specifying the mechanism of homeostatic regulation). We compared the times to multiple mutation acquisition in our constant and variable stem cell population size models and found that without dedifferentiation both models yield similar results. With dedifferentiation, however, we found that the two models differ substantially. We explicitly considered different ways that homeostasis can be maintained in the stem cell population, and showed that these assumptions can lead to very different results. Our results suggest that if homeostasis is controlled through division asymmetry and if de-differentiated cells do not depend on openings to colonize the stem cell niche, then for de-differentiation rate larger than a critical threshold, the cancer stem cell will most likely originate in a progenitor cell that has undergone dedifferentiation. This is a prediction of our model that can be experimentally tested using inducible genetic labeling, the same technique that permitted lineage-tracing experiments allowing quantification of symmetric versus asymmetric divisions [29,30]. A similar method was previously used to identify oligodendrocyte precursor cells as the tumor cell of origin in glioma [76].
Our model contributes to the existing literature on the trade-offs between symmetric and asymmetric divisions of the stem cell population in stem cell-driven [35,51,77]. To our knowledge, our model is the first to quantify the effects of dedifferentiation on the time to carcinogenesis. There are a number of important aspects of homeostasis our model does not consider, such as spatial aspects of stem cell position in the niche (we assume the cells are wellmixed) or negative feedback to stem cell divisions from the progenitor population (we assume progenitors only influence the stem cells by dedifferentiating). Lander et al. have previously shown that negative feedback control is needed for homeostasis, and that feedback regulating replication probabilities is more effective than feedback regulating cell cycle lengths [71]. In ongoing work, we are investigating what effect including spatial structure and feedback from progeny will have on dedifferentiation times. The effect of spatial structure on mutation acquisition is still not fully resolved. Some groups argue that time to acquire M~2 mutations is actually decreased in a spatial model compared to the space-free model [78]. Other groups find that time to multiple mutations is increased when space is considered [79,80]. Like other mathematical models, our model suggests that eradication of cancer is dependent on eradication of cancer stem cells [81][82][83]. The potential for progenitor cells to dedifferentiate and repopulate the stem cell compartment, however, may complicate successful treatment. Our work suggests that further progress in understanding initiation and treatment of cancer requires a more detailed understanding dedifferentiation and of stem cell homeostasis. (S7), and (B,D) global competition between subpopulations given by Eq. (S8). Bottom: Corresponding plots of total cell density. Basal dynamics are constant death rate m~1 and sigmoidal birth rate with maximal growth rate b 0~2 , b i~( 1zs)b i{1 for i~1,2. The same carrying capacity is used for all simulations: N 1~2 00N sc , N 2~2 50N sc , N 3~3 00N sc . Note that there is a sharp transition zone at which mutant cells go from nearly zero fraction of total population to majority of the differentiating cell population. However, the mutation rate u and proliferative advantage s at which this is observed is unreasonably high, just as for the model without progenitor competition (Fig. 2). (TIF) Text S1 Analytic solutions and derivations, alternative models, and Matlab code.

(PDF)
Author Contributions