## Abstract

With rare exceptions, human tumors arise from single cells that have accumulated the necessary number and types of heritable alterations. Each such cell leads to dysregulated growth and eventually the formation of a tumor. Despite their monoclonal origin, at the time of diagnosis most tumors show a striking amount of intratumor heterogeneity in all measurable phenotypes; such heterogeneity has implications for diagnosis, treatment efficacy, and the identification of drug targets. An understanding of the extent and evolution of intratumor heterogeneity is therefore of direct clinical importance. In this article, we investigate the evolutionary dynamics of heterogeneity arising during exponential expansion of a tumor cell population, in which heritable alterations confer random fitness changes to cells. We obtain analytical estimates for the extent of heterogeneity and quantify the effects of system parameters on this tumor trait. Our work contributes to a mathematical understanding of intratumor heterogeneity and is also applicable to organisms like bacteria, agricultural pests, and other microbes.

HUMAN cancers frequently display substantial intratumor heterogeneity in genotype, gene expression, cellular morphology, metabolic activity, motility, and behaviors such as proliferation rate, antigen expression, drug response, and metastatic potential (Fidler and Hart 1982; Heppner 1984; Nicolson 1984; Campbell and Polyak 2007; Dick 2008). For example, a molecular and phenotypic analysis of breast cancer cells has revealed defined subpopulations with distinct gene expression and (epi)genetic profiles (Shipitsin*et al.* 2007). Heterogeneity and the existence of subpopulations within single tumors have also been demonstrated via flow cytometry in cervical cancers and lymph node metastases (Nguyen*et al.* 1993) as well as in leukemias (Wolman 1986). Virtually every major type of human cancer has been shown to contain distinct cell subpopulations with differing heritable alterations (Heppner 1984; Merlo*et al.* 2006; Campbell and Polyak 2007). Heterogeneity is also present in premalignant lesions; for instance, genetic clonal diversity has been observed in Barrett's esophagus, a condition associated with increased risk of developing esophageal adenocarcinoma (Maley*et al.* 2006; Lai*et al.* 2007).

Tumor heterogeneity has direct clinical implications for disease classification and prognosis as well as for treatment efficacy and the identification of drug targets (Merlo*et al.* 2006; Campbell and Polyak 2007). The degree of clonal diversity in Barrett's esophagus, for instance, is correlated with clinical progression to esophageal adenocarcinoma (Maley*et al.* 2006). In prostate carcinomas, tumor heterogeneity plays a key role in pretreatment underestimation of tumor aggressiveness and incorrect assessment of DNA ploidy status of tumors (Wolman 1986; Haggarth*et al.* 2005). Heterogeneity has long been implicated in the development of resistance to cancer therapies after an initial response (Geisler 2002; Merlo*et al.* 2006) as well as in the development of metastases (Fidler 1978). In addition, tumor heterogeneity hampers the precision of microarray-based analyses of gene expression patterns, which are widely used for the identification of genes associated with specific tumor types (O'Sullivan*et al.* 2005). These issues underscore the importance of obtaining a more detailed understanding of the origin and temporal evolution of intratumor heterogeneity.

To study the dynamics of intratumor heterogeneity, we construct and analyze a stochastic evolutionary model of an expanding population with random mutational fitness advances. Evolutionary models of populations with random mutational advances have been studied in the context of fixed-size Wright–Fisher processes for both finite and infinite populations (Gerrish and Lenski 1998; Park and Krug 2007; Park*et al*. 2010); Gerrish and Lenski (1998) studied the speed of evolution in a Wright–Fisher model with random mutational advances in the context of finite but large populations while Park*et al*. (2010) obtained accurate asymptotic approximations for the evolutionary dynamics of the population, following ideas presented in Park and Krug (2007). The latter work and references therein constitute a substantial exploration of the effects of random mutational advances in fixed-size populations. Our present work complements this research by exploring the effects of random mutational advances in the context of exponentially expanding populations and in particular the implications of these mutational advances on population heterogeneity. Models of exponentially expanding populations are appropriate for the study of situations arising during tumorigenesis, but are also applicable to other organisms undergoing binary replication such as bacteria, agricultural pests, and other microbes and pathogens. Bacterial populations, for instance, are diverging quickly in both genotype and phenotype, as studied by Lenski and colleagues. These investigators examined the dynamics of phenotypic evolution in populations of *Escherichia coli* that were propagated by daily serial transfer for 1500 days, yielding 10,000 generations of binary fission (Lenski*et al.* 1991; Lenski and Travisano 1994). The fitness of the bacteria improved on average by 50% relative to the ancestor, and other phenotypic properties, such as cell size, also underwent large changes. Similarly, single malaria isolates have been found to consist of heterogeneous populations of parasites that can have varying characteristics of drug response, from highly resistant to completely sensitive (Foley and Tilley 1997). These findings have implications for treatment strategies, as not all pathogen populations are sensitive to therapeutic interventions, and necessitate the study of diversity dynamics in growing populations of cells.

In this article, we consider an exponentially expanding population of tumor cells in which (epi)genetic alterations confer random fitness changes to cells. This model is used to investigate the extent of genetic diversity in tumor subpopulations as well as its evolution over time. The mathematical framework is based on the clonal evolution model of carcinogenesis, which postulates that tumors are monoclonal (*i.e.*, originating from a single abnormal cell) and that over time the descendants of this ancestral cell acquire various combinations of mutations (Merlo*et al.* 2006; Campbell and Polyak 2007). According to this model, genetic drift and natural selection drive the progression and diversity of tumors. Our work complements studies of the effects of random mutational fitness distributions on the growth kinetics of tumors (Durrett*et al.* 2010) and contributes to the mathematical investigation of intratumor heterogeneity (Coldman and Goldie 1985, 1986; Michelson*et al.* 1989; Kansal*et al.* 2000; Komarova 2006; Haeno*et al.* 2007; Schweinsberg 2008; Bozic*et al.* 2010; Durrett and Moseley 2010).

## MATERIALS AND METHODS

We consider a multitype branching process model of tumorigenesis in which (epi)genetic alterations confer an additive change to the birth rate of the cell. This additive change is drawn according to a probability distribution *v*, which is referred to as the mutational fitness distribution. Cells that have accumulated *i* ≥ 0 mutations are denoted as type-*i* cells. Initially, the population consists entirely of type-0 cells, which divide at rate *a*_{0}, die at rate *b*_{0}, and produce type-1 cells at rate *u*. The initial population, whose size is given by *V*_{0}, is considered to be sufficiently large so that its growth can be approximated by , where *λ*_{0} *= a*_{0} *− b*_{0} and time *t* is measured in units of cell division. Type-1 cells divide at rate *a*_{0} + *X*, where *X* ≥ 0 is drawn according to the distribution *v*, and give rise to type-2 cells at rate *u*. All cell types die at rate *b*_{0} < *a*_{0}. In general, a type-(*k* − 1) cell with birth rate *a _{k−}*

_{1}produces a new type-

*k*cell at rate

*u*, and the new type-

*k*cell type divides at rate

*a*

_{k}_{−1}+

*X*, where is drawn according to

*v*. Each type-

*k*cell produced by a type-(

*k*− 1) cell initiates a genetically distinct lineage of cells, and the set of all of its type-

*k*descendants is referred to as its family. The total number of type-

*k*cells in the population at time

*t*is given by and the set of all type-

*k*cells is called the

*k*th wave or generation

*k*. The total population size at time

*t*is given by .

Note that in our model, mutations are not coupled to cell divisions but instead occur at a fixed rate per unit of time. This model can be modified to exclusively consider (epi)genetic alterations that occur during cell division by assuming that type-0 individuals divide at rate *α*_{0} and during each division, there is a chance *μ* that a mutation occurs, producing a type-1 individual. These two model versions are equivalent for wave-1 cells if *a*_{0} = *α*_{0}(1 − *μ*) and *u* = *α*_{0}*μ*. For later waves, the model must be altered so that the mutation rate is dependent upon the genetic constitution of the cells, since accumulation of alterations modifies the fitness of the cell and hence the rate at which further changes are accumulated. Analysis of this model would then require the mutation rate term to be inside the integral over the support of the fitness distribution; however, this modification does not alter the limiting results significantly. Our model and results can also be modified to account for generation-dependent mutation rates (Coldman and Goldie 1985; Park*et al.* 2010).

The mutational fitness distribution *v* determines the effects of each (epi)genetic change that is accumulated in the population of cells. We consider fitness distributions concentrated on [0, *b*] for some *b* > 0. We discuss two distinct classes of distributions: (i) *v* is discrete and assigns mass *g _{i}* to a finite number of values ; (ii)

*v*is continuous with a bounded density that is continuous and positive at

*b*. Figure 1 shows a snapshot of the population decomposition in a sample simulation for case ii and illustrates the complex genotypic composition of the population of cells generated by our model.

The determination of the distribution *v* in the context of bacteria and viruses has been the subject of several experimental studies (Imhof and Schlotterer 2001; Sanjuan*et al.* 2004), which have generally produced results leading to the conclusion that *v* has an exponential distribution. However, more recently Rokyta and colleagues (Rokyta*et al.* 2008) presented studies of bacteriophages, in which the distribution of beneficial mutational effects appears to have a truncated right tail. One possible explanation for this result is that the experiments were done at an elevated temperature that might have led to a limited number of available beneficial mutations. A similar scenario might arise during tumorigenesis when only a limited number of (epi)genetic alterations allow a cell to progress to a more aggressive phenotype.

## RESULTS

There are two sources of heterogeneity present in the population: variability in the number of mutations per cell (heterogeneity between generations) and genotypic variation between members of the same generation (heterogeneity within a generation). We investigate these two sources of heterogeneity and derive analytic results that quantify the relationship between model parameters—*e.g.*, mutation rate and mutational fitness distribution—and the amount of genotypic variation present in the population over time.

#### Between-generation heterogeneity:

Asymptotic results for the size of generation *k* were obtained in (Durrett*et al.* 2010; Durrett and Moseley 2010); see Equations A1 and A2 in the appendix for restatements of the relevant results from these articles. Equation A1 implies that in case i, for example, we have the approximation(2)when *t* is large, where is the maximum growth rate that can be attained by generation *k* mutants,and is a positive random variable with known Laplace transform. In case ii there is an additional term in (2) of the form . Dividing both sides of (2) by and speeding up time by a factor of *L*, we note that the log size of generation *k* approaches a deterministic, linear limit as the mutation rate becomes very small. In particular, as , we have(3a)in probability, where(3b)The limiting process depends on *λ*_{0}, the growth rate of type-0 cells, and *b*, the maximum attainable fitness increase, but is otherwise independent of the particular choice of fitness distribution. An example of the limiting process is displayed in Figure 2, a and b.

As a consequence of Equation 3, we obtain the following insight regarding the birth time of type-*k* cells: if is the first time a type-*k* individual is born, then, as ,(4)in probability for all *k* ≥ 0. From the definition of , we have that is decreasing so that the increments between the birth times for successive generations decrease as *k* increases. This effect leads to an acceleration in the rate at which new mutations are accumulated (Figure 2c). This acceleration occurs regardless of the choice of fitness distribution, assuming that i or ii holds. Since is inversely proportional to *b*, distributions that allow for larger fitness increases tend to exhibit shorter increments.

Bozic*et al.* (2010) observed a similar acceleration of waves in their model of mutation accumulation. On the basis of approximations by Beerenwinkel*et al.* (2007), they concluded that this acceleration was an artifact of the presence of both passenger and driver mutations and that it does not occur when only driver mutations conferring a fixed selective advantage are considered—*i.e.*, when the fitness increments are deterministic. In contrast to these conclusions, we find that the acceleration of waves occurs regardless of the choice of mutational fitness distribution and is due to the difference in growth rates between successive generations: type-*k* cells arise when generation *k* − 1 reaches size and since the asymptotic growth rate of generation *k* is larger than the asymptotic growth rate of generation *k* − 1, generation *k* + 1 reaches size faster than generation *k*. We observed another example of this phenomenon earlier during our study of a related Moran model for tumor growth in which the total population of cells grows at a fixed exponential rate (Durrett and Mayberry 2010). In this model, the cause of acceleration was similarly related to growth rates—later generations take longer to achieve dominance in the expanding population of cells, and hence new types are born with a higher fitness advantage compared to the population bulk, allowing them to grow more rapidly.

As a second application of Equation 3, we derive an analytic expression for the time at which type-*k* cells become dominant in the population. Let for all be the first time that type-*k* cells become the most frequent cell type in the population. Then, as , we have(5)in probability for all The limit is the solution to , *i.e.*, the time when first overtakes . Figure 2d demonstrates how the index of the largest generation at time *t*, defined as changes over time. The transitions between periods of dominance are sharp only in the small mutation limit. At any given time, the population consists primarily of members of the current dominant generation; *i.e.*, (1/*L*)log *Z*(*Lt*)→ *z _{m}*

_{(t)}(

*t*) as . Therefore, for small mutation rates, the amount of genetic heterogeneity present in the population is determined by the amount of heterogeneity present in the dominant generation.

To determine the accuracy of our results when the mutation rate is small, we compare our limiting approximation in Equation 3 with the results of numerical simulations, using a mutational fitness distribution corresponding to a point mass at *b* [*i.e.*, ]. Given that there are ∼3 billion base pairs in the human genome and the mutation rate per base pair is (Seshadri*et al.* 1987; Oller*et al.* 1989; Kunkel and Bebenek 2000), point mutations occur at a rate of 0.3–30 per cell division. However, since advantageous mutations constitute only a fraction of all possible mutations, the mutation rate per cell division for advantageous mutations is smaller than the overall mutation rate per cell division. In the following example, we use a mutation rate per cell per time unit of and a cell division rate of *a* = 0.2. In Figure 3, we compare the average size of the *k*th generation in simulations (Figure 3a) with the approximation given by Equation 3 (Figure 3b). Although the behavior is qualitatively similar to the small mutation limit, the approximation consistently underestimates the times at which new waves appear. To explain the source of this bias, we use the alternative approximation given by the right-hand side of Equation 2, which can be rewritten as,(6)where . Using the expression for the Laplace transform of and the numerical algorithm presented in Ridout (2008), we sample 1000 variates from the distribution of . Table 1 shows the sample mean and standard deviation of log , for . The distribution of log has a positive mean and is skewed to the right (Figure 4), implying that the limit in Equation 3 in general underestimates the size of generation *k* for positive mutation rates. The approximation obtained by replacing log with the sample mean of log is displayed in Figure 3c. After an initial period in which the number of type-*k* cells is small, the behavior of the process closely resembles the one shown in Figure 3a. The variance of log increases with *k* and hence, we expect an increasing amount of variability in the simulations in the time when type-*k* cells arise. Figure 3d displays the right-hand side of Equation 6, replacing log with the value two standard deviations above its mean to illustrate an extreme scenario.

While the limit in (3) depends on the fitness distribution only through the maximum attainable fitness increase *b*, the distribution of also depends on the fitness distribution through the probability of attaining a fitness advance of *b* if *v* is discrete and the value of the probability density function at *b* if *v* is continuous (see Equation 4.9 in Durrett*et al.* 2010). As a consequence, our finite time approximation (6) takes into consideration the shape of the fitness distribution near *b* and the corrector term log accounts for variations in the likelihood of attaining the maximum possible fitness advance.

#### Within-generation heterogeneity:

We begin our investigation of within-generation heterogeneity by examining the extent of diversity present in the first generation of cells. We use two statistical measures to assess heterogeneity: (i) Simpson's index, which is given by the probability that two randomly chosen cells from the first generation stem from the same clone, and (ii) the fraction of individuals in the first generation that stem from the largest family of cells. To obtain these results, we derive an alternate formulation of the limit in Equation 1 that shows the limit is the sum of points in a nonhomogeneous Poisson process (see the appendix for more details). Each point in the limiting process represents the contribution of a different mutant lineage to so that it suffices to calculate i and ii for the points in the limiting process.

#### Simpson's index:

Let us introduce some terminology by defining to be the largest point in the limiting point process and by setting . Then Simpson's index for the point process is defined byWe may also consider Simpson's index for a random walk , where the are independent random variables with a tail probability , and . Then the limit as *n* goes to infinity of is (Fuchs*et al.* 2001), where denotes the ratio of the growth rate of type-0 cells to the maximal growth rate of type-1 cells. Furthermore, the expected value of converges to the expected value of *R* so we have(7)See the appendix for details of this calculation. Equation 7 shows that the average amount of heterogeneity present in the first generation depends only on *α*. Figure 5 displays the sample mean of Simpson's index of *Z* as a function of time for different values of *α*. Initially, the sample mean tends toward the expected value of Simpson's index for the limiting point process. Although the sample mean is greater than the limiting value for larger values of *b*, our theory guarantees that eventually, the values of the sample mean converge. However, it is impossible to simulate the process for such long times since the population size and number of cell types become too large.

Expressions for the density and higher moments of *R* can be obtained in a similar manner by the comparison techniques we used in the proof of Equation 7 (see appendix and Logan*et al*. 1973; Shao 1997). In addition, near the origin, the density *g*(*x*) of *R* has the form(8)as _{.} Figure 6 displays simulations of Simpson's index for the first wave mutants in the branching process *Z*. Note the convergence of the empirical distribution of Simpson's index to the distribution for the limiting point process.

#### Largest clones:

To further investigate heterogeneity properties of the point process, we examine the fraction of cells descended from the largest family of first-generation mutants defined as . This quantity reveals the degree of dominance of the largest clone in the first wave of mutants. For large *n*, values of *V _{n}* near one indicate that the population is largely dominated by a single clone, while values near zero indicate a highly heterogeneous population where no single clone contributes significantly to the total. Using a similar approach to that in the previous section, we again consider this calculation in the context of a random walk. Consider a sequence of independent, identically distributed random variables

*Y*with partial sums

_{i}*W*. Define the maximum value of this sequence from 1 to

_{n}*n*to be

*Y*

_{(1)}. Then classical results regarding one-dimensional random walks characterize the limiting characteristic function of

*W*/

_{n}*Y*

_{(1)}(Darling 1952). In the appendix, we demonstrate that these results can be applied to the study the largest clone contributions in the limiting point process.

We show that 1/*V _{n}* converges in distribution to a nontrivial limit

*W*and obtain an explicit formula for the characteristic function of the limit: as , , where

*W*has characteristic function

*ψ*satisfying

*ψ*(0) = 1 and(9a)with(9b)Interestingly, the characteristic function of

*W*is nonintegrable since its density has a singularity at 1. This finding implies that there is a disproportionately large chance that a single clone dominates the population. Further details are shown in the appendix.

Differentiating *ψ* then leads to simple expressions for the mean and variance of the limit,(10)Figure 7a suggests that the rate of convergence is slow for *α* close to 1. These formulas provide us with the first two moments of the diversity measure *W* and reveal its dependence on *α*.

Equation 9 implies that *V _{n}* converges to a nontrivial limit

*V*=

*W*

^{−1}and Jensen's inequality applied to the strictly convex function 1/

*x*implies that This result provides a lower bound on the expected value of the limit of

*V*and indicates that for values of

_{n}*α*close to zero, the population is eventually dominated by a single clone. Even though this result indicates only a lower bound, simulations suggest that deviations of the mean from 1−

*α*are small, as illustrated in Figure 7b.

#### Extensions to generation *k*:

The results obtained in the previous two sections for the first generation of cells can easily be extended to later generations by noting that each mutation to a type-(*k −* 1) individual starts a new family regardless of the mutated cell's family tree. Therefore, the amount of heterogeneity within generation *k* depends on the relative growth rates of type-*k* and type-(*k* − 1) individuals in the same way that the amount of heterogeneity within generation 1 depends on the relative growth rates of type-1 and type-0 individuals. In particular, letting *R _{k}* denote the limiting value of Simpson's index for generation

*k*, we have(11)where is the relative growth rate of type-(

*k*−1) individuals compared to type-

*k*individuals (see also Equation 7 above). Note that the relative growth rates

*α*are increasing functions of

_{k}*k*. Therefore, Equation 11 shows that the mean of Simpson's index is a decreasing function of the generation number.

#### Total population heterogeneity:

We conclude this section by demonstrating how heterogeneity measures for the population as a whole can be calculated using information about both intra- and interwave heterogeneity. First, define the total collection of all genotypes present at time *t* as *N*(*t*). For a particular genotype, *x*, the number of cells at time *t* that have exactly this genotype is given by *Z*^{(x)}(*t*). Two cells have the same genotype if they contain exactly the same collection of mutations. Then Simpson's index for the entire population is given byTo show how this expression depends on the contributions of different waves, define the total collection of genotypes present in the wave-*k* population at time *t* by *N _{k}*(

*t*), and let

*Z*

_{k}^{(x)}be the population of cells in wave

*k*that have genotype

*x*. By defining

*K*(

*t*) to be the number of waves present at time

*t*, we obtain the alternate expressionThis decomposition expresses Simpson's index for the whole population in terms of Simpson's index for each wave and the contribution of each wave to the total population. Combining the result in (5), which gives the dominant wave as a function of time, with (11), which describes wave-

*k*heterogeneity, we obtain a description of how the extent of heterogeneity changes in time. However, more refined results are needed to describe the transitions between the dominance of successive waves.

## DISCUSSION

In this article, we investigated the evolution of intratumor heterogeneity in a stochastic model of tumor cell expansion. Our model incorporates random mutational advances conferred by (epi)genetic alterations and our analysis focused on the extent of heterogeneity present in the tumor. We first considered heterogeneity between tumor subpopulations with varying numbers of alterations and obtained limiting results, as the mutation rate approaches zero, for the contribution of each wave of mutants to the total tumor cell population. We showed that in the limit, this intergeneration heterogeneity depends on the maximum attainable fitness advance conferred by (epi)genetic alterations, but not on the specific form of the mutational fitness distribution. Our analysis also led to analytical expressions for the arrival time of the first cell with *k* mutations and showed that the rate of accumulation of new genetic alterations accelerates over time due to the increasing growth rates of successive generations. We demonstrated with stochastic simulations that for small but positive mutation rates, our limiting approximations provide good predictions of the model behavior (see Figure 3). These simulations also suggest that as time increases, multiple waves of mutants coexist without a single, largely dominant wave. For large *t*, the mean growth rate of the *k*th wave is given by , showing that variation in fitness within a particular wave is a transient property. The extent to which this variation affects tumor dynamics at small times is the subject of ongoing work.

We also investigated the genotypic diversity present within the *k*th generation of mutants by considering two measures of diversity: Simpson's index, which is given by the probability that two randomly selected cells stem from the same family, and the fraction of individuals in generation *k* that stem from the largest family of individuals. We obtained limiting expressions for the mean of Simpson's index as well as the form of its density near the origin. Interestingly, the limiting mean of Simpson's index is given by the quantity 1−*α*, where *α* is the ratio between the maximum attainable fitness values of type-(*k* − 1) and type-*k* individuals. We then observed that, as time increases, the mass of the distribution of Simpson's index moves closer to 0, indicating higher levels of diversity in the tumor at later times (see Figure 6). This behavior was also observed via direct numerical simulation of the branching process—the distribution and mean of Simpson's index converged to the predicted limiting values.

Finally, we investigated the ratio between the total population size of the *k*th wave of mutants and the size of the largest family. We showed that this ratio can be approximated by a random variable with mean (1−*α*)^{−1}. An explicit formula for the characteristic function of this random variable was also obtained (see Equation 9). Note that as α approaches 1, the mean of the ratio grows to infinity—*i.e.*, the largest family of cells constitutes a vanishing proportion of the total population of wave-*k* cells as the maximum possible fitness advance goes to zero.

In the context of tumorigenesis, where a tumor originally starting from a single cell reaches cell numbers of ≥10^{12}, the limit as *t* goes to infinity is indeed the appropriate regime of study. We have compared our results regarding Simpson's index with numerical simulations at finite times (see, *e.g.*, Figures 5 and 6) and have found good qualitative agreement with the asymptotic limits. Estimates of overall mutation rates in evolving cancer cell populations range from 10^{−9} to 10^{−2}. In Figure 3 we demonstrate good agreement between our results regarding interwave heterogeneity in the small mutation rate limit and numerical simulations when the mutation rate is ; thus, analysis in this limiting regime captures the dynamics for mutation rates < per cell division.

In conclusion, our analysis indicates that tumor diversity is strongly dependent upon the age of the tumor and the maximum attainable fitness advance of mutant cells. If only small fitness advances are possible, then the tumor population is expected to have a larger extent of diversity compared to situations in which fitness advances are large. The acceleration of waves observed in our studies of intergeneration heterogeneity provides evidence that an older tumor has a higher level of diversity than a young tumor. In addition, we have shown that the mean of Simpson's index for generation *k* is a decreasing function of the generation number (see Equation 11), indicating a larger extent of diversity in later generations and suggesting a further increase in the total extent of heterogeneity present in the tumor at later times.

Possible extensions of our model include spatial considerations and the effects of tissue organization on the generation of intratumor heterogeneity as well as the inclusion of other cell types, such as immune system cells and the microenvironment. Furthermore, alternative growth dynamics should be considered to test the extent of heterogeneity arising in populations that follow logistic, Gompertzian, or other growth models. We have neglected these aspects in the current version of the model to focus on the dynamics of tumor diversity in an exponentially growing population of cells. Our model provides a rational understanding of the extent and dynamics of intratumor heterogeneity and is useful for obtaining an accurate picture of its generation during tumorigenesis.

## Acknowledgments

The authors thank Theresa Edmonds for exceptional research assistance. This work is supported by National Cancer Institute grant U54CA143798 (to J.F., K.L., and F.M.), National Science Foundation (NSF) grant DMS 0704996 (to R.D.), and NSF RTG grant DMS 0739164 (to J.M.).

## APPENDIX

#### Convergence in distribution of *Z*_{k}(*t*):

_{k}

Define *λ _{k}* =

*λ*

_{0}+

*kb*to be the maximum growth rate that can be attained by a generation

*k*mutant, and letIf

*ν*is discrete and assigns mass

*g*to a finite number of values

_{i}*b*

_{1}<

*b*

_{2}< … <

*b*=

_{N}*b*, then(A1)where “⇒” denotes convergences in distribution, and

*V*

_{d,k}has Laplace transformfor all

*θ*≥ 0.

If instead *ν* is a continuous distribution on [0, *b*] with a bounded density *g* that is continuous and positive at *b*, then(A2)where *V*_{c,k} has Laplace transform for all *θ* ≥ 0. Here *d _{k}*(

*λ*

_{0},

*b*) and

*c*(

_{k}*λ*

_{0},

*b*) are constants that depend on the model parameters (Durrett

*et al.*2010), and “d” and “c” in the constants and subscripts stand for discrete and continuous. See Durrett and Moseley (2010) for a proof of Equation A1 and Durrett

*et al.*(2010) for a proof of (A2). To show the dependence on the density

*g*, we write the value of

*c*

_{1}(

*λ*

_{0},

*b*) asThe convergence in (A2) was also numerically investigated in Figures 1 and 2 of Durrett

*et al.*(2010).

#### Point process limit:

In this section, we discuss the point process representation for the limit in Equations A1 and A2 in the case where *k* = 1. The limit is the sum of points in a nonhomogeneous Poisson process. Before stating this result, we introduce some terminology. Here and in what follows, we use to denote the number of points in the set *A*. We say that Λ is a Poisson process on (0, ∞) with mean measure *μ* if Λ is a random set of points in (0, ∞) with the following properties:

For any

*A*⊂(0, ∞),*N*(*A*) = is a Poisson random variable with mean*μ*(*A*).For any

*k*≥ 1, if*A*_{1}, … ,*A*are disjoint subsets of (0, ∞), then_{k}*N*(*A*), 1 ≤_{i}*i*≤*k*are independent.

We also let *α* = *λ*_{0}/*λ*_{1} ∈ (0, 1) denote the ratio of the growth rate of type-0 cells to the maximal growth rate of type-1 cells and note that 1 + *p*_{1} = 1/*α*.

Theorem 1. *Let* Λ *be a Poisson process on* (0, ∞) *with mean measure**and let S denote the sum of the points in* Λ*. Then positive constants A*_{d}, *A*_{c} *= A*_{d}(*λ*_{0}, *b*), *A*_{c}(*λ*_{0}, *b*) *exist that depend on the indicated parameters so that in case i as t* →∞*and in case ii as t* →∞

For more details, see Durrett*et al.* (2010), Theorem 3, and Durrett and Moseley (2010), Corollary to Theorem 3. Note that the mean measure for Λ has tail *m*(*x*, ∞) = *x*^{−α}.

Let *X _{n}* denote the

*n*th largest point in Λ, and let denote the sum of the

*n*largest points. To determine the dependence of

*X*on

_{n}*n*we first note that if we define Λ′ =

*f*(Λ) where

*f*(

*x*) =

*x*

^{−α}, then Λ′ is a Poisson process and after making the change of variables

*y*=

*x*

^{−α}, we can see that the mean measure isIn other words, Λ′ is a homogeneous Poisson process with constant intensity and hence the spacings between points are independent exponentials with mean 1. If we let

*T*denote the time of the

_{n}*n*th arrival in Λ′, then the law of large numbers implies that

*T*∼

_{n}*n*as

*n*→ ∞. Since , we obtain

*X*∼

_{n}*n*

^{−1/α}as

*n*→ ∞. In addition, we have the following Lemma:

Lemma 1.

*Furthermore*, *if we define* , *then*

*Proof*. Since *T _{n}* has a Gamma(

*n*, 1) distribution, we have Stirling's approximation implies that Γ(

*n*− 1/

*α*)/ Γ(

*n*) ∼

*n*

^{−1/α}and the second conclusion follows. ■

#### Simpson's index:

To prove Equation 7 in the text, we use a result in Fuchs*et al.* (2001) that shows that(A3)Herewhere *Y _{i}* are iid random variables in the domain of attraction of a stable law with index

*α*and

*S*=

_{n}*Y*

_{1}+ … +

*Y*. To explain the connection between the two results, note that if we have

_{n}*P*(

*Y*>

_{i}*x*) =

*x*

^{−α}, for

*x*≥ 1 and letting

*Y*

_{n}_{,i}=

*Y*/

_{i}*n*

^{1/α}, thenThis implies that if we let Δ

*= {*

_{n}*Y*

_{n}_{,i}:

*i*≤

*n*} be the point process associated with the

*Y*

_{n}_{,i}and define the measures and then we haveso that we should expect

*ER*to agree with lim

*ER*.

_{n}To make this argument rigorous, letdenote the truncated value of Simpson's index for Λ* _{n}* anddenote the truncated value for Λ. Then for any

*ε*> 0, we have(A4)We complete the proof by deriving appropriate bounds for each of the three terms on the right-hand side of (A4). For the first term, we have the following:Lemma 2.

*where h*→ 0

_{ε}*as ε*→ 0.

*Proof.* Let *ε* > 0 and write and for *k* = 1, 2. Sincefor *k* = 1, 2, we have the bound(A5)After noting that and we have for any *δ* > 0,(A6)where we have used (A5) in the last line. To control the third term on the right, let *φ* denote the Laplace transform of *Y*_{1}. Thenas *t* → 0 since 1 − *e*^{−x} ∼ *x* as *x* → 0 implies that We can conclude thatas *n* → ∞. In particular, *A _{n,}*

_{1}⇒

*A*

_{1}, where

*A*

_{1}has the above Laplace transform. Sinceas

*t*→ ∞, we have

*P*(

*A*

_{1}= 0) = 0 so that taking

*δ*=

*ε*

^{(1−α)/2}in (A6) yields the result. ■

To bound the second term on the right-hand side of (A4), we need some notation. Let *M _{p}* denote the class of all point measures on (0, ∞). In a slight abuse of notation, we write

*v*∈

*ν*when

*ν*∈

*M*and

_{p}*v*∈ supp(

*ν*). We equip

*M*with the topology of vague convergence (see, for example, Section 3.4 in Resnick 1987) and take as our

_{p}*σ*-algebra the one generated by open sets in this topology. Associated with any random set of points, we can associate a measure

*ξ*that is a random variable with values in

*M*

_{p}_{.}We write Λ

*⇒ Λ to mean that the associated random measures*

_{n}*ξ*⇒

_{n}*ξ*.Lemma 3. Λ

*⇒ Λ*

_{n}*and if we define the maps F*

_{k}_{,ε}:

*M*→ [0, ∞)

_{p}*by*

*for k =*1, 2,

*then*

*Proof*.

^{ }Sincefor all Borel sets

*A*, the first claim follows from Proposition 3.21 in Resnick (1987). The second claim follows from the continuous mapping theorem (see, for example, Resnick 1987, p. 152) the fact that

*F*

_{k}_{,ε}is continuous away from measures

*ν*with

*ε*∉

*ν*and the fact that the random measure associated with Λ has no point masses with probability 1. ■

As a consequence of this lemma, the fact that *R _{n}*(

*ε*) ≤ 1, and the bounded convergence theorem, we have the following:

Corollary 1.*as n* → ∞ *for any ε* > 0.

It thus remains to establish the following:

Lemma 4.

*Proof*. We can establish this result using the same results as in the proof of Lemma 2, in particular if we define and for *k* = 1, 2. Then, following the display in Equation A6, we have for any *δ* > 0It is obvious that *P*(*A*_{1} = 0) = 0 and for *ε* < 1, so it remains only to establish thatas *ε* → 0. This result follows immediately from Lemma 1. Therefore, taking completes the proof. ■

We can now complete the proof of Equation 7 by letting *n* → ∞ and then *ε* → 0 in (A4) and applying Lemmas 2 and 4 and Corollary 1.

To prove Equation 8, we use a result in Logan*et al.* (1973) that establishes that as *n* → ∞,has a limiting distribution with a density *f* that satisfiesfor some constants *a*, *b* > 0 (see Logan*et al.* 1973, Equation 5.7, and Shao 1997, Theorem 6.1). Making the change of variables *x* = *y*^{−1/2} yields Equation 8. ■

#### Largest clones:

Equation 9 is a consequence of the following theorem:Theorem 2.* As n* → ∞, *where W has characteristic function ψ satisfying ψ* (0) = 1 *and**for all t* ≠ 0 *with*

The form of the characteristic function is the same as the characteristic function for lim_{n}_{→∞}*T _{n}*/

*Y*

_{(1)}, where the

*Y*are iid random variables with power law tails,

_{i}*Y*

_{(1)}= max

_{i}_{≤n}

*Y*, and (see, for example, Darling 1952). Again, this agreement is a consequence of the previously discussed connection between Δ

_{i}*and the limiting Poisson point process.*

_{n}To prove Theorem 2, we need the following notation. For a real number *t*, we define the functionFor a complex number *z* we denote the real part of *z* by Re[*z*] and its imaginary part by Im[*z*].

*Proof of Theorem 2.* Theorem 5.1 in Darling (1952) implies that we haveas *n* → ∞ whereas in the *Simpson's index* section, *Y*_{(1)} = max_{i}_{≤n}*Y _{i}* and To conclude that

*T*/

_{n}*Y*

_{(1)}⇒

*V*, we need to show that

*y*is continuous at 0. To establish this fact, we make the change of variables

*v*=

*tu*to conclude that(A7)Since 1 − exp(

*iv*) ∼ −

*iv*as

*v*→ 0, the integral on the right-hand side of (A7) is finite and henceas

*t*→ 0. Since

*T*/

_{n}*Y*

_{(1)}⇒

*V*, the fact that

*S*/

_{n}*X*

_{1}⇒

*V*follows from the arguments in the previous section. ■

It is interesting to note that the characteristic function in Theorem 2 is not integrable. The problem is that the density of blows up near 1. As an explanation for this, we note that with probabilitythere is a point in the process bigger than x and no points in [1, x). When this happens,and soIf we had *F*_{n}(*y*) ? (*y* − 1)^{α}, then the density would blow up like (*y*−1)^{α−1} as *y* → 1. We confirm that this gives the right asymptotic by providing an explicit formula for the density of W.

Corollary 2. *W has a density on* (1, ∞) *given by*

Note that integral expression above does not converge absolutely so part of the proof consists of showing that the limit exists. If we apply the change of variable *s* = *t*(*y* − 1) in the definition of *f*, we see thatthus confirming the intuition that the density blows up like (*y* − 1)^{α}^{−1} as *y* approaches 1.

*Proof of Corollary 2*. We first establish that there are no point masses in the distribution of *V*. By the inversion formula we have for any *a* ∈ ℝ,If we focus on the positive axis and use the change of variable *s* = *t*/*T*,From display (A7) it follows that for every *s* ∈ (0, 1) we have *e ^{isT}*

^{(1−a)}/

*f*(

_{α}*sT*) → 0 as

*T*→∞. Note thatwhich implies |

*f*(

_{α}*t*)| ≥ 1 for all

*t*. Therefore |

*e*

^{isT}^{(1−a)}/

*f*(

_{α}*t*)| ≤ 1 for all

*t*and it follows via the dominated convergence theorem thatA similar result holds for the integral on the negative axis and we conclude thatWe can therefore conclude for

*x*> 1 and

*h*> 0 via the inversion formula (see Durrett 2005, Equation 3.2) and Fubini's theorem thatTherefore, to establish the result we need to show thatThis follows if we show that is a convergent integral and that a bounded function

*h*exists defined on (

*x*,

*x*+

*h*) such thatWe first use integration by parts to seeRecalling that |

*f*(

_{α}*T*)| → ∞ as

*T*→ ±∞, it follows that if we establish that is integrable on (−∞, ∞), then the convergence of the integral and the existence of a bounded dominating function will be established. Since

*f*is bounded away from 0, it suffices to check that the function decays fast enough. Recalling the definition of

_{α}*f*,which follows by passing the derivative inside the integral in the definition of

_{α}*f*. We can establish thatby observing,which can be found in many places,

_{α}*e.g.*, Loya (2005). Thus,We can similarly establish that for

*t*sufficiently largefor a positive finite constant

*C*

_{1}. Thus for

*t*sufficiently largeestablishing the result. ■

We conclude this section with the proof of Equation 10. Using the Taylor series expansion of exp(*iu*) ∼ 0 in (A7) above implies thatand thereforeso that in particular,for all *k* ≥ 1. Let *S*(*t*) = log*ψ*(*t*) = *it* − log *f _{α}*(

*t*). Then dropping the

*α*subscript on

*f*, we havewhich yields the desired result for the mean:Nowsocompleting the proof. ■

_{α}## Footnotes

Communicating editor: M. W. Feldman

- Received December 7, 2010.
- Accepted March 1, 2011.

- Copyright © 2011 by the Genetics Society of America