An evolutionary explanation for the presence of cancer nonstem cells in neoplasms

Contrary to conventional views that assume all cells in a neoplasm can propagate the tumor, the cancer stem cell hypothesis posits that only a fraction of the cells (the cancer stem cells) can act as tumor-propagating cells, while most of the tumor is composed of cells with limited replication potential. Here, we offer an evolutionary approach to this controversy. We used several evolutionary, computational models to investigate cancer cell dynamics and conditions consistent with the stem cell hypothesis. Our models predict that if selection acts at the cell level, neoplasms should be primarily comprised of cancer stem cells, in contrast to experimental data indicating that neoplasms contain large fractions of cancer nonstem cells. We explore several solutions explaining the paradoxical existence of cancer nonstem cells in neoplasms, including the possibility that selection acts at the level of multicellular proliferative units.


Supplementary Figures
. The effect of varying the probability for symmetric division and the number of differentiation stages on the final percent of tumor-propagating cells and the proportion of tumor-propagating cell symmetric divisions. (A) In the absence of evolution, increasing the probability for symmetric division is sufficient to increase the percent of tumor-propagating cells in a neoplasm. Simulations in which the probability of symmetric division could not evolve show that high levels of symmetric division result in high numbers of tumorpropagating cells. Values of less than 0.5 for the probability of symmetric division led to neoplasm extinction (n=[0,0,0,50,50,50,50] for probability of symmetric division= [0.05,0.2,0.35,0.5,0.65,0.8,0.95]). The probability of symmetric division and the final percent of tumor-propagating cells is highly correlated (Pearson's correlation R=0.97). (B) In the absence of evolution, increasing the number of differentiation stages decreases the percent of tumor-propagating cells in the neoplasm for populations with a fixed symmetric division probability of 0.5 (n= [49,50,49,26] for number transient cell stages= [2,4,8,16]). The number of transient stages and the final percent of tumor-propagating cells are negatively correlated (Pearson's correlation R=0.88). (C) In evolving populations, the number of differentiation stages did not affect the final percent of tumor-propagating cells or (D) the proportion of tumor-propagating cell symmetric divisions reached by natural selection (n=[50,50,47,29] for number transient cell stages= [2,4,8,16]). The proportion of tumor-propagating cell symmetric divisions was significantly higher than control simulations conducted in the absence of evolution (p=[3x10 -164 ,2x10 -167 ,1x10 -153 ,3x10 -79 ] for number transient cell stages= [2,4,8,16]). Error bars show standard error of the mean, 50 simulations were initialized per condition. Replicate simulations from each condition were similar and led to small error bars; in some cases the error bars are indistinguishable from the top of the data bar. There was sufficient variability for the evolution of (A) high tumor-propagating cell and (B) symmetric division levels for the standard deviation of the symmetric division trait. The proportion of tumorpropagating cell symmetric divisions was significantly higher than control simulations conducted in the absence of evolution (p=[5x10 -28 ,2x10 -36 ,6x10 -120 ,2x10 -169 ,6x10 -166 ,4x10 -155 ,3x10 -169 ] for standard deviation=[0.01,0.02,0.04,0.08,0.16,0.32,0.64]). While the proportion of symmetric divisions did evolve under conditions of low trait variability, its rate of evolution was slow enough that there was insufficient time for the population to reach the high levels observed once it exceeded 0.04 (n=[48,47,47,48,47,48,48] for standard deviation=[0.01,0.02,0.04,0.08,0.16,0.32,0.64]). There was insufficient variability for the evolution of (C) high tumor-propagating cell and (D) symmetric division levels at low mutation rates as compared to control simulations conducted in the absence of evolution using a t-test with significance levels adjusted using the Bonferroni correction for multiple-testing (p=[0.08,0.002], n= [48,48] for mutation rate=[1x10 -7 ,1x10 -6 ]). . However, the mutation rate did not affect the final percent of tumor-propagating cells or the proportion of tumor-propagating cell symmetric divisions reached by natural selection once it exceeded 1x10 -4 mutations in the asymmetric / symmetric division gene set per cell division (p=[4x10 -21 ,1x10 -112 ,3x10 -161 ] n= [48,48,50] for mutation rate=[1x10 -5 ,1x10 -4 ,1x10 -3 ]). Together, the mutation rate and the standard deviation of new mutations in the probability of symmetric division can define whether the symmetric division trait can evolve quickly enough for tumor-propagating cells to take over the neoplasm. Due to computational limitations, we did not run the simulation at realistic population sizes (10 9 -10 12 ), however, with larger populations, lower mutation rates result in the evolution of high symmetric division probabilities and tumor-propagating cells will continue to dominate the neoplasm. Error bars show standard error of the mean, 50 simulations were initialized per condition.    S5. A graphical representation of the transitions between stages associated with the analytical stagestructured population model. Tumor-propagating cells with high symmetric division rates, p h , are labeled 0H and those with lower symmetric division rates, p l , are labeled 0L. Transient cells, labelled 1-3, differentiate at cell division until they terminally differentiate and reach stage 4. The apoptosis rate p a and the cell division rate p c are constant across all stages. Note that one cell in transient stage s-1 divides to become 2 cells in stage s, which is not depicted graphically but is included in the population projection matrix A in the Materials and Methods.

Stage structured population model
In our analytical, deterministic approach, we constructed a stage-structured population model (Lefkovitch 1965;Caswell 2001) to characterize the stable stage distribution of differentiation states. To do so, we split the tumor-propagating cells into two groups: stages 0H and 0L and assigned stage 0H a higher symmetric division rate p h than stage 0L p l , such that p h > p l . We then observed the outcome on the asymptotic population composition of "competition" between clones of tumor-propagating cells with different symmetric division rates. Transient and differentiated cells were defined as depicted in Fig. S5.
We constructed a stage-classified matrix model, in which n(t+1) = A n(t), where n(t) is a vector of stage abundances at time t and A is the population projection matrix. We depict our model both graphically (Fig.  S5) and mathematically using the population projection matrix A, .
This accounts for both the rate of transition between stages and the number of individuals affected. A matrix entry a ij can be interpreted as the contribution of cell abundances from stage j at time t to the cell abundances in stage i at time t + 1. For example, the number of cells in stage 1 at time t+1 is given by, n 1 (t+1) = n oH (t) p c (1-p h ) + n oL (t) p c (1-p l ) + n 1 (t) (1 -p c -p a ). The rest of the parameters in A are as previously defined (Table S2), where p a is the rate of apoptosis per day and p c is the rate of cell division per day.
The dominant eigenvalue λ 1 of matrix A corresponds to the change in population size over time and is constant when λ 1 = 1. When λ 1 > 1, the tumor size increases and when λ 1 < 1, the tumor shrinks. The dominant eigenvalue λ 1 of the characteristic equation of A is λ 1 = 1 -p a -p c p h , since p h > p l . Thus, tumor growth depends on the apoptosis and cell division rates and only the higher symmetric division rate. Using our default parameter values for apoptosis and cell division (Table S2), the tumor size is constant when p h = 0.41176, grows exponentially when p h > 0.41176, and shrinks to extinction when p h < 0.41176. This is consistent with the observation in our agent-based model that small initial probabilities of symmetric divisions could not sustain a tumor population (Fig. S1).
The right eigenvector of λ 1 gives the asymptotic distribution of cells in the stage classes and in the general case is .
From this, we observe that the asymptotic population distribution value for stage 0L is 0 for all values of p a , p c , p h , and p l , and thus cells with the lower symmetric division rates do not survive. This is consistent with the evolutionary simulations from our individual-based model, in which mutant clones with lower symmetric division rates eventually reached extinction. The actual dynamics of this cell population depend on the initial conditions, in that there must be cells in stage 0H.
We can calculate the proportion of the neoplasm that comprises tumor-propagating cells from the aymptotic population distribution by dividing the fraction of the cells that are in stage 0H by the sum of the fraction of cells in any stage. Doing this, we find that the fraction of tumor-propagating cells in the neoplasm depends only on the symmetric division rate for fast-dividing tumor-propagating cells and is given by . Thus, the proportion of tumor-propagating cells in the neoplasm increases as the tumor-propagating cell's rate of symmetric division, p h , increases (Fig. S6).

Individual-based model design and implementation
These methods are presented in ODD protocol standard format for individual-based models (Grimm et al. 2006). Briefly, tumor-propagating cells and transient cells were modeled using NetLogo 4.0.4 (Wilensky 1999). Each day in the simulation, a cell had the opportunity to die, divide, or do nothing. Transient cells could divide symmetrically a parameterized number of times before undergoing apoptosis or senescence. Except where noted, transient cells could not dedifferentiate into tumor-propagating cells. The probability that a tumor-propagating cell divided symmetrically or asymmetrically was determined by a mutable cell intrinsic trait that could evolve.

Model overview
Purpose This individual-based model tests how the probability for symmetric division evolves in tumor-propagating cells, or cancer stem cells, given the cancer stem cell hypothesis of cancer.

Entities, state variables, and scales
Each cell in a simulation has values for the state variables defined in Table S2.

Process overview and scheduling
The model is broken into timesteps, each of which is equivalent to 1 day. The flow chart in Fig. S4 outlines the decisions that a cell makes each timestep. Details of the implementation are described in the Submodels below.
Adaptation Cells in this model are not explicitly encoded with adaptive behavior. Instead, increasing propensity of symmetric division implicitly increases the fitness of the cell by decreasing the rate at which cells become fully differentiated and are removed from the population.

Sensing
Cells are not encoded with any sensing ability.

Interaction
There are no direct interactions or communication between cells in the model.

Stochasticity
Cell division, cell death, and mutation are stochastic in the model. They are constrained by parameterized probabilities.

Observation
For all simulations, summary data are recorded at the end of the simulation. For certain simulations, summary data are written to a file every 10 days. These data include the number of tumor-propagating cells, the number of transient cells, the observed average propensity for symmetric division by tumorpropagating cells, and the time at which simulation crossed the 90% tumor-propagating cells threshold.

Model details
Initialization See Table S2 for  propensity for cell division p c and apoptosis p a . Tumor-propagating cells can divide symmetrically or asymmetrically. Both daughter cells from a symmetric division remain tumor-propagating cells. When a tumor-propagating cell divides asymmetrically, one daughter cell remains a tumor-propagating cell and the other differentiates into a transient cell (Fig. 1A-C). Initially, each cell i has probability p sd,i (0) = p sd of symmetric division if it is a tumor-propagating cell. Each cell's probability of symmetric division p sd,i (t) is allowed to evolve over time. Because a tumor-propagating cell can only divide symmetrically or asymmetrically, allowing the probabilities of symmetric division or asymmetric division to evolve is equivalent.

Submodels
Each submodel corresponds to a decision in the flowchart (Fig. S4).

Each cell can die
For most of the simulations, a cell dies from background mortality if it draws a random number from a uniform distribution between 0 and 1 that is less than the probability for apoptosis, p a .
In a subset of the simulations that test the effects of feedback from transient cells to stem cells, transient cells die as before while stem cells with higher numbers of clonal transient cells die less often than those with fewer. The number of clonal transient cells for a stem cell i at time t, or n c,i (t), is defined as the number of transient cells that arise from the most recent common stem cell ancestor (as such, the stem cell and its clonal transient cells have common mutational state). Specifically, a stem cells dies if it draws a random number from a uniform distribution between 0 and 1 that is less than p af (t), the probability for apoptosis under feedback at time t, where p af (t) is obtained from the sigmoid function , and r is an arbitrary constant (r = -0.0001). Varying negative values of r had no effect on the outcome of the simulations.
A cell that is fully differentiated dies. The number of differentiation stages is a parameter s to the model, and each cell i has a state variable that tracks its differentiation level at time t, s i (t). Thus, cells with s i (t) = s die. The qualitative results of the model do not change if terminally differentiated cells are allowed to remain as quiescent cells in the model and are slowly cleared by background apoptosis.

Each cell can divide
A cell divides if it draws a random number from the uniform distribution between 0 and 1 that is less than the probability for cell division, p c .

Allow dividing cells to mutate
A cell i acquires a mutation in its probability for its tumor-propagating symmetric cell division trait if it draws a number from a uniform distribution between 0 and 1 that is less than the mutation rate, µ. Once it gets a new mutation, the current probability for symmetric division for its daughters, cells j and k, are updated to p sd,j (t+1) = p sd,k (t+1) = N(p sd,i (t),σ), where the mean of the normal distribution is the parental probability of symmetric division, p sd,i (t), and σ is the standard deviation parameter. The updated probabilities for symmetric division are bounded by 0 and 1.

Decide daughters' differentiation level
Recall that the differentiation stage for a cell i at time t is stored in its local variable s i (t). Here, we determine the differentiation stage for daughter cells j and k. In most simulations, a tumor-propagating cell can divide symmetrically or asymmetrically. It divides symmetrically if it draws a random number from a uniform distribution between 0 and 1 that is less that its current probability for symmetric division p sd,i (t), and asymmetrically otherwise. When a tumor-propagating cell divides symmetrically, both of its daughters remain tumor-propagating cells, s j (t+1) = s k (t+1) = 0. When a tumor-propagating cell divides asymmetrically, one daughter remains a tumor-propagating cell, s j (t+1) = 0, and the other differentiates into a transient cell, s k (t+1) = 1.
A transient cell divides symmetrically. Most of the simulations were run under a normal, or forward-only differentiation regime, in which the progeny from transient cells were always more differentiated, s j (t+1) = s k (t+1) = s i (t) +1. In a subset of the simulations, transient cells were allowed to dedifferentiate according to a parameterized probability, such that both daughters became tumor-propagating cells, s j (t+1) = s k (t+1) = 0.

Impose limits on neoplasm size
If the total number of cells in the neoplasm N exceeds the maximum possible carrying capacity of the tissue n max , then N -n max cells are chosen at random for death until the total population is n max . The total number of cells in a simulation fluctuates over the course of each timestep, typically decreasing when the cell death submodels are executed and increasing when the cell birth submodels is executed.

Niche-determination of tumor-propagating cells
We extended the basic model to explore the effects of microenvironmentally-determined "stemness", in which the differentiation structure remains, and any cell that moves into the niche becomes a tumorpropagating cell. Tumor-propagating cell self-renewal is defined as a symmetric division in which both daughter cells remain tumor-propagating cells, and has a probability of occurring in cell i at time t of p sr,i (t).
Here, we implement the mutation and evolution of self-renewal just as we implemented mutation of symmetric division in the basic model. All cells are placed on the neoplasm at (0,0) in a 32x32 torus, and the niche is located at (0,0). At cell division, one daughter remains in place and the other moves 1 nichewidth in a random direction. Any transient cell i that moves into the niche becomes a tumor-propagating cell, such that s i (t+1) = 0. Both daughters for a tumor-propagating cell in the niche remain tumorpropagating cells. Tumor-propagating cells not in the niche can only divide symmetrically, and so both daughters j and k are either tumor-propagating cells s j (t+1) = s k (t+1) = 0 or transient cells s j (t+1) = s k (t+1) = 1. A tumor-propagating cell not in the niche self-renews if it draws a random number from a uniform distribution between 0 and 1 that is less that its current probability for self-renewing division p sr,i (t), and differentiates otherwise. The initial probability for tumor-propagating cell self-renewing division is zero for all cells. All other transient cells differentiate as previously defined.
Supplementary Tables   Table S1. Summary data for all of the experimental conditions simulated.
Mean final probability for symmetric division and percent of tumor-propagating cells. Fifty simulations were initiated for each experimental condition, including the default parameter values (see Table S2), control experiments in which the default parameters were used but mutation was not allowed, and varying each of the parameter values. The percent of tumor-propagating cells, the population's mean probability for symmetric division, and the time to 90% tumor propagating cells were recorded. The mean and standard error are reported for these metrics. Most parameter values resulted in symmetric division rates that were statistically greater than the control simulation cases (marked *). Significance was obtained using a t-test with the Bonferroni multiple-testing correction (significant at p < 0.0016). Due to the nature of stochastic simulation modeling, sometimes populations went extinct before the simulation completed. The number of simulations per experiment that survived for 10 years is reported.   The probability that a tumor-propagating cell i divides symmetrically at time t, once it divides

Supplementary Text
Natural selection continues to favor a high tumor-propagating cell symmetric division probability, regardless of the number of differentiation stages.
In normal hematopoiesis, there are at least 8 distinct maturation stages, or levels of differentiation, in the neutrophil lineage (Khanna-Gupta and Berliner 2005). It is not known if the levels of differentiation in AML mimic that of hematopoiesis, or how many times a transient cell in AML may divide before it terminally differentiates. Thus, we varied the number of times a cell could divide before terminal differentiation. In experiments with constant probability for tumor-propagating symmetric divisions (i.e., no evolution), we found that increasing the number of differentiation stages decreases the percent of tumorpropagating cells in the neoplasm (Fig. S1B). Showing the robustness of our model by varying the number of differentiation stages, we found that natural selection still favored high chances of tumor-propagating cell symmetric divisions (Fig. S1D), and thus high levels of tumor-propagating cells (Fig. S1C).
Natural selection continues to favor a high tumor-propagating cell symmetric division probability as long as there is sufficient variability in the tumor-propagating cell symmetric division trait. Mutation in the probability of symmetric division was modeled at cell division as a two-step process. Each cell could either mutate or not, and if it did then its current probability for symmetric division was updated by drawing a random number from a normal distribution centered on the parental cell's symmetric division rate with a parameterized standard deviation (see Supplementary Methods). Because the standard deviation of the tumor-propagating cell symmetric division rate is unknown, we tested a range of values. We also varied the mutation rate across physiologically reasonable ranges. We found that when the mutation rate or the effects of mutations are extremely small, the tumor-propagating cell's symmetric division probability trait does not evolve during the 10-year experiment. However, with sufficient variability in the symmetric division trait, natural selection continues to favor a high tumor-propagating cell symmetric division probability and high levels of tumor-propagating cells (Fig. S2).
Natural selection continues to favor a high tumor-propagating cell symmetric division probability, regardless of the maximum neoplasm size. To keep the model computationally tractable, default in silico neoplasms had resources to maintain 100,000 cells. When a neoplasm reached this carrying capacity, cells were randomly selected to undergo apoptosis until the carrying capacity was again met. We varied the carrying capacity of the neoplasm and found that the population size of the neoplasm has no effect on either the evolution of the probability of symmetric division, or the final proportion of tumor-propagating cells in the neoplasm (Fig. S3).