Intratumor Heterogeneity in Evolutionary Models of Tumor Progression

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₀, die at rate b₀, and produce type-1 cells at rate u. The initial population, whose size is given by V₀, is considered to be sufficiently large so that its growth can be approximated by $Z_0\left(t\right)=V_0\exp \left[{\lambda }_{0}t\right]$⁠, where λ₀  = a₀  − b₀ and time t is measured in units of cell division. Type-1 cells divide at rate a₀ + 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₀ < a₀. 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 $X\ge 0$ 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 $Z_k\left(t\right)$ and the set of all type-k cells is called the kth wave or generation k. The total population size at time t is given by $Z\left(t\right)$⁠.

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 α₀ 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₀ = α₀(1 − μ) and u = α₀μ. 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 $b_1<b_2<...<b_N=b$⁠; (ii) v is continuous with a bounded density $g\left(x\right)$ 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.

Figure 1.—

A sample cross-section of tumor heterogeneity. (A) The composition of a tumor cell population at $t=150$ time units after tumor initiation. Each “wave” of cells, defined as the set of cells harboring the same number of (epi)genetic alterations, is represented by a different color. Individual clones of cells with identical genotype are shown as circles, positioned on the horizontal axis according to fitness (i.e., birth rate) and on the vertical axis according to clone size. Note that later waves tend to have larger fitness values. (B) The time at which individual clones of cells were created during tumor progression. The color scale depicts the time of emergence of each clone. In A and B, parameters are $a_0=0.2$⁠, $b_0=0.1$⁠, $v\sim U\left(\right[0,0.05\left]\right)$⁠, and $u=\mathrm{0.001.}$

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

$$\log Z_k\left(t\right)\approx {\lambda }_{k}t-(k+p_k)\log (1/u)+\log V_{\text{d},k},$$

(2)

when t is large, where ${\lambda }_k={\lambda }_0+kb$ is the maximum growth rate that can be attained by generation k mutants,

$$p_k=-k+\frac{{\displaystyle {\sum }_{j=0}^{k-1}{\lambda }_k}}{{\lambda }_j},$$

and $V_{\text{d},k}$ is a positive random variable with known Laplace transform. In case ii there is an additional term in (2) of the form $(k+p_k)\log t$⁠. Dividing both sides of (2) by $L=\log (1/u)$ 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 $u\to 0$⁠, we have

$$(1/L){\log }^{+}{Z}_k\left(Lt\right)\to z_k\left(t\right)={[{\lambda }_{k}t-(k+p_k\left)\right]}^{+}={\lambda }_k{(t-{\beta }_k)}^{+}$$

(3a)

in probability, where

$${\beta }_k=\frac{k+pk}{{\lambda }_k}={\displaystyle \sum _{j=0}^{k-1}\frac{1}{{\lambda }_j}}.$$

(3b)

The limiting process depends on λ₀, 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.

Figure 2.—

The process for the small mutation limit. The limiting process is shown for the case in which the mutation rate goes to zero. (A) The time, t, on the horizontal axis vs. the log number of cells of wave k, $z_k\left(t\right)$⁠, as given by Equation 3. Both time and space are given in units of L = log(1/u). This plot shows the first 20 waves started at $t=10=1/{\lambda }_0,$  i.e., the time that type-1 cells begin to be born. (B) A closer look at the first 7 waves from a, showing the changes in the dominant type. (C) The birth times of the first 20 generations as a function of the generation number. (D) The dominant type in the population as a function of time. The index $m\left(t\right)$ of the largest generation at time t is defined as $z_{m\left(t\right)}\left(t\right)=\max \left\{z_k\right(t):k\ge 0\}$⁠. In A–D, parameters are ${\lambda }_0=0.1$ and $b=0.05$⁠.

As a consequence of Equation 3, we obtain the following insight regarding the birth time of type-k cells: if $T_k=\inf \{t\ge 0:Z_k(t)>0\}$ is the first time a type-k individual is born, then, as $u\to 0$⁠,

$$\frac{T_k}{L}\to {\beta }_k$$

(4)

in probability for all k ≥ 0. From the definition of ${\beta }_k$⁠, we have that ${\beta }_k-{\beta }_{k-1}=1/{\lambda }_k$ 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 $1/{\lambda }_k$ 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 $O(1/u)$ 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 $O(1/u)$ 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 $S_k=\inf \{t\ge 0:Z_k(t)>Z_j(t),$ for all $j\ne k\}$ be the first time that type-k cells become the most frequent cell type in the population. Then, as $u\to 0$⁠, we have

$$\frac{S_k}{L}\to t_k=b^{-1}+{\beta }_k$$

(5)

in probability for all $k\ge 1.$ The limit $t_k$ is the solution to ${\lambda }_k(t_k-{\beta }_k)={\lambda }_{k-1}(t_k-{\beta }_{k-1})$⁠, i.e., the time when $z_k\left(t\right)$ first overtakes $z_{k-1}\left(t\right)$⁠. Figure 2d demonstrates how the index $m\left(t\right)$ of the largest generation at time t, defined as $z_{m\left(t\right)}\left(t\right)=\max \left\{z_k\right(t):k\ge 0\},$ 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 $u\to 0$⁠. 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., $v\sim \delta \left(b\right)$]. Given that there are ∼3 billion base pairs in the human genome and the mutation rate per base pair is $O\left({10}^{-8}\right)-O\left({10}^{-10}\right)$ (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 $u={10}^{-5}$ and a cell division rate of a = 0.2. In Figure 3, we compare the average size of the kth 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

$${\hat{z}}_k^L\left(t\right)=L{z}_k(t/L)+\log V_{\text{d},k}$$,

(6)

where $L=\log (1/u)$⁠. Using the expression for the Laplace transform of $V_{\text{d},k}$ and the numerical algorithm presented in Ridout (2008), we sample 1000 variates from the distribution of $V_{\text{d},k}$⁠. Table 1 shows the sample mean and standard deviation of log $V_{\text{d},k}$⁠, for $k=1,2,3$⁠. The distribution of log $V_{\text{d},k}$ 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 $V_{\text{d},k}$ with the sample mean of log $V_{\text{d},k}$ 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 $V_{\text{d},k}$ 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 $V_{\text{d},k}$ with the value two standard deviations above its mean to illustrate an extreme scenario.

Figure 3.—

The size of the first four generations of cells. The log size of generations 1–4 is shown as a function of time t. Both time and space are plotted in units of $L=\log (1/u)$⁠. (A) The average values of the log generation sizes over ${10}^6$ sample simulations. (B) The limiting approximation from Equation 3 for the log size of the generations. (C) The approximation from Equation 6 using the mean of $\log V_{\text{d},k}$⁠. (D) The approximation from Equation 6 using a value two standard deviations above the mean of $\log V_{\text{d},k}$ to demonstrate an extreme scenario. Parameters are $u={10}^{-5}$⁠, $v\sim \delta \left(b\right)$⁠, $a_0=0.2$⁠, $b_0=0.1$⁠, and $b=\mathrm{0.05.}$

Figure 4.—

The distribution of log V_d,1. The relative frequency histogram of 1000 random samples from the distribution of $\log V_{\text{d},1}$ is shown, as specified by Equation 1.

TABLE 1

Means and standard deviations for log V_d,k for k = 1, 2, 3, as specified by Equation 1

Generation	Mean	Standard deviation
1	4.7638	1.3738
2	10.6010	2.2434
3	17.0519	3.0282

While the limit in (3) depends on the fitness distribution only through the maximum attainable fitness increase b, the distribution of $V_{\text{d},k}$ 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 $V_{\text{d},k}$ 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 $Z_1\left(t\right)$ 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 $X_n$ to be the $n\text{th}$ largest point in the limiting point process and by setting $S_n={\displaystyle {\sum }_{i=1}^n{X}_i}$⁠. Then Simpson's index for the point process is defined by

$$R=\frac{{\displaystyle {\sum }_{i=1}^{\mathrm{\infty }}{X}_i^2}}{{\left(S_{\mathrm{\infty }}\right)}^2}={{\displaystyle \sum _{i=1}^{\mathrm{\infty }}\left(\frac{X_i}{S_{\mathrm{\infty }}}\right)}}^2.$$

We may also consider Simpson's index for a random walk $R_n={{\displaystyle {\sum }_{i=1}^n\left(Y_i/W_n\right)}}^2$⁠, where the $Y_i$ are independent random variables with a tail probability $P(Y_i>x)=x^{-\alpha }$⁠, and $W_n={\displaystyle {\sum }_{i=1}^n{Y}_i}$⁠. Then the limit as n goes to infinity of $R_n$ is $1-\alpha$ (Fuchs  et al. 2001), where $\alpha ={\lambda }_0/{\lambda }_1\in (0,1)$ 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 $R_n$ converges to the expected value of R so we have

$$ER=1-\alpha .$$

(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 $1-\alpha ,$ 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.

Figure 5.—

The expected value of Simpson's index for the first wave of cells. The sample mean of Simpson's index (dots) over time t and the expected value of Simpson's index for the limiting point process (line) are shown, for two different values of $\alpha$⁠. Parameters are ${\lambda }_0=0.1$⁠, $a_0=0.2$⁠, and $v\sim U\left(\right[0,b\left]\right)$⁠, where $b=0.01$ in a and $b=0.05$ in b.

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

$$g\left(x\right)\sim a{x}^{-3/2}\exp [-b{x}^{-1}]$$

(8)

as $x\to 0$_.  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.

Figure 6.—

The empirical distribution of Simpson's index for the first wave of cells. Individual plots of Simpson's index are shown for the branching process at times $t=70$⁠, 90, 110, and 130 along with Simpson's index for the limiting point process $(t=\mathrm{\infty })$⁠. The histograms show the average over 1000 sample simulations. For the limiting point process, we approximate Simpson's index by examining the largest ${10}^4$ points in the process. Parameters are ${\lambda }_0=0.1$⁠, $a_0=0.2$⁠, and $v\sim U\left(\right[0,0.01\left]\right)$⁠.

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 $V_n=X_1/S_n$⁠. 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_i with partial sums W_n. Define the maximum value of this sequence from 1 to n to be Y₍₁₎. Then classical results regarding one-dimensional random walks characterize the limiting characteristic function of W_n/Y₍₁₎ (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 $n\to \mathrm{\infty }$⁠, $V_n^{-1}\Rightarrow W$⁠, where W has characteristic function ψ satisfying ψ(0) = 1 and

$$\psi \left(t\right)=\frac{\exp \left[it\right]}{f_{\alpha }\left(t\right)} \mathrm{for} \mathrm{all} t\ne 0,$$

(9a)

with

$$f_a\left(t\right)=1+\alpha {\displaystyle {\int }_0^1(1-e^{itu})u^{-(\alpha +1)}du.}$$

(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,

$$EW=\frac{1}{1-\alpha } \mathrm{and} \mathrm{var}\left(W\right)=\frac{2}{{(1-\alpha )}^2(2-\alpha )}.$$

(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 α.

Figure 7.—

The largest clones in the population of cells. (A) A comparison between Monte Carlo estimates for $E{V}_n^{-1}$ and the limit ${(1-\alpha )}^{-1}.$ (B) A comparison of the Monte Carlo estimates for $E{V}_n$ and the curve ${(1-\alpha )}^{+}$⁠. The Monte Carlo estimates are averaged over 100 sample simulations.

Equation 9 implies that V_n converges to a nontrivial limit V = W⁻¹ and Jensen's inequality applied to the strictly convex function 1/x implies that $E(\lim X_1/S_n)>1-\alpha .$ This result provides a lower bound on the expected value of the limit of V_n and indicates that for values of α 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

$$E{R}_k=1-{\alpha }_k,$$

(11)

where ${\alpha }_k={\lambda }_{k-1}/{\lambda }_k$ 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 α_k are increasing functions of 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 by

$$\text{SI}\left(Z\right(t\left)\right)={\displaystyle \sum _{x\in N\left(t\right)}{\left(\frac{Z^{\left(x\right)}\left(t\right)}{Z\left(t\right)}\right)}^2}.$$

To 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 expression

$$\text{SI}\left(Z\right(t\left)\right)={{\displaystyle \sum _{k=0}^{K\left(t\right)}{\displaystyle \sum _{x_k\in N_k\left(t\right)}\left(\frac{Z_k^{\left(x_k\right)}\left(t\right)}{Z_k\left(t\right)}\right)}}}^2{\left(\frac{Z_k\left(t\right)}{Z\left(t\right)}\right)}^2.$$

This 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 kth wave is given by ${\lambda }_0+kb$⁠, 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 kth 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 kth 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−α)⁻¹. 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¹², 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⁻⁹ to 10⁻². In Figure 3 we demonstrate good agreement between our results regarding interwave heterogeneity in the small mutation rate limit $(u\to 0)$ and numerical simulations when the mutation rate is ${10}^{-5}$⁠; thus, analysis in this limiting regime captures the dynamics for mutation rates <${10}^{-5}$ 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.

Acknowledgements

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.).

Footnotes

Communicating editor: M. W. Feldman

APPENDIX

Convergence in distribution of Z_k(t):

Define λ_k = λ₀ + kb to be the maximum growth rate that can be attained by a generation k mutant, and let

$$p_k=-k+\stackrel{k-1}{{\displaystyle \sum }_{j=0}}\frac{{\lambda }_k}{{\lambda }_j}.$$

If ν is discrete and assigns mass g_i to a finite number of values b₁ < b₂ < … <b_N = b, then

$${\left(1/u\right)}^{k+p_k}{e}^{-{\lambda }_{k}t}{Z}_k\left(t\right)\Rightarrow V_{\text{d},k},$$

(A1)

where “⇒” denotes convergences in distribution, and V_d,k has Laplace transform

$$\exp \left[-d_k\left({\lambda }_0,b\right)V_0{\theta }^{{\lambda }_0/{\lambda }_k}\right]$$

for all θ ≥ 0.

If instead ν is a continuous distribution on [0, b] with a bounded density g that is continuous and positive at b, then

$${\left(t/u\right)}^{k+p_k}{e}^{-{\lambda }_{k}t}{Z}_k\left(t\right)\Rightarrow V_{\text{c},k},$$

(A2)

where V_c,k has Laplace transform $\exp \left[-c_k\left({\mathit{\lambda }}_{\mathrm{0}}\mathrm{,} \mathit{b}\right)V_0{\mathit{\theta }}^{{\mathit{\lambda }}_{\mathrm{0}}/{\mathit{\lambda }}_{\mathit{k}}}\right]$ for all θ ≥ 0. Here d_k(λ₀, b) and c_k(λ₀, 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₁(λ₀, b) as

$$c_1\left({\lambda }_0,b\right)=g\left(b\right)\frac{{\lambda }_0+b}{{\lambda }_0}\cdot \frac{1}{{\lambda }_0+b}{\left(\frac{a_0+b}{{\lambda }_0+b}\right)}^{-b/\left({\lambda }_0+b\right)}\text{Γ}\left(\frac{{\lambda }_0}{{\lambda }_0+b}\right)\text{Γ}\left(1-\frac{{\lambda }_0}{{\lambda }_0+b}\right).$$

The 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 $\left|A\right|$ 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) = $\left|\Lambda \cap A\right|$ is a Poisson random variable with mean μ(A).

• For any k ≥ 1, if A₁, … , A_k are disjoint subsets of (0, ∞), then N(A_i), 1 ≤ i ≤ k are independent.

We also let α = λ₀/λ₁ ∈ (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/α.

Theorem 1. Let Λ be a Poisson process on (0, ∞) with mean measure

$$\mu \left(A\right)=\underset{A}{{\displaystyle \int }} \alpha z^{-\left(\alpha +1\right)}dz$$

and let S denote the sum of the points in Λ. Then positive constants A_d, A_c  = A_d(λ₀, b), A_c(λ₀, b) exist that depend on the indicated parameters so that in case i as t →∞

$${\left(A_{\text{d}}u{V}_0\right)}^{-\left(1+p_1\right)}{e}^{-{\lambda }_{1}t}{Z}_1\left(t\right)\Rightarrow S,$$

and in case ii as t →∞

$${\left(A_{\text{c}}u{V}_0\right)}^{-(1+p_1)}{t}^{1+p_1}{e}^{-{\lambda }_{1}t}{Z}_1\left(t\right)\Rightarrow S.$$

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 nth largest point in Λ, and let $S_n={\displaystyle {\sum }_{i=1}^n{X}_i}$ denote the sum of the n largest points. To determine the dependence of X_n on 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 is

$$\mu \prime \left(A\right)={\displaystyle {\int }_{f^{-1}\left(A\right)}\alpha x^{-\left(\alpha +1\right)}dx={\displaystyle {\int }_{A}dy=\left|A\right|}}.$$

In 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_n denote the time of the nth arrival in Λ′, then the law of large numbers implies that T_n ∼ n as n → ∞. Since $X_n=T_n^{-1/\alpha }$⁠, we obtain X_n ∼ n^−1/α as n → ∞. In addition, we have the following Lemma:

Lemma 1.

$$E{X}_n=\frac{\text{Γ}\left(n-1/\alpha \right)}{\text{Γ}\left(n\right)}.$$

Furthermore, if we define  $S_{\mathrm{\infty }}\text{=}{\displaystyle {\sum }_{\mathit{i}\text{=1}}^{\mathrm{\infty }}{X}_i}$⁠, then

$$E{S}_{\mathrm{\infty }}<\mathrm{\infty }.$$

Proof. Since T_n has a Gamma(n, 1) distribution, we have $E{X}_n=E{T}_n^{-1/\alpha }=\text{Γ}\left(n-1/\alpha \right)/\text{Γ}\left(n\right).$ 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

$$\underset{n\to \mathrm{\infty }}{\lim }E{R}_n=1-\alpha .$$

(A3)

Here

$$R_n=\stackrel{n}{{\displaystyle \sum }_{i=1}}{\left(\frac{Y_i}{S_n}\right)}^2,$$

where Y_i are iid random variables in the domain of attraction of a stable law with index α and S_n = Y₁ + … + Y_n. To explain the connection between the two results, note that if we have P(Y_i > x) = x^−α, for x ≥ 1 and letting Y_n,i = Y_i/n^1/α, then

$$nP\left(Y_{n,i}\in A\right)\to \mu \left(A\right).$$

This 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 ${\xi }_n\equiv {\displaystyle {\sum }_{x\in {\text{Λ}}_n}{\delta }_x}$ and $\xi ={\displaystyle {\sum }_{x\in \text{Λ}}{\delta }_x},$ then we have

$${\xi }_n\Rightarrow \xi ,$$

so that we should expect ER to agree with lim ER_n.

To make this argument rigorous, let

$$R_n \left(\mathrm{\epsilon }\right)=\frac{{{\displaystyle \sum }}_{i=1}^{\mathit{n}} Y_{n,i}^2{\mathbf{1}}_{Y_{n,i}>\mathrm{\epsilon }}}{{\left({{\displaystyle \sum }}_{i=1}^{\mathit{n}}{Y}_{n,i}^2{\mathbf{1}}_{Y_{n,i}>\mathrm{\epsilon }}\right)}^2}$$

denote the truncated value of Simpson's index for Λ_n and

$$R\left(\mathrm{\epsilon }\right)=\frac{{{\displaystyle \sum }}_i^{\mathrm{\infty }} X_i^2{\mathbf{1}}_{X_i>\mathrm{\epsilon }}}{{\left({{\displaystyle \sum }}_{i=1}^{\mathrm{\infty }} X_i{\mathbf{1}}_{X_i>\mathrm{\epsilon }}\right)}^2}$$

denote the truncated value for Λ. Then for any ε > 0, we have

$$\left|E{R}_n-ER\right|\le E\left|R_n-R_n\left(\mathrm{\epsilon }\right)\right|+E\left|R_n\left(\mathrm{\epsilon }\right)-R\left(\mathrm{\epsilon }\right)\right|+E\left|R\left(\mathrm{\epsilon }\right)-R\right|.$$

(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.

$$\underset{n\to \mathrm{\infty }}{\lim }\sup E\left|R_n-R_n\left(\mathrm{\epsilon }\right)\right|\le h_{\mathrm{\epsilon }},$$

where h_ε → 0 as ε → 0.

Proof. Let ε > 0 and write $A_{n,k}={\displaystyle {\sum }_{i=1}^n{Y}_{n,i}^k},$  $A_{n,k}\left(\mathrm{\epsilon }\right)={\displaystyle {\sum }_{i=1}^n{Y}_{n,i}^k{\mathbf{1}}_{Y_{n,i}>\mathrm{\epsilon }}},$ and ${\overline{A}}_{n,k}\left(\mathrm{\epsilon }\right)=A_{n,k}-A_{n,k}\left(\mathrm{\epsilon }\right)$ for k = 1, 2. Since

$$E{Y}_1^k{\mathbf{1}}_{Y_1\le \mathrm{\epsilon }{n}^{1/\alpha }}={\displaystyle {\int }_1^{\mathrm{\epsilon }{n}^{1/\alpha }}k{y}^{k-1}{y}^{-\alpha }dy\le C{\mathrm{\epsilon }}^{k-\alpha }{n}^{k/\alpha -1}}$$

for k = 1, 2, we have the bound

$$E{\overline{A}}_{n,k}\left(\mathrm{\epsilon }\right)\le C{\mathrm{\epsilon }}^{k-\alpha }.$$

(A5)

After noting that $A_{n,2}\le A_{n,1}^2,$  $A_{n,2}\left(\mathrm{\epsilon }\right)\le A_{n,1}^2\left(\mathrm{\epsilon }\right)\le A_{n,1}^2,$ and ${\overline{A}}_{n,1}\left(\mathrm{\epsilon }\right)+A_{n,1}\left(\mathrm{\epsilon }\right)=A_{n,1},$ we have for any δ > 0,

$$\begin{array}{l}E\left|R_n-R_n\left(\mathrm{\epsilon }\right)\right|=E\left|\frac{{\overline{A}}_{n,2}\left(\mathrm{\epsilon }\right)}{A_{n,1}^2}+R_n\left(\mathrm{\epsilon }\right)\left(\frac{A_{n,1}^2\left(\mathrm{\epsilon }\right)-A_{n,1}^2}{A_{n,1}^2}\right)\right|\\ \le {\delta }^{-2}E{\overline{A}}_{n,2}\left(\mathrm{\epsilon }\right)+P\left(A_{n,1}\le \delta \right)\\ +{\delta }^{-1}E\left(\frac{\left|2A_{n,1}\left(\mathrm{\epsilon }\right){\overline{A}}_{n,1}\left(\mathrm{\epsilon }\right){\overline{A}}_{n,1}^2\left(\mathrm{\epsilon }\right)\right|}{A_{n,1}}\right)+P\left(A_{n,1}\le \delta \right)\\ ={\delta }^{-2}E\left|{\overline{A}}_{n,2}\left(\mathrm{\epsilon }\right)\right|+{\delta }^{-1}E\left(\left|{\overline{A}}_{n,1}\left(\mathrm{\epsilon }\right)\left(2+\frac{A_{n,1}\left(\mathrm{\epsilon }\right)}{A_{n,1}}\right)\right|\right)+2P\left(A_{n,1}\le \delta \right)\\ \le {\delta }^{-2}E{\overline{A}}_{n,2}\left(\mathrm{\epsilon }\right)+3{\delta }^{-1}E{\overline{A}}_{n,1}\left(\mathrm{\epsilon }\right)+2P\left(A_{n,1}\le \delta \right)\\ \le \frac{{\mathrm{\epsilon }}^{2-\alpha }}{{\delta }^2}+\frac{3{\mathrm{\epsilon }}^{1-\alpha }}{\delta }+2P\left(A_{n,1}\le \delta \right),\end{array}$$

(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₁. Then

$$\begin{array}{c}1-\phi \left(t\right)=\alpha {\displaystyle {\int }_1^{\infty }\left(1-e^{-ty}\right)y^{-\left(\alpha +1\right)}}dy\\ =\alpha t^{\alpha }{\displaystyle {\int }_t^{\mathrm{\infty }}\left(1-e^{-x}\right)x^{-\left(\alpha +1\right)}}dx\sim C{t}^{\alpha }\end{array}$$

as t → 0 since 1 − e^−x ∼ x as x → 0 implies that ${\displaystyle {\int }_0^{\infty }\left(1-e^{-x}\right)x^{-\left(\alpha +1\right)}}dx<\mathrm{\infty }.$ We can conclude that

$$E\exp \left(-t{A}_{n,1}\right)=(1-{\left(1-\phi \frac{t}{n^{1/\alpha }}\right)}^n)\to \exp \left(-C{t}^{\alpha }\right)$$

as n → ∞. In particular, A_n,1 ⇒ A₁, where A₁ has the above Laplace transform. Since

$$1-\exp \left(-C{t}^{\alpha }\right)\to 1$$

as t → ∞, we have P(A₁ = 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_p and v ∈ supp(ν). We equip M_p with the topology of vague convergence (see, for example, Section 3.4 in Resnick 1987) and take as our σ-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 Λ_n ⇒ Λ to mean that the associated random measures ξ_n ⇒ ξ.Lemma 3. Λ_n ⇒ Λ and if we define the maps F_k,ε: M_p → [0, ∞) by

$$F_{k,\mathrm{\epsilon }}\left({\mu }_n\right)={\displaystyle \sum }_{x\in {\mu }_n}{x}^k{\mathbf{1}}_{x>\mathrm{\epsilon }}$$

for k = 1, 2, then

$$\left(F_{1,\mathrm{\epsilon }}\left({\text{Λ}}_n\right),F_{2,\mathrm{\epsilon }}\left({\text{Λ}}_n\right)\right)\Rightarrow \left(F_{1,\mathrm{\epsilon }}\left(\text{Λ}\right),F_{2,\mathrm{\epsilon }}\left(\text{Λ}\right)\right).$$

Proof.Since

$$nP\left(Y_{n,i}\in A\right)=n{\displaystyle {\int }_{n^{1/\alpha }A}\alpha y^{-\left(\alpha +1\right)}}dy={\displaystyle {\int }_A\alpha x^{-\left(\alpha +1\right)}}dx=\mu \left(A\right)$$

for 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.

$$E\left|R_n\left(\mathrm{\epsilon }\right)-R\left(\mathrm{\epsilon }\right)\right|\to 0$$

as n → ∞ for any ε > 0.

It thus remains to establish the following:

Lemma 4.

$$\underset{\mathrm{\epsilon }\to 0}{\lim }\sup E\left|R\left(\mathrm{\epsilon }\right)-R\right|=0.$$

Proof. We can establish this result using the same results as in the proof of Lemma 2, in particular if we define $A_k={\displaystyle {\sum }_{n=1}^{\mathrm{\infty }}{X}_n^k}$ and ${\overline{A}}_k\left(\mathrm{\epsilon }\right)={\displaystyle {\sum }_{n=1}^{\mathrm{\infty }}{X}_n^k{\mathbf{1}}_{X_n\mathrm{<\epsilon E}}}$ for k = 1, 2. Then, following the display in Equation A6, we have for any δ > 0

$$E\left|R-R\left(\mathrm{\epsilon }\right)\right|\le {\delta }^{-2}E{\overline{A}}_2\left(\mathrm{\epsilon }\right)+3{\delta }^{-1}E{\overline{A}}_1\left(\mathrm{\epsilon }\right)+2P\left(A_1\ge \delta \right).$$

It is obvious that P(A₁ = 0) = 0 and $E{\overline{A}}_2\left(\mathrm{\epsilon }\right)\le E{\overline{A}}_1\left(\mathrm{\epsilon }\right)$ for ε < 1, so it remains only to establish that

$$E{\overline{A}}_1\left(\mathrm{\epsilon }\right)\to 0,$$

as ε → 0. This result follows immediately from Lemma 1. Therefore, taking $\delta ={\left(E{\overline{A}}_1\left(\mathrm{\epsilon }\right)\right)}^{1/4}$ 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 → ∞,

$$S_n\left(2\right)=R^{-1/2}=\frac{{{\displaystyle \sum }}_{i=1}^n Y_i}{{\left({{\displaystyle \sum }}_{i=1}^n Y_i^2\right)}^{1/2}}$$

has a limiting distribution with a density f that satisfies

$$f\left(y\right)\sim a{e}^{-b{y}^2}, \mathrm{as}y\to \mathrm{\infty },$$

for 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 → ∞, $V_n^{\mathrm{-}\text{1}}\Rightarrow W,$  where W has characteristic function ψ satisfying ψ (0) = 1 and

$$\mathit{\psi }\left(t\right)=\frac{e^{it}}{f_{\alpha }\left(t\right)}$$

for all t ≠ 0 with

$$f_{\alpha }\left(t\right)=1+\alpha {\displaystyle {\int }_0^1\left(1-e^{itu}\right)u^{-\left(\alpha +1\right)}}du.$$

The form of the characteristic function is the same as the characteristic function for lim_n→∞T_n/Y₍₁₎, where the Y_i are iid random variables with power law tails, Y₍₁₎ = max_i≤nY_i, and $T_n={\displaystyle {\sum }_{i=1}^n{Y}_i}$ (see, for example, Darling 1952). Again, this agreement is a consequence of the previously discussed connection between Δ_n and the limiting Poisson point process.

To prove Theorem 2, we need the following notation. For a real number t, we define the function

$$sgn\left(t\right)=\{\begin{array}{l}-1, t<0\\ 0, t=0\\ 1, t>0.\end{array}$$

For 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 have

$$E\exp \left(it{T}_n/Y_{\left(1\right)}\right)\to \mathit{\psi }\left(t\right)$$

as n → ∞ whereas in the Simpson's index section, Y₍₁₎ = max_i≤nY_i and $T_n={\displaystyle {\sum }_{i=1}^n{Y}_i}.$ To conclude that T_n/Y₍₁₎ ⇒ 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

$$f_{\alpha }\left(t\right)=1+\alpha {\displaystyle {\int }_0^1\left(1-e^{itu}\right)u^{-\left(\alpha +1\right)}}du=1+\alpha {\left|t\right|}^{\alpha }{\displaystyle {\int }_0^{\left|t\right|}\left(1-e^{iv \mathrm{sgn}\left(t\right)}\right)v^{-\left(\alpha +1\right)}}dv.$$

(A7)

Since 1 − exp(iv) ∼ −iv as v → 0, the integral on the right-hand side of (A7) is finite and hence

$$\mathit{\psi }\left(t\right)=e^{it}{f}_{\alpha }^{-1}\left(t\right)\to 1$$

as t → 0. Since T_n/Y₍₁₎ ⇒ V, the fact that S_n/X₁ ⇒ 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 $V_n^{-1}$blows up near 1. As an explanation for this, we note that with probability

$$\exp \left(-\left(1-x^{-\alpha }\right)\right)\exp \left(-x^{-\alpha }\right)x^{-\alpha }=e^{-1}{x}^{-\alpha }$$

there is a point in the process bigger than x and no points in [1, x). When this happens,

$$V_n^{-1}=\frac{S_n}{X_1}\le 1+\frac{n}{x}$$

and so

$$F_n\left(y\right)=P\left(V_n^{-1}\le y\right)\ge e^{-1}{n}^{-\alpha }{\left(y-1\right)}^{\alpha }.$$

If 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

$$f\left(y\right)=\underset{M\to \mathrm{\infty }}{\lim }{\displaystyle {\int }_{-M}^M\frac{e^{it\left(1-y\right)}}{f_{\alpha }\left(t\right)}}dt.$$

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 that

$$f\left(y\right)={\left(y-1\right)}^{\alpha -1}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-it}}{{\left(y-1\right)}^{\alpha }+{{\displaystyle \int }}_0^{{\left(y-1\right)}^{-1}}\left(\right(1-e^{iut})/u^{\alpha +1)})du}},$$

thus 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 ∈ ℝ,

$$\begin{array}{c}P\left(V=a\right)=\underset{T\to \mathrm{\infty }}{\lim }\frac{1}{2T}{\displaystyle {\int }_{-T}^T{e}^{-iat}\mathit{\psi }\left(t\right)}dt\\ =\underset{T\to \mathrm{\infty }}{\lim }\frac{1}{2T}{\displaystyle {\int }_{-T}^T\frac{e^{it\left(1-a\right)}}{f_{\alpha }\left(t\right)}}dt.\end{array}$$

If we focus on the positive axis and use the change of variable s = t/T,

$$\frac{1}{2T}{\displaystyle {\int }_0^T\frac{e^{it\left(1-a\right)}}{f_{\alpha }\left(t\right)}} dt=\frac{1}{2}{\displaystyle {\int }_0^1\frac{e^{isT\left(1-a\right)}}{f_{\alpha }\left(sT\right)}}ds.$$

From display (A7) it follows that for every s ∈ (0, 1) we have e^isT(1−a)/f_α(sT) → 0 as T→∞. Note that

$$Re\left[f_{\alpha }\left(t\right)\right]=1+\alpha {\displaystyle {\int }_0^1\frac{1-\cos ut}{u^{\alpha +1}}}du>1,$$

which 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 that

$$\underset{T\to \mathrm{\infty }}{\lim }\frac{1}{2}{\displaystyle {\int }_0^1\frac{e^{isT\left(1-a\right)}}{f_{\alpha }\left(sT\right)}}ds=0.$$

A similar result holds for the integral on the negative axis and we conclude that

$$P\left(V=a\right)=0.$$

We can therefore conclude for x > 1 and h > 0 via the inversion formula (see Durrett 2005, Equation 3.2) and Fubini's theorem that

$$\begin{array}{c}P\left(V\in \left(x,x+h\right)\right)=\underset{T\to \mathrm{\infty }}{\lim }\frac{1}{2\pi }{\displaystyle {\int }_{-T}^T{\displaystyle {\int }_x^{x+h}{e}^{-ity}\mathit{\psi }\left(t\right)}}dydt\\ =\underset{T\to \mathrm{\infty }}{\lim }\frac{1}{2\pi }{\displaystyle {\int }_x^{x+h}{\displaystyle {\int }_{-T}^T{e}^{-ity}\mathit{\psi }\left(t\right)}}dtdy.\end{array}$$

Therefore, to establish the result we need to show that

$$\underset{T\to \mathrm{\infty }}{\lim }\frac{1}{2\pi }{\displaystyle {\int }_x^{x+h}{\displaystyle {\int }_{-T}^T{e}^{-ity}\mathit{\psi }\left(t\right)}}dtdy=\frac{1}{2\pi }{\displaystyle {\int }_x^{x+h}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-ity}\mathit{\psi }\left(t\right)}}dtdy.$$

This follows if we show that ${\lim }_{T\to \mathrm{\infty }}{\displaystyle {\int }_{-T}^T{e}^{-ity}\mathit{\psi }\left(t\right)}dt$ is a convergent integral and that a bounded function h exists defined on (x, x + h) such that

$$\left|h_T\left(y\right)\right|=\left|{\displaystyle {\int }_{-T}^T{e}^{-ity}\mathit{\psi }\left(t\right)}dt\right|\le h\left(y\right).$$

We first use integration by parts to see

$$h_T\left(y\right)={\displaystyle {\int }_{-T}^T\frac{e^{it\left(1-y\right)}}{f_{\alpha }\left(t\right)}}dt=\frac{i}{1-y}\left(\frac{e^{iT\left(1-y\right)}}{f_{\alpha }\left(T\right)}-\frac{e^{-iT\left(1-y\right)}}{f_{\alpha }\left(-T\right)}+{\displaystyle {\int }_{\mathit{-T}}^T\frac{e^{it\left(1-y\right)}{f}_{\alpha }^{\prime }\left(t\right)}{f_{\alpha }{\left(t\right)}^2}}dt\right).$$

Recalling that |f_α(T)| → ∞ as T → ±∞, it follows that if we establish that $f_{\alpha }^{\prime }\left(t\right)/f_{\alpha }{\left(t\right)}^2$ 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_α,

$${f^{\prime }}_{\alpha }\left(t\right)=-i\alpha t^{\alpha -1}{\displaystyle {\int }_0^t\frac{e^{iv}}{v^{\alpha }}}dv,$$

which follows by passing the derivative inside the integral in the definition of f_α. We can establish that

$$\underset{T<\mathrm{\infty }}{\sup }\left|{\displaystyle {\int }_0^T\frac{e^{iv}}{v^{\alpha }}}dv\right|<\mathrm{\infty }$$

by observing

$${\displaystyle {\int }_0^{\infty }\frac{e^{iv}}{v^{\alpha }}}dv=e^{-i\pi \left(1-\alpha \right)/2}\text{Γ}\left(1-\alpha \right)$$,

which can be found in many places, e.g., Loya (2005). Thus,

$$\left|f_{\alpha }^{\prime }\left(t\right)\right|\le \alpha t^{\alpha -1}\underset{T<\mathrm{\infty }}{\sup }\left|{\displaystyle {\int }_0^T\frac{e^{iv}}{v^{\alpha }}}dv\right|\le C_0{t}^{\alpha -1}.$$

We can similarly establish that for t sufficiently large

$${\left|f_{\alpha }\left(t\right)\right|}^2\ge C_1{t}^{2\alpha },$$

for a positive finite constant C₁. Thus for t sufficiently large

$$\left|\frac{e^{it\left(1-y\right)}{f}_{\alpha }^{\prime }\left(t\right)}{f_{\alpha }{\left(t\right)}^2}\right|\le \frac{C}{t^{\alpha +1}},$$

establishing 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 that

$$1+f_{\alpha }\left(t\right)=1-\stackrel{\mathrm{\infty }}{{\displaystyle \sum }_{n=1}}\frac{\alpha {\left(it\right)}^n}{\left(n-\alpha \right)n\text{!}}$$

and therefore

$$f_{\alpha }^{\left(k\right)}\left(t\right)=\stackrel{\mathrm{\infty }}{{\displaystyle \sum }_{n=k}}\frac{\alpha i^n{t}^{n-k}}{\left(n-\alpha \right)\left(n-k\right)\text{!}}$$

so that in particular,

$$f_{\alpha }^{\left(k\right)}\left(0\right)=\frac{i^k\alpha }{k-\alpha }$$

for all k ≥ 1. Let S(t) = logψ(t) = it − log f_α(t). Then dropping the α subscript on f_α, we have

$$S\prime \left(t\right)=\left(i-\frac{f\prime \left(t\right)}{f\left(t\right)}\right)=\left(i-\left(\log f\left(t\right)\right)\prime \right),$$

which yields the desired result for the mean:

$$EY=iS\prime \left(0\right)=i\left(i-f\prime \left(0\right)\right)=\frac{1}{1-\alpha }.$$

Now

$$S\prime \prime \left(t\right)=-\left(\log f\left(t\right)\right)\prime \prime =-\frac{f^{\prime \prime }\left(t\right)f\left(t\right)-{\left(f\prime \left(t\right)\right)}^2}{f^2\left(t\right)},$$

so

$$\text{var}\left(Y\right)=S\prime \prime \left(0\right)=-f\prime \prime \left(0\right)+{\left(f\prime \left(0\right)\right)}^2=\frac{\alpha }{2-\alpha }+\frac{{\alpha }^2}{{\left(1-\alpha \right)}^2}=\frac{\alpha }{{\left(1-\alpha \right)}^2\left(2-\alpha \right)},$$

completing the proof. ■