Traveling waves of selective sweeps

Traveling waves of selective sweeps

\fnmsRick \snmDurrett\thanksreft1label=e1]rtd@math.duke.edu [    \fnmsJohn \snmMayberry\corref\thanksreft2label=e2]jmayberry@pacific.edu [ Department of Mathematics
Duke University
Box 90320
Durham, North Carolina 27708-0320
USA
Department of Mathematics
University of the Pacific
3601 Pacific Avenue
Stockton, California 95211
USA
Duke University and University of the Pacific
\smonth10 \syear2009\smonth6 \syear2010
\smonth10 \syear2009\smonth6 \syear2010
\smonth10 \syear2009\smonth6 \syear2010
Abstract

The goal of cancer genome sequencing projects is to determine the genetic alterations that cause common cancers. Many malignancies arise during the clonal expansion of a benign tumor which motivates the study of recurrent selective sweeps in an exponentially growing population. To better understand this process, Beerenwinkel et al. [PLoS Comput. Biol. 3 (2007) 2239–2246] consider a Wright–Fisher model in which cells from an exponentially growing population accumulate advantageous mutations. Simulations show a traveling wave in which the time of the first -fold mutant, , is approximately linear in and heuristics are used to obtain formulas for . Here, we consider the analogous problem for the Moran model and prove that as the mutation rate , , where the can be computed explicitly. In addition, we derive a limiting result on a log scale for the size of the number of cells with mutations at time .

[
\kwd
\doi

10.1214/10-AAP721 \volume21 \issue2 2011 \firstpage699 \lastpage744 \newproclaimAssumptionsAssumptions

\runtitle

Traveling waves of selective sweeps

{aug}

and

\thankstext

t1Supported in part by NSF Grant DMS-07-04996 from the probability program. \thankstextt2Supported in part by NSF RTG Grant DMS-07-39164.

class=AMS] \kwd[Primary ]60J85, 92D25 \kwd[; secondary ]92C50.

Moran model \kwdselective sweep \kwdrate of adaptation \kwdstochastic tunneling \kwdbranching processes \kwdcancer models.

1 Introduction

Recent studies have sought to identify the mutations that give rise to common cancers by sequencing protein-coding genes in common tumor types, including breast and colon cancer SJW (), WPJ (), pancreatic cancer Jetal () and glioblastoma Petal (), TCGA (). The last study is part of a 100 million dollar pilot project of the NIH, which could lead to a 1.5 billion dollar effort. These studies have rediscovered genes known to play a role in cancer (e.g., APC, KRAS and TP53 in colon cancer), but they have also found that tumors contain a large number of mutations. Analysis of 13,023 genes in 11 breast and 11 colorectal cancers in Sjoblom et al. SJW () revealed that individual tumors accumulate an average of 90 mutated genes, but only a subset of these contribute to the development of cancer.

Follow-up work in Wood et al. WPJ () studied 18,191 distinct genes in the same 22 samples. Any gene that was mutated in a tumor but not normal tissue was analyzed in 24 additional tumors, and selected genes were further analyzed in 96 colorectal cancers. Statistical analysis suggested that most of the 80 mutations in an individual tumor were harmless and that 15 were likely to be responsible for driving the initiation, progression or maintenance of the tumor. These two types of mutations are commonly referred to as “drivers” and “passengers.” The latter provide no selective advantage to the growing cancer mass, but are retained by chance during repeated rounds of cell division and clonal expansion (exponential growth).

The results of SJW () and WPJ () are in contrast with the long-held belief that most cancers are the end result of a handful of mutations. Armitage and Doll AD () constructed log–log plots of cancer mortality versus age and found slopes of 5.18 and 4.97 for colon cancer in men and women, respectively. From this they predicted that the occurrence of colon cancer was the result of a six-stage process. In essence their argument is that the density function of the sum of six exponentials with rates is

 ≈μ1⋯μ6t5/5!for small t.

This result yields the density of the well-known gamma distribution when all the are equal, but only readers with well developed skills in calculus (or complex variables) will succeed in deriving this result for unequal .

The work in AD () and the subsequent work of Knudson K0 (), who used statistics to argue that retinoblastoma was the end result of two mutations, gave rise to a large amount of work; see KNRC () and the books by Wodarz and Komarova WK () and Frank Fr () for surveys. From this large body of work on multistage carcinogenesis, we will cite only two sources. Luebeck and Moolgavakar LM () used multistage models to fit the age-specific incidence of colorectal cancers in the SEER registry (which covers 10% of the US population) to conclude that a four-stage model gives the best fit. Calabrese et al. Cetal () used data for 1022 colorectal cancers to argue that “sporadic” cancers developed after six mutations, but that in the subgroup of individuals with strong familial predispositions, only five mutations were required.

There is good reason to doubt some of the conclusions of SJW () and WPJ (). First, the statistical methods of SJW () have been criticized (see letters on pages 762–763 in the February 9, 2007 issue of Science). Furthermore, in WPJ (), a follow-up study on 40 of the 119 highest scoring genes, chosen because they were in pathways of biological interest, showed that 15 of the 40 genes (37.5%) were not mutated in any of the 96 tumors, casting doubt on the claimed 10% false discovery rate. However, the more recent studies Jetal (), Petal (), TCGA () using well-known and trusted statistical methods have found similar patterns: an average of 63 mutations in pancreatic cancers and 47 in glioblastoma.

To better understand this process by which an exponentially growing cell mass accumulates driver and passenger mutations, and, in particular, to understand the data in SJW (), Beerenwinkel et al. BN () considered a Wright–Fisher model with selection and mutation in an exponentially growing population. They assumed that there were 100 potential driver genes and asked for the waiting time until one cell has accumulated mutations. Simulations (see their Figure 3) showed that a traveling wave developed, in which the time until the first -fold mutant was approximately linear in and the authors used heuristic arguments to obtain quantitative predictions for the first time that a cell with mutations appears.

Here we will consider this problem for the analogous Moran model and prove asymptotic results as the mutation rate for the behavior of the number of cells with mutations at time . A cell with mutations will be referred to as a type- individual. Our main result is Theorem 2, which allows for an exponentially growing population of individuals. The process of fixation of advantageous mutations in a population of constant size has been the subject of much theoretical work (see, e.g., Chapter 6 of RD ()), so it is natural to ask how the behavior changes in an exponentially growing population. A second difference from the standard theory of the fixation of a single mutation is that we consider a situation in which new mutations arise before older ones have gone to fixation, a process often referred to as “stochastic tunneling.” The resulting “Hill–Robertson” interference (see, e.g., Section 8.3 in RD ()) can be analyzed here because only the newest mutation is stochastic while the older mutations behave deterministically. This idea was used by Rouzine et al. in RWC () (and later developed in more detail in BRW (), RBW ()) as a heuristic, but here it leads to rigorous results.

The rest of the paper is organized as follows. In Section 1.1 we begin with a fixed population size of individuals and state Theorem 1, which says that when time is scaled by , the sizes of , divided by , converge to a limit that is deterministic and piecewise linear and so the time the first type- individual appears is . Since we have assumed the population size is , this time scale agrees with (i) results in YE (), YEC (), which show that the rate of adaptation (defined as the change in the mean fitness of the population) for a related fixed-population-size Moran model is and (ii) simulations in DF () which suggest that the speed of adaptation depends logarithmically on both the mutation rate and the population size. In Section 1.2 we return to the growing population scenario and state our main result, Theorem 2, which generalizes Theorem 1. Section 2 contains some examples elucidating the nature of the limit in Theorem 2 and illustrating the traveling-wave-like behavior of the limit. We state the main tools used to prove Theorem 2 in Section 3, and Sections 4 and 5 contain the technical details.

1.1 Fixed population size: Main result

We begin by considering our Moran model in a fixed population of individuals and return to our analysis of the exponentially growing population in Section 1.2. We assume that:

{longlist}

[(iii)]

initially, all individuals are of type ;

type- individuals mutate to individuals of type at rate ;

all individuals die at rate 1 and, upon death, are replaced by an individual of type with probability

 (1+γ)jXμj(t)Wμ(t),

where is the relative fitness of type- individuals compared to type-0, and

 Wμ(t)=∞∑i=0(1+γ)iXμi(t)

is the “total fitness” of the population. We assume throughout that is fixed (i.e., mutations are advantageous). Approximations of the time the first type- individual appears have been carried out for the neutral case () in INM (), HIM (), DSS (), JS () (and applied to regulatory sequence evolution in DS ()). The case is of interest in studying Muller’s ratchet Mul (), but since deleterious mutations behave very differently from advantageous mutations, we will not consider this case here.

We will suppose throughout that , that is, . If , then the 1’s arise and go almost to fixation before the first mutation to a 2 occurs, so the times between fixations are independent exponentials. We will not here consider the borderline scenario, although we note that in the case when and , the limiting behavior of the system has been well studied and we obtain a diffusion limit describing the evolution of type- frequencies (see, e.g., Sections 7.2 and 8.1 in RD ()). Let and for define

 Tμj =inf{t≥0\dvtxXμj(t)≥1},

that is, is the time of the first appearance of a type- individual. In order to study the birth times we will prove a limit theorem for the sizes of the on a scale. Let , and define

 γj=(1+γ)j−1

for all . In what follows, we shall use to denote the greatest integer less than or equal to and let .

Theorem 1

Suppose that with for some and suppose that , the set of generic parameter values defined in (2). Then, as

 Yμj(t)≡1Llog+(Xj(Lt/γ))→yj(t)in probability

uniformly on compact subsets of where is defined in (1). The limit , which depends on the parameters , is deterministic and piecewise linear and will be described by (a) and (b) below. Furthermore, if we define

 tj=tj(α,γ)=inf{t\dvtxyj(t)=1}

for , then

 Tμj+1L/γ→tjin probability

as for all .

1. [(a)]

2. Initial behavior. . The convergence only occurs on because we have for all by assumption, so a discontinuity is created at time 0.

3. Inductive step. Let and suppose that at some time the following conditions are satisfied: {longlist}[(ii)]

4. and both exist and for all

5. for all so that, in particular, .

Let if , if and define

For all , we then have

and conditions (i) and (ii) are satisfied at time .

Our description of the limiting dynamical system can be understood as follows. If type is the most fit of the dominant types in the population at time , then the , , grow linearly with slope while the , , decrease linearly with slope until they hit zero and stays constant. These rates remain valid until either reaches level for some and there is a change in the most fit dominant type, or reaches level 1 and individuals of type are born. These two events correspond to and , respectively. The condition guarantees that after birth, the growth of type- individuals is driven by selection and not by mutations from type- individuals. If this condition failed, we would encounter a discontinuity in the limiting dynamics like the one at time 0. We have rescaled time by since, in most cases of interest, is small (e.g., ) and when is small, we have so that the limit process described above is almost independent of .

Parts (a) and (b) together describe the limiting dynamical system for all times

 t

since by part (a), (i) and (ii) in (b) hold at time and it is easy to see from the form of that if (i) and (ii) hold at time , then they also hold at time . Note that the form of the limit implies that at times , there is always a unique such that , that is, a unique dominant type. This observation will prove useful on a number of occasions. In Section 2 we shall see examples in which , but in Section 2.4 we will prove the following result.

Lemma 1

For any , there exists such that whenever .

Therefore, in general, our construction cannot be extended up to arbitrarily large times. We prove this lemma by showing that if is large, then an infinite number of types will be born before any existing type- individual achieves dominance. However, since it is easy to see by construction that we have for any , this is the only way that blow-up can occur. Hence, the limit process will accumulate an infinite number of mutations before time which means our approximation is valid on any time interval of practical importance.

The generic set of parameter values is

 G(α)≡{γ∈(0,∞)\dvtxδn,j≠δn,i % for all i≠j,n≥0}. (2)

In other words, when , there is always a unique such that . For any given value of , is clearly countable, so our result is good enough for applications. If for some , then either (i) type ’s and type ’s reach level at the same time or (ii) type ’s reach level at the same time that type ’s reach level 1. It is tempting to argue that since generic parameters are dense, the result for general parameters follows, but proving this is made difficult by the fact that the growth rates are not continuous functions of the parameters since they depend on .

Theorem 1 is very general, but not very transparent. In Section 2 we will give some examples in which more explicit expressions for the birth times are available. Figure 1 shows examples in the first three “regimes” of behavior that we will consider. In the th regime, type arises (but not type ) before type  “fixates,” that is, is of order . These regimes closely correspond to the different scenarios considered in BRW (), in which the “stochastic edge,” that is, the class of the most fit mutant present at positive quantities, is always assumed to be fitness classes ahead of the population bulk. is referred to as the “lead.” In the notation of Theorem 1, the lead is always on the interval and in regime , the lead is always . In all three examples in Figure 1, we see the traveling-wave-like behavior observed in the simulations of Beerenwinkel et al. BN () (see also RWC ()). The time between successive waves is constant in the example from regime 1, while in the examples from regimes 2 and 3, the time between successive waves is not constant, but converges to a constant as the number of waves goes to infinity.

1.2 Growing population

We now consider a growing population of individuals , with a random initial population in distributed according to some measure . At time , all individuals are of type 0 and we suppose that, in addition to the previously imposed Moran dynamics, at rate , , new individuals are added and their type is chosen to be with probability

 (1+γ)kXμk(t)Wμ(t).

As in the case of fixed population size, we are able to derive a limiting, piecewise linear approximation to

 Yμk(t)≡(1/L)log+Xμk(Lt/γ).

To determine the correct growth rates, suppose that there are individuals of type  and that the population size is . We then have the jump rates

 xj ↦ xj+1rate: [(1+ρ)N−xj](1+γ)jxj∑i≥0(1+γ)ixi+μxj−1, xj ↦ xj−1rate: xj∑i≠j(1+γ)ixi∑i≥0(1+γ)ixi−μxj.

If mutations can be ignored, then the growth rate of type ’s is

 ∑i≥0[(1+ρ)(1+γ)j−(1+γ)i]xixj∑i≥0(1+γ)ixi≈[(1+ρ)(1+γ)j−m−1]xj

if for (recall that in the limiting dynamical system from Theorem 1, there is a unique dominant type at time for all but a countable number of times). This yields the expression

 λj−m≡(1+ρ)(1+γ)j−m−1

for the limiting growth rate of type ’s in a population dominated by type .

If type- individuals have size at time 0, are growing at rate for some and the initial total population size is , then type ’s will achieve fixation at the approximate time satisfying

 (1/μ)xeλℓ(j)t=(1/μ)zeρtort=z−xλℓ(j)−ρlog(1/μ).

This leads to the following result. Theorem 1 corresponds to the special case .

Theorem 2

Let and suppose that in probability for some . Then, for all , the set of generic parameter values given in (3), as we have and in probability uniformly on compact subsets of and respectively, where is given in (4). The limits , which depend on the parameters , are deterministic and piecewise linear and described by (a) and (b) below. Furthermore, if we define

 tj=tj(α,ρ,γ)=inf{t\dvtxyj(t)=1}

for , then

 Tμj+1L/γ→tjin probability

as for all .

1. [(a)]

2. Initial behavior. .

3. Inductive step. Let and suppose that at some time the following conditions are satisfied: {longlist}[(ii)]

4. and both exist and for all

5. for all so that, in particular, .

Let , if , if and define

For all , we then have

and conditions (i) and (ii) are satisfied at time .

The generic set of parameter values is

 G(α,ρ)≡{γ∈(0,∞)\dvtxδn,j≠δn,i for all i≠j,n≥0}\vspace∗−2pt (3)

and, of course,

 t∗≡∞∑n=1Δn.\vspace∗−2pt (4)

The argument which we use to prove Lemma 1 implies that whenever since as .

An example is given in Figure 2. Since the population size is growing, we progress through the different “regimes” of behavior defined earlier for the fixed population size, and the time between successive waves of sweeps decreases. This behavior can also be seen in Figure 3 of BN (). Here we are dealing with the small mutation limit so that our waves have sharp peaks.

Motivated by the statistical analysis of cancer data in SJW (), Beerenwinkel et al. BN () were interested in the time at which a cell first accumulates 20 mutations. Their choice of the number 20 was inspired by data from SJW (). Using heuristics they obtained the approximation

 Tμj≈sj=j(log(γ/μ))2γlog(N(0)N(Tμ20))\vspace∗−2pt (5)

for . Their model evolves in discrete time, but the heuristics only use the fact that the drift in the Wright–Fisher diffusion limit (ignoring mutations) is given by

 bj(x)≈γxj(j−⟨j⟩),

where (see RD (), page 253). To get the same drift in continuous time, we need to rescale time by , as opposed to and hence we should replace by and by to obtain the analogous approximations for the Moran model. The important point to note is that the approximation in (5) is linear in and hence yields constant estimates for the increments , whereas we can see that in the limiting dynamical system, the increments are not constant, but decrease in length as the population size increases. Figure 3 shows a plot of , the constants from Theorem 2, as a function of and illustrates the nonlinearity in .

2 Examples

In this section, we discuss some examples with constant population size () in which the explicit computation of the times is possible. To do so, it is more convenient to study the increments

 τμj≡Tμj−Tμj−1.

Theorem 1 then implies that

 τμjL/γ→βj,

where for all if we use the convention that . We begin with the first regime of behavior where the lead is always 1 and we have for all . In Section 2.2 we move on to discuss regime 2 where the lead is always 2. In this case, we will see that depends on , but we have as so that asymptotically, the times between successive waves is constant. Section 2.3 contains some conjectures on the parameter range for the regime . In all these scenarios, we conjecture that and our limiting result holds for all time. In Section 2.4, we will prove Lemma 1, showing that for any , there exists such that whenever .

2.1 Results for regime 1

Let . The first regime occurs for . If is small, then and the condition is roughly . If , then so throughout regime 1.

Table 1 summarizes the situation. To explain the entries, we note that applying part (a) of the limit description implies that and part (b) then implies that

 y1(s)=(α−1)+s

for . Since we have assumed that , we have

 Δ1=δ1,1∧δ1,2=(α−1)∧γ2γ=α−1,

and applying part (b) tells us that we have for all . Another application of (b) then yields which gives the additional amount of time needed for to hit . Since the relative sizes of 1’s, 2’s and 3’s at time are the same as the relative sizes of 0’s, 1’s and 2’s at time , we obtain the following result giving the limiting coefficients of .

Corollary 1

Suppose that for some . Then, as

 τμ1L/γ→(2−α),and for % all j≥2,τμjL/γ→βin % probability,

where

Figure 1 illustrates the limiting dynamical system in the case where and . We can see that in regime 1, the system is characterized by a “traveling wave of selective sweeps” in type space, that is, the growth and decay of types occur translated in time by a fixed amount. In Figure 4 we show the distributions of types at the times when type-5, -9, -13 and -17 individuals are born (the times , , and from Theorem 1). As we move from time to , the distribution is shifted by a constant amount.

2.2 Results for regime 2

Regime 2 occurs for with . When is small, so this regime is roughly . In general, so we have throughout this regime. As in the previous section, it is easiest to explain the conclusions of Theorem 1 with a table; see Table 2.

Since the first two rows are the same as in regime 1, and we again have

 τμ1∼(2−α)L/γ.

However, we now have so that

 Δ1=δ1,1∧δ1,2=(α−1)∧(γ/γ2)=γ/γ2,

and hence the 2’s reach level before the 1’s fixate. This yields

 τμ2∼γγ2⋅Lγ.

Now , so the additional time it takes to reach level is . Since , we have and hence , that is, the 1’s will fixate before the 3’s reach level . To show that the 1’s fixate before the 2’s and conclude that , we need to show that which holds if and only if

 α<2+γ1+γ. (6)

However, comparing with the upper bound in (6), we can see that

 1+γ/γ2+γ/γ3<2+γ1+γ ⟺ (3+3γ+γ2)+(2+γ)(2+γ)(3+3γ+γ2)<11+γ ⟺ 5+9γ+5γ2+γ36+9γ+5γ2+γ3<1.

The last inequality is always true and therefore (6) holds throughout regime 2 and , justifying the fourth line in Table 2. Finally, to check that the 2’s have not yet fixated when the 3’s reach level and prove

 Δ3=δ3,3=γγ2(1−Δ2γ3γ),

we note that the size of is

 1+Δ2γ2/γ+δ3,3 = 1+γ2γ(α−r2)+γγ2−γ3γ2(α−r2) = 1+γ/γ2+ℓ(α−r2)

with

 ℓ≡γ2/γ−γ3/γ2=(2+γ)2−(3+3γ+γ2)2+γ=1+γ2+γ∈(0,1),

and hence . This justifies the final line of Table 2, and we conclude that

 τμ3L/γ→Δ2+Δ3=α−r2+γ(1−γ3(α−r2)/γ)γ2.

In contrast to regime 1, the relative sizes of types when the 3’s reach are not exactly the same as the relative sizes when the 2’s reach level . To describe this more complicated situation, suppose that type- individuals have size at the time type- individuals reach level . Then, if we assume:

{longlist}

[(2a)]

type reaches fixation before type ;

type reaches fixation before ’s reach ;

type reaches level before type reaches fixation, we can repeat the arithmetic leading to Table 2 to yield Table 3, where here with as before.

Since the density of 2’s is when the 3’s have reached size , we see that when type reaches size the density of type is . This leads to the statement of our next result.

Corollary 2

Suppose for some . Then, as

 Tμ1L/γ→(2−α),and for all % j≥2,τμjL/γ→βjin % probability,

where and if we let then for all we have

Furthermore, the coefficients as where

 β∞=α−r∗+1−(3+3γ+γ2)(α−r∗)2+γ

with .

{pf}

We need to show that conditions (2a)–(2c) above are satisfied for any and that converges. The latter follows from the fact that has slope with , so, as

 fj(r2)→r∗=r2+ℓα1+ℓ,

the unique fixed point of . It is easy to see that implies that

 r2≤fj(r2)<α (7)

for all and (2c) immediately follows. Since , we have , which, along with (7), tells us that (2b) holds for all as well. Finally, (2a) is equivalent to

 α−1γ2>α−fj−3(r2)γ

and so (7) implies that to prove (2a), we need only show that

 α−1γ2>α−r2γ.

Rearranging terms, we obtain the equivalent condition which holds by (6), completing the proof.

Again, the behavior of the limits can be read from Tables 2 and 3. The formulas are messy, but it is easy to compute for a fixed value of . As Figure 1 shows, after a short transient phase, the increments between the appearance of successive types settle down into the steady-state behavior guaranteed by Corollary 2. Figure 4 shows the distribution of types at various times throughout the evolution of the system, which agree with simulations given in Figure 1 in the Appendix of BN ().

2.3 Regime j, j≥3

Regime 3 occurs for with . When is small, , so this regime is roughly . If then the initial phases are similar to regime 2, but now type 3 reaches before the 1’s fixate; see Table 4.

If we now assume that: {longlist}

type reaches fixation before types and ;

type reaches fixation before type ’s reach ;

type reaches level before types and reaches fixation, then the recursion in Table 3 becomes a pair of equations (see Table 5). To imitate the proof in regime 2 we would have to show that (3a)–(3c) hold for and , and for all of the iterates where . Figure 5 shows that this is true when and ; however, verifying this algebraically is difficult because may fail to satisfy the conditions when does.

In general, we conjecture that if we define

 rj=j∑i=1(γ/γi),

then we are in regime if and we have as . In particular, this would imply that as long as

 α

However, as Figure 6 shows, the converse is not true. There, we have so that , but it appears that we still approach a constant increment between waves.

2.4 Blow-up in finite time

In this section we prove Lemma 1, which shows that for any , we can choose large enough to make finite. To this end, let and for all define

 Sj=Sj(γ)=∞∑i=0γjγj+i

and let . Note that since

 Sj=∞∑i=01−(1+γ)−j(1+γ)i−(1+γ)−j≤∞∑i=01(1+γ)i⋅11−(1+γ)−j−i,

we have . With this notation to hand, we can prove the lemma.

{pf*}

Proof of Lemma 1 Fix , define as above and choose large enough so that and

 γa/2/γa<1/S,

where . We will show that and for so that

 t∗=∞∑n=1Δn≤γ∞∑n=0[(1+γ)a+n−1]−1<∞.

To prove that has the desired bound, we will show that for all , hits level 1 before hits level for any . Since this implies that the growth rate of type ’s is , we have, in the notation of Theorem 1, for all so that and