The pace of evolution across fitness valleys

Chaitanya S. Gokhale 1, Yoh Iwasa 2, Martin A. Nowak 3, and Arne Traulsen 1

1 Max-Planck-Institute for Evolutionary Biology, August-Thienemann-Str. 2,24306 Plon, Germany, 2 Department of Biology,

Faculty of Sciences, Kyushu University, Fukuoka 812-8581, Japan,3 Program for Evolutionary Dynamics, Department of Mathematics,

Department of Organismic and Evolutionary Biology,Harvard University, Cambridge MA 02138, USA

How fast does a population evolve from one fitness peak to another? We study the dynamics ofevolving, asexually reproducing populations in which a certain number of mutations jointly confera fitness advantage. We consider the time until a population has evolved from one fitness peak toanother one with a higher fitness. The order of mutations can either be fixed or random. If theorder of mutations is fixed, then the population follows a metaphorical ridge, a single path. If theorder of mutations is arbitrary, then there are many ways to evolve to the higher fitness state. Weaddress the time required for fixation in such scenarios and study how it is affected by the order ofmutations, the population size, the fitness values and the mutation rate.

I. INTRODUCTION

Evolutionary dynamics is based on natural selection,mutation and genetic drift [1]. It can be illustrated as thedynamics of a population in an abstract, typically high-dimensional fitness landscape. Since individuals withhigher fitness produce more offspring, the average den-sity of individuals is highest close to the fitness maxima.Many such features as the stationary population den-sity in the fitness landscape or the mutation rate underwhich a population can still be concentrated around a fit-ness maximum have been addressed [2–5]. An importantquestion is how a population evolves towards a fitnesspeak via several intermediate states. If the intermedi-ate states have the same fitness as the initial state, thenevolution to higher fitness states is neutral at first andthus poses no significant problems [6]. If the intermedi-ate states have lower fitness than the initial state, then afitness valley has to be overcome and it is more difficultto reach the fitness peak [7, 8]. In this case, populationstuck on a local peak cannot escape by natural selectionalone, because there is no evolutionary trajectory withsuccessively advantageous mutations. Instead, neutralgenetic drift becomes important.

Here, we consider the dynamics of these systems froma different perspective. We address the average time apopulation needs to transfer from one peak to anotherone. For small mutation rates and finite populations, wecalculate this average time analytically. When mutationrates are high, we can describe the system by a set ofdifferential equations and obtain the relevant times froma numerical integration of the differential equations. Inthis framework, the relevant question is how fast a pop-ulation evolves [9].

In particular, we can address the question whether apopulation evolves faster from one peak to another via dmutations if

(i) mutations have to occur in a certain order, i.e. onlya single evolutionary trajectory is available or

(ii) the order of the mutations does not matter, i.e.there are d! evolutionary paths.

In the simplest case the intermediate fitness values areidentical in both the cases and equal to that of the initialstate. Thus the only difference remaining is the numberof available paths. When the order of mutations is notfixed then multiple paths are available and the evolution-ary dynamics will be faster when compared to a singlepath. We can then ask the question: Does a populationevolve faster on a narrow ridge or a broad valley? Thisimplies that we move away from the simplest case men-tioned above and decrease the fitness in the intermediatestates of the multiple paths compared to the fitness inthe intermediate states of the single path. We show howthe pace of evolution depends on the depth of the val-ley, the number of intermediate states and the size of thepopulation.

In general, evolutionary dynamics depends crucially onthe size of the population. In a small population a singlemutation will typically reach fixation or extinction be-fore another mutation can arise. The population movesas a whole step by step on the fitness landscape. Forlarge populations, even for small mutation rates usuallymultiple types arise at the same time. This results in anon-zero population density in many states at the sametime. For intermediate mutation rates, the populationcan either move stepwise across the fitness landscape ormove several stpdf without getting concentrated in oneof the intermediate states. This phenomenon has beentermed stochastic tunneling [10]. If the mutation ratesare too small, tunneling does not occur because it is un-likely that a second mutation arises before the first onehas reached fixation or has gone extinct. If the mutationrates are high, tunneling occurs trivially, because the sys-tem can be approximated by differential equations for thedensities in the different states. These different scenariosincluding the limiting cases of stepwise evolution (typi-cal for small populations) and continuous evolution (typ-ical for large populations) can also be observed when the

arX

iv:1

003.

5828

v1 [

q-bi

o.PE

] 3

0 M

ar 2

010

2

population size is kept constant, but the mutation ratesare increased. For computer simulations increasing themutation rate is more convenient than simulating hugepopulations for moderate mutation rates.

One important example for an evolutionary process inwhich the timescales are of crucial importance is the so-matic evolution of cancer [11]. Cancer progression hasbeen investigated mathematically since the 1950s [12–14]. Of special interest are the tumor suppressor genes[15, 16]. In a normal cell, there are two alleles of the tu-mor suppressor gene. The mutation in the first allele isneutral if the second wild-type allele can sufficiently per-form the function. Inactivation of both the alleles confersa selective advantage to the cell and can lead to cancerprogression. This is an example in which the order ofmutations does not matter. Many cancers also requirecertain particular mutations that initiate cancer growthand pave the way for the accumulation of further muta-tions [17]. Recently, it has been shown that after cancerinitiation, a large number of different mutations may beinvolved in cancer progression [18–21]. So far, it is un-clear if the mutations have to occur in a specific order orif there is more variation in the order [22, 23].

For simplicity, we consider only very simple fitnesslandscapes here in which the fitness in all the intermedi-ate states is identical. In natural systems, these fitnessvalues will differ and also the mutation rate may not beconstant. In addition, sometimes the order of mutationswill matter and sometimes, it will not. Thus, sometimesa particular mutation will be a prerequisite to obtain anew function, but sometimes new mutations do not re-quire any prerequisites. For example, this is the case inthe evolution of resistance to β lactam antibiotics studiedby [24] and [7]. However, here we focus on a very simplemodel to highlight the general aspects of the dynamicsby analytical and numerical considerations.

This paper is organized as follows. We begin with thedescription of the two ways to order the mutations, thesingle path and the hypercube. We then derive analyticalapproximations of the fixation times for small mutationrates and discuss the effect of the different parameterson the fixation times. Next, we address the dynamicsfor intermediate and high mutation rates. Finally, weexplore biological examples which can be modeled usingthis approach.

II. MODEL

To model evolutionary dynamics in a haploid popu-lation of size N , we use the Moran process [1, 25]. Ineach time step, one individual is selected at random,but proportional to fitness. It produces one offspring,which replaces a randomly chosen individual. In onegeneration, each individual reproduces on average once.During reproduction, mutations occur with probabilityµ. We are interested in the time it takes until d mu-tations reach fixation in the population, starting from a

homogeneous population in the initial state without anymutants. Moreover, we aim to explore the dynamicalfeatures of this process. We restrict ourselves to two dif-ferent cases that allow the derivation of some analyticalresults.

A. Single path

If the mutations can occur only in a particular order,we have a single evolutionary path, see Fig. 1 for an illus-tration. Individuals in the initial state have fitness r0 = 1and individuals in the final state have fitness rd > 1. Itis instructive to characterize an individual by a string ofd sites, which can either be wild-type or mutated. If theorder of mutations is fixed, then a particular mutationrequires another particular mutation as a prerequisite.For simplicity, we assume that all the d− 1 intermediatestates have the same fitness rj = s < rd (j = 1, . . . , d−1).For s < 1, the joint effect of the set of mutations make upfor the loss of fitness caused by the individually deleteri-ous mutations. This can be considered as a very specialcase of epistasis [26].

B. Hypercube

If the order of mutations does not matter, evolutionarydynamics takes place on a hypercube in d dimensions cf.Fig. 1. Thus, there are 2d different types of individuals.In the initial state, we have d possible mutations. In thenext step, d − 1 mutations are available. Consequently,we have d! possible paths to fixation. Again, we assumer0 = 1 and rd > 1. Further, all individuals with some,but not all mutations have fitness s < rd.

111

0001

100

110

rd

s s

011

0001

100

010

001

101

110

111rd

FIG. 1: The order of mutations determines the geometryfor evolutionary dynamics, shown here for d = 3 sites (e.g.genes, nucleotide sites etc.). If mutations can only occur in aparticular order, only a single path is available (left). If theorder of mutations is arbitrary, evolutionary dynamics occurson a hypercube (right). The initial states have fitness 1 andthe final states fitness r > 1. All intermediate states areassumed to have the same fitness s < r. States are labeled bybit-strings, 0 is an wild-type site and 1 is a mutated site.

3

If the mutation probability is small, we do not needto make specific assumptions on the mutation process.But when the mutation probability increases, we can nolonger be certain that only a single mutation occurs dur-ing reproduction. For simplicity, we do not consider thepossibility of backward mutations. Although back muta-tions are often relevant, especially to escape from evolu-tionary dead ends [27], it is not straightforward to definethe speed of evolution in a system with backward muta-tions. This is due to the fact that for sufficiently highmutation rates, fixation in the final state might neveroccur. Other definitions of the end state of the systembecome arbitrary to a certain extend. The probabilityum→m+k that the offspring of an individual with m mu-tations has m+ k mutations (m ≤ m+ k ≤ d) is

um→m+k =

(d−mk

)µk(1− µ)d−m−k. (1)

This equation is valid for the hypercube, where the orderof mutations does not matter. Here,

(d−mk

)is the number

of different types of mutants with k additional mutations,µk is the probability that mutations occur at k sites and(1− µ)d−m−k is the probability that no mutation occursat the remaining d − m − k sites. For the single path,there is only one possibility to arrange the m+ k muta-tions. Thus, for k > 0, um→m+k is identical to Eq. (1),except that the binomial factor has to be dropped. Theprobability um→m that no mutation occurs follows from

normalization, um→m = 1 −∑d−mk=1 um→m+k. Our ana-

lytical calculations for small mutation rates as well as theconsiderations for high mutation rates are independent ofthe precise form of the mutation rates. However, we needto specify the form of the mutation probabilities to per-form our numerical simulations for intermediate and highmutation rates.

III. SMALL MUTATION RATES

A. The pace of evolution for small mutation rates

For small mutation probabilities, double mutations canbe neglected. Since mutations occur rarely, we can cal-culate the average time until d mutations are fixed in thepopulation analytically. Let us first address the evolu-tionary dynamics when mutants with fitness rm are al-ready present in a resident population of fitness rw, butno new mutations occur. This scenario is relevant whenmutation rates are sufficiently small. The probability toincrease the number of mutants from j to j + 1 is

T+j =

rmj

rmj + rw(N − j)N − jN

. (2)

Similarly, the number of mutants decreases from j to j−1with probability

T−j =rw(N − j)

rmj + rw(N − j)j

N. (3)

The probability that k mutants take over the entire pop-ulation is given by [1, 28–30]

φk

(rmrw

)=

1 +∑k−1i=1

∏ij=1

T−j

T+j

1 +∑N−1i=1

∏ij=1

T−j

T+j

=1−

(rwrm

)k1−

(rwrm

)N . (4)

If a mutant reaches fixation, the average number of gen-erations this process takes is given by [31, 32]

τfix

(rmrw

)=

1

N

N−1∑k=1

k∑l=1

φl

T+l

k∏m=l+1

T−mT+m. (5)

For a neutral process with rm = rw, this reduces toτfix = N−1. For sufficiently largeN , this is the maximumconditional fixation time of a mutant. Even for disadvan-tageous mutants (rm < rw) the conditional fixation timeis smaller than N −1 [32]. Since there are µN mutationsper generation, the time between two mutations is 1

µN .

Thus, for µ� N−2 a mutant reaches fixation before thenext one arises and mutations will not occur when a mu-tant is already present. Thus the population evolves bya process where the mutations occur one after the other,which has been termed periodic selection [33] and theo-retically described as the strong-selection weak-mutationregime [34, 35].

The total time τ until a mutation reaches fixation ina population is the sum of the waiting time until a suc-cessful mutant occurs and the fixation time of the mutantτ = τwait+τfix. The waiting time is the inverse of the mu-tation rate divided by the probability that a particularmutant is successful,

τwait

(rmrw

)=

1

µN

1

φ1

(rmrw

) . (6)

For µ→ 0, we have τwait →∞, but τfix remains approx-imately constant. Thus, τ ≈ τwait for small mutationrates. In principle, we could calculate τfix in the presenceof mutations. But since our approximation is only validfor small mutation rates, this will be a minor correction.

For µ� N−2, the population is homogeneous most ofthe time. Only occasionally, a mutant arises and reachesfixation or goes to extinction. The total time until dmutations are fixed in the population is the sum of thewaiting times for the successful mutants plus the time ofthe d fixation events. For a single path with initial fitness1, intermediate fitness s and final fitness r, we find forthe total time τS

τS = τwait (s) + (d− 2)τwait (1) + τwait (r/s) (7)

+ τfix (s) + (d− 2)τfix (1) + τfix (r/s) .

For small µ, we have τfix � τwait and hence the totaltime can be approximated by

τS =1

µ

[1

Nφ1(s)+ d− 2 +

1

Nφ1(r/s)

](8)

4

Consider now a “fitness valley”, in which the intermedi-ate states have fitness s < 1, but the final state has fitnessr > 1. To move from the fitness peak in the initial stateto the fitness state in the final state, first a disadvanta-geous mutation has to be fixed in the population. Sinceφ1(s < 1) � 1

N , the waiting time of such an event isvery long. The waiting time for the neutral mutations,τwait (1) = 1

µ and the waiting time for a successful mu-

tation, τwait (r/s) are significantly shorter. Thus, τS isdominated by 1

µNφ1(s) for s < 1 and sufficiently high N in

a fitness valley. Fig. 2 shows a good agreement betweenexact numerical simulations and our analytical approxi-mation for small mutation rates Eq. (8).

If the order of mutations is arbitrary, evolutionary dy-namics occurs on a hypercube. In this case, the wholeprocess will be faster, as we have d! possible paths insteadof a single one. Now, the waiting times in the differentstates depend on the number of mutations that are stillavailable. For the total time τH , we obtain,

τH =1

dτwait (s) +

d−2∑k=1

1

d− kτwait (1) + τwait

(rs

)(9)

+ τfix (s) + (d− 2)τfix (1) + τfix

(rs

)Note that the time of the fixations alone is identical forthe hypercube and the single path. Neglecting these fix-ation times for small µ ( as τfix � τwait) yields

τH =1

µ

[1

dNφ1(s)+

d−2∑k=1

1

d− k+

1

Nφ1( rs )

]. (10)

Since 1dNφ1(s) <

1Nφ1(s) and

∑d−2k=1

1d−k <

∑d−2k=1 1 = d−2

it is obvious that τH < τS , i.e. evolutionary dynamics isfaster if the order of mutations is arbitrary. For fitnessvalleys with s < 1 and a large population size, τH is dom-inated by 1

dµNφ1(s) . As d more mutations are available,

this is always faster than the corresponding equation fora single path, see Fig. 2.

B. Thresholds of the waiting times

Next, we derive expressions for some interestingthresholds of the waiting times in the limit of small muta-tion rates. Since evolutionary dynamics is always fasterif many paths are available, we now compare a fitnessvalley in which many paths are available to a single pathin which the order of mutations is important, but fitnessdoes not decrease in the course of evolution. The basicquestion we address here is, whether it is faster to cross abroad valley or a narrow ridge in fitness space. In otherwords, we compare τS(s = 1) to τH(s < 1). Since weconsider only small mutation rates µ, we neglect the fix-ation times τfix here, although they will not be identicalin the two scenarios. For s = 1, the single path is neutral.

105

106

0.9 1.0 1.1

Tim

e in

gen

erat

ions

(τ)

Intermediate fitness (s)

Single pathHypercube

FIG. 2: Fixation time for a single path (squares) and a hy-percube (circles) with small mutation rates (µ � N−2) fordifferent intermediate fitness values. Evolutionary dynamicsis always faster in the hypercube. Solid lines show the analyti-cal approximation for small mutation rates, Eqs. (8) and (10).Numerical simulations shown by symbols agree well with theanalytical approximation (population size N = 100, mutationrate µ = 10−5, d = 5, r = 1.1, simulations averaged over a1000 realizations).

We decrease s in the hypercube until we have identicalwaiting times. This yields an implicit expression for s,

d− 1 +1

N

1− 1rN

1− 1r

=1

dN

1− 1sN

1− 1s

+

d−2∑k=1

1

d− k+

1

N

1− ( sr )N

1− sr

(11)

From this equation, we can numerically determine s forany given N . For large N , Eq. (11) simplifies to

d(d− 1)−d−2∑k=1

d

d− k=

1

N

1− 1sN

1− 1s

≈ eN(1−s)−1

N(1− s), (12)

where we used (1−x/N)−N → ex for large N . Thus, thequantity N(1− s) becomes constant for large N , see Fig.3. Thus, we can now say how broad and deep a fitnessvalley has to be to lead to the same cumulative waitingtime as a single neutral path.

Next, we address the effect of the intermediate fitnesss, which has an important influence on the cumulativewaiting time τ . Fig. 2 shows how the waiting time de-creases with increasing fitness in the intermediate statess. If s comes very close to the fitness in the final state,the waiting time increases. This increase is seen both inthe single path and the hypercube. An increase in theintermediate state fitness will not always lead to a reduc-tion in waiting times. Instead, the fixation times reacha minimum when the fitness growth is constant betweenany two consecutive states [9, 36]. For the hypercube,the fastest trajectory will be steeper than on a singlepath: at first, many mutations are available and a big

5

0.0

1.0

2.0

3.0

4.0

5.0

6.0

100 101 102 103 104

N(1

-s)

Population size (N)

d=10d=5d=2

FIG. 3: The figure shows the threshold values for which evo-lution on a hypercube with fitness s < 1 in the intermediatestates proceeds as fast as on a single neutral path with s = 1.Full lines show N(1 − s) based on the numerical solution ofEq. 11. Above the lines, evolution proceeds faster on thehypercube. Below them, the neutral single path is faster.For N → ∞, the lines converge to a constant, see Eq. (12).The symbols show the results from numerical simulations forN = 5 and d = 10 (triangles), N = 20 and d = 5 (circles)as well as N = 80 and d = 2 (squares). Symbols are open(red) when the single path is faster and filled (blue) when thehypercube is faster (µ = 10−5, r = 1.1, simulations averagedover a 1000 realizations).

fitness increase is not necessary. Later, fewer mutationsare available and thus, the fitness should increase faster.The precise form of the trajectory will in this case dependon the number of mutations d and the population size N .We note that a similar reasoning can be applied to con-struct a fitness landscape that allows to cross a fitnessvalley fastest: the fastest trajectory has the same formregardless if a fitness peak is approached (r > 1) or a fit-ness minimum is approached (r < 1). Thus, the fastestway to cross a fitness valley is to descend to the minimumwith exponentially decreasing fitness and to increase fromthe minimum again with exponentially increasing fitness.

Now, we turn to the effect of the intermediate fitness son the individual waiting times. Eqs. (8) and (10) bothconsist of three terms each. The first term denotes thetime required to leave the initial state. The second termis the time spent in moving through all the intermediatestates. This second term is independent of s, because thetransitions are neutral. The last term denotes the timerequired to reach the ultimate state from the penultimatestate. For small values of s, the probability to fixate thedisadvantageous mutation is very small. Thus, the totaltime is dominated by the first term. When s is increasedto a threshold value s1, then the time for leaving the firststate is identical to the waiting time in the intermediatestates. For the hypercube, s1 is given by 1

dτwait (s1) =

∑d−2k=1

1d−k τwait (1), which reduces to

1− 1sN1

1− 1s1

= dN

d−2∑k=1

1

d− k. (13)

This equation can be solved numerically for specific val-ues of N and d. For the single path, the right hand side ofthis equation has to be replaced by N(d−2). For s > s1,the time to cross the intermediate states is larger thanthe waiting time in the first state. On the hypercube,we can define a second threshold for which the waitingtime in the first state is the same as the time requiredto reach the final state from the penultimate state. Thisarises because the effective mutation rate in state 0 is dtimes larger than the effective mutation rate in state d−1.

The threshold s2 is given by 1dτwait (s2) = τwait

(rs2

)or

1

d

1− 1sN2

1− 1s12

=1−

(s2r

)N1− s2

r

(14)

Again, s2 has to be determined numerically. For a singlepath, the factor d−1 in Eq. (14) has to be dropped. Thus,the threshold s2 occurs for s > 1 and is simply given bys2 =

√r.

The fixation time is also strongly influenced by thenumber of mutations d. A larger d increases the lengthof the path and usually also the fixation times. For thesingle path, this increase results only from the increasein the time required to cross the intermediate states, be-cause the time for leaving the initial state and the timeto reach the final state from the penultimate state are in-dependent of d. The time required to reach the ultimatestate from the penultimate state is also independent ofd for the hypercube, but the time required to leave theinitial state decreases with increasing d. This is becauseas d increases, there are more states available in the firsterror class and thus the effective rate of mutation out ofthe initial state increases. As for the single path, the timeto cross the intermediate states increases with d in thehypercube. For the hypercube, this interplay can lead toa non-monotonic dependence of the fixation time on d.For example, for N = 100 and s = 0.95, the fixation timeτH decreases with d for d < 31, but it increases with d ford > 31. In contrast, the fixation time always increasesmonotonically with d for the single path.

Increasing the fitness of the final state r increases theadvantage of the final state over the intermediate states.This will result in the decrease in the time required forthe population to make the last move. Increasing r hasno effect on the time required to cross the intermediatestates or the time required to move away from the initialstate. As a result, those two times remain constant evenas r increases, both in the single path and the hypercube.

6

0.0

0.2

0.4

0.6

0.8

1.0

10-3 10-2

Pro

babi

lity

of tu

nnel

ing

Probability of mutation (µ)

Single path s=1.0Single path s=0.9Hypercube s=1.0Hypercube s=0.9

FIG. 4: The probability of tunneling across the hypercube(circles, blue) is larger than in the single path (squares, red)due to the higher effective mutation rate. The tunnelingacross the valley denoted by the filled symbols (s = 0.9) isalways larger than the probability of tunneling across a flatfitness landscape denoted by open symbols (s = 1.0). Thisarises from the fact that the stepwise accumulation of muta-tions would involve the fixation of a disadvantageous mutationin the first step. All symbols show the probabilities that thepopulation tunnels at least across one state for the single pathor at least one error class for the hypercube (N = 100, d = 5,r = 1.1, averaged over a 1000 realizations).

IV. INTERMEDIATE MUTATION RATES

The analytical approach is only valid as long as themutation rate is small, µ � N−2. For higher mutationrates, the population does not have to consist of at mosttwo different types at any time. Instead, d mutations canbe fixed in the population without sequentially fixing oneafter the other. This process has been termed stochastictunneling and is of great importance in the context ofcancer initiation [10, 16, 22, 37]. Tunneling across fit-ness valleys is more likely than tunneling across a flatfitness landscape (see Fig. 4). Even for d = 2, the evolu-tionary dynamics is characterized by a doubly stochasticprocess, which makes analytical approaches tedious [10].As discussed above, for µN2 � 1 the population usuallycontains at most two different types. In this case, theprobability of stochastic tunneling will be very small. Onthe other hand, for µN > 1, at least one mutant is pro-duced per generation. Thus, the probability of stochastictunneling approaches 1. For N−2 < µ < N−1, the mu-tations are sometimes fixed sequentially and sometimesvia stochastic tunneling. Fig. 5 shows how the tunnelingprobability increases from 0 to 1 in this interval.

For intermediate mutation rates, it is likely that thepopulation contains more than two different types. Thetypes with beneficial mutations will compete for fixation.This process is termed clonal interference [30, 38–41].Clonal interference has been considered to slow down

adaptation, but recently it has been shown that it canhave a positive influence on a rugged fitness landscape[40, 42, 43].

The states in a single path can be characterized by thenumber of mutations. In the hypercube, the states arecharacterized by the types of mutations that have alreadyoccurred. Thus, there are many different types that haveundergone a specific number of mutations. However, alltypes that have already accumulated k mutations can bepooled into the error class k. The number of differenttypes in the error class k is given by

(dk

)= d!

k!(d−k)! .

In a single path, a population is said to tunnel acrossa state if it passes through a state without ever reach-ing fixation in that state. Analogously, in a hypercube apopulation said to tunnel across an error class if it passesthrough that error class without ever reaching fixation init. Within the error classes, tunneling can occur acrossindividual states, but also across several states at once.This means that the whole population passes only acrossthat particular state and not across any other, withoutever reaching fixation in that particular state. Tunnelingacross an error class can also occur in a second way: thepopulation can use all of the available states in the er-ror class, but the total number of individuals in the errorclass never reaches N . Thus, the probability of tunnelingvia the individual states is always lower than the proba-bility of tunneling across the error classes. Fig. 6 showsthe relation between the different probabilities of tunnel-ing in the hypercube with respect to the rate of mutation

0.0

0.2

0.4

0.6

0.8

1.0

10-4 10-3

Pro

babi

lity

of tu

nnel

ing

Probability of mutation (µ)

Single pathHypercube

FIG. 5: The probability of tunneling across a neutral hyper-cube (circles) is always higher than the probability of tunnel-ing across a neutral single path (squares). Here, the proba-bility that the system tunnels across at least one state or oneerror class is shown for d = 5 and N = 1000. As expected,for Nµ > 1 the probability of tunneling approaches 1. Incontrast to our conservative estimate that tunneling can beneglected only as long as the mutation rate is below N−2,even for mutation rates as large as 100N−2 the probability oftunneling remains close to zero (s = 1.0, r = 1.1, averagedover a 1000 realizations).

7

µ for the special case d = 2.

0

0.2

0.4

0.6

0.8

1

10-6 10-5 10-4 10-3 10-2 10-1 100

Pro

babi

lity

of tu

nnel

ing

Probability of mutation (µ)

Single pathHypercube

FIG. 6: For d = 2, the probability of tunneling across the in-termediate state is slightly higher in the single path (squares)than in the hypercube (open circles), shown for N = 100 here.This is because the effective mutation rate into the intermedi-ate state is twice as big in the hypercube, leading to a higherprobability of fixation. Filled circles show the probability totunnel across individual states of the hypercube. For µN > 1,the system always tunnels. In the hypercube, both states areused for this. As expected, we need µ� N−2 for the tunnel-ing probability to vanish (s = 1.0, r = 1.1, averaged over a1000 realizations).

Due to higher effective rates of mutation, the prob-ability of tunneling across a hypercube is expected tobe greater than or equal to the probability of tunnelingacross a single path. However, numerical simulations re-veal that for d = 2 the probability of tunneling in a singlepath is higher than in the hypercube. This is a specialcase: for d = 2 in a hypercube, the number of states intowhich the initial state can mutate into is 2. The effectiverate of mutations is thus twice as much as in the singlepath. The number of states which can be mutated intonext is one, both in the single path and the hypercube.Thus the rate at which the individuals are pushed intothe first state is higher in hypercube than in the singlepath while the rate of individuals being pushed out isthe same. Thus there is a higher probability of reachingfixation in the first error class in a hypercube (see Fig.6). We only observe this effect for d = 2, for d > 2, theprobability of tunneling is higher in a hypercube than ina single path, as expected (see Fig. 4).

V. HIGH MUTATION RATES

For µN > 1, the stochastic features of the dynamicsbecome less important. In this case, the system can bedescribed by a set of d+1 deterministic differential equa-tions for the fraction xk(t) of the population that has

k mutations [43]. Obviously, we have∑dk=0 xk(t) = 1.

Transitions out of state 0 occur with probability T0→ =(1− x0

φ u0→0)x0, where φ = x0+(1−x0−xd) s+xd r is the

average fitness of the population. This includes all thereproductive events except for the one where a type 0 isproduced. Transitions into state 0 occur with probabilityT→0 = x0

φ u0→0(1−x0). Thus, the fraction of individuals

in the initial state follows the differential equation

x0(t) =1

N

[x0

φu0→0(1− x0)− (1− x0

φu0→0)x0

]. (15)

The probability that an offspring is of type k is given

by λk =∑kj=0

xjrjφ uj→k. The difference between the

hypercube and the single path only occurs in the quantityuj→k, which is given above for both cases. The sum inλk is over all individuals with k or less mutations and rjis the fitness of individuals with j mutations. This leadsto the differential equation for the fraction of individualswith k mutations,

xk(t) =1

N[λk(1− xk)− (1− λk)xk] , (16)

where k = 0, . . . , d. Of course, the special case k = 0recovers Eq. (15). This set of d + 1 differential equa-tions describes how the system moves from state k = 0to state k = d. In general, only a numerical solution ofthis system of equations is feasible. While this allowsus to infer details of the dynamics, our main interest isthe time required for fixation of d number of mutations.Thus, we solve the differential equation numerically us-ing a standard Runge-Kutta algorithm [44]. To find anequivalent to the fixation time in a stochastic simulation,we average between fixation (xd = 1) and the time whenthere are on average less than 1 individuals outside thefinal state (xd = 1 − 1

N ). Thus, the fixation time is thetime when the solution of the differential equation crossesxd = 1− 1

2N .Fig. 7 shows an overview of the fixation times, cover-

ing the full range of mutation rates. For small mutationrates, we have sequential fixation of mutations and thetime can be well approximated by Eqs. (8) and (10). Forhigh mutation rates, the numerical solution of Eq. (16)leads to a good approximation for the fixation times.

VI. DISCUSSION

We have determined the average time during which apopulation moves from a certain initial state to a finalstate of higher fitness. The initial and the final statesare separated by a fixed number of mutations d. Themutations jointly confer a fitness advantage to the finalmutant which can be represented by a peak in the fitnesslandscape. If the intermediate mutations need to occur ina specific order for the evolution of the final mutant thenit corresponds to the single path. Otherwise, evolutionoccurs on a hypercube and there are d! ways of reachingthe final state.

8

101

102

103

104

105

106

10-6 10-4 10-2 100

Fix

atio

n tim

e in

gen

erat

ions

(τ)

Probability of mutation (µ)

Single pathHypercube

FIG. 7: The fixation times decrease with increasing mutationrate. Fixation always occurs faster on the hypercube (circles)than in the single path (squares). For small mutation rates,mutations fixate sequentially and the fixation time can be wellapproximated by Eqs. (8) and (10). Here, the fixation timesdecrease as µ−1. For high mutation rates, the system can beapproximated by a set of deterministic differential equationsand the simulation results for the fixation times can be ap-proximated based on the numerical solution of Eq. (16). Inthis case, fixation times decrease in general slower than µ−1

with increasing mutation rate (population size N = 1000,d = 5, s = 1, r = 1.1, averages over 1000 realizations).

We have explored the simplest system in which thefitness in all intermediate states is the same. As ex-pected, the fixation times on a hypercube are shorterthan on a single path, due to the presence of multiplepaths available in a hypercube. This observation leadsto the question: for which parameters does the hyper-cube show shorter fixation times than the single path,even with an added disadvantage? The fitness in the in-termediate states was then set to lower values than theones in the single path. Up to a certain threshold valueof the fitness of the intermediate states, the hypercubeshows shorter fixation times than in the single path. Thevalue of the threshold depends on the population size,total number of required mutations and the fitness in thefinal state.

The fixation times for large populations largely dependon the fitness function and are qualitatively independentof the order of mutations. Let us first focus on a flat land-scape: when the intermediate states have a fitness equalto the fitness of the initial wild-type, then for small mu-tation rates large populations have shorter fixation timesthan small populations. This is because the neutral rateof evolution does not depend on the population size. Butthe waiting time for fixation of the last mutation becomesshorter with larger population size. For intermediate mu-tation rates, tunneling starts earlier in larger populations.This leads to a marked decrease in the fixation time withlarger population size. For high mutation rates, the timeto fixation is no longer dominated by the time for the first

mutant to reach the final state, but by the time until allindividuals are in that state. Due to this, for high mu-tation rates the time required for fixation can be shorterin smaller population as compared to larger populations.Next, we focus on fitness valley: If the fitness landscapeconsists of a valley with reduced fitness of the intermedi-ate states, small populations have an advantage for smallmutation rates, as they can easily leave the initial stateand enter the valley. But for high mutation rates, largepopulations reach fixation faster, because they can ex-plore states within and beyond the fitness valley moreeasily.

Our numerical simulations reveal that tunneling canbe neglected even when the mutation rate exceeds N−2,at least by one order of magnitude. Thus, Eqs. (8) and(10) provide good estimates for the fixation times evenin relatively large populations. Concrete values for fix-ation times are collected in Table 1. They reveal thateven in long-term studies of experimental evolution, itis difficult to observe the consecutive fixation of neutralmutants [45]. Consecutive fixation of advantageous mu-tants, however, is significantly faster. For example, whileTable 1 reveals a fixation time of ∼ 1011 generations ona single path for d = 10, s = 1 and N = 106, an optimalchoice of the intermediate fitness values [9] would lead toa fixation time of ∼ 107 generations.

N d = 3 d = 10Single Path Hypercube Single Path Hypercube

102 2.10999 0.943325 9.10999 2.03896104 2.0011 0.834433 9.0011 1.93007106 2.00001 0.833344 9.00001 1.92898

TABLE I: The time required for fixation of d mutations inunits of 1010 generations for a mutation rate of µ = 10−10

based on Eqs. (8) and (10). The intermediate mutations areneutral, s = 1. For small mutation rates, the fixation timesscale linearly with µ−1. For N →∞, the fixation time on thesingle path approaches µ−1(d − 1) and the fixation time on

the hypercube approaches µ−1 ∑d−2k=0(d− k)−1. However, the

mutation rates have to decrease with increasing N to makethe approximation for the fixation times valid. (initial fitness1.0 and final fitness r = 1.1).

While we have focused on the simplest possible sys-tem which allows analytical approximations, experimen-tal studies reveal of course a much higher complexity.[7] studied experimentally the point mutations in theβ-lactamase gene of bacteria. β lactam antibiotics arecommonly used, but the bacteria can develop resistanceto the drugs. Five point mutations in a particular alleleof the β-lactamase gene increases the resistance of thebacteria to cefotaxime by a factor of ∼ 100, 000. Theo-retically the mutations leading from the wild-type alleleto the resistant allele can occur in 5! = 120 ways. Thesecan be represented by a hypercube of d = 5. But in only18 of the 120 trajectories, the intermediate mutations areeither neutral with respect to the initial state or bene-ficial. Weinreich and colleagues have shown that these

9

have the highest probability of realization. For all bene-ficial intermediates the fastest way to reach the final statewould be when the relative fitness increase between anytwo consecutive mutations is constant [9], but usually innature several different mutations are available and thepopulation first evolves to states that provide the highestselective advantage.

In another experimental study the sequence space ofthe 5s rRNA of a marine bacterium, Vibrio proteolyticuswas explored [46]. The sequences from Vibrio proteolyti-cus and Vibrio alginolyticus differ in only four positions.All the possible intermediates were constructed by theauthors and the fitness of each was calculated [47]. Twoof the valid intermediates have a fitness lower than theinitial wild-type. We have shown how such fitness valleyscan be crossed by exploring the phenomenon of tunnelingor multiple mutations (for high mutation rates). Thus,the population does not need to move in a Wrightianfashion (the whole population moving as a whole acrossthe valley).

The theory discussed herein deals with basic evolution-ary concepts which are important to the kind of biologicalexamples described above. More complex properties ofthe experimental studies like more general cases of epis-tasis and compensatory mutations can easily be incorpo-rated, but there is a huge number of possibilities. Even

if we are only interested in the ordering of fitness values,we can have up to 2d! distinct epistatic patterns. Thus,one should rather focus on concrete systems instead. Forexample, one could simulate the dynamics in a systemwith experimentally derived fitness values and mutationrates. Not all the paths of a hypercube might be acces-sible for selection, but still some of them might proveto be significant depending upon the particular valuesof the parameters, such as fitness values and populationsize. Our goal here was to characterize the simplest fea-tures of the dynamics of a population crossing a fitnessvalley. This approach can be helpful when more realisticscenarios are addressed.

VII. ACKNOWLEDGMENT

We thank Andrew Murray for initiating this project.C.S.G. and A.T. are grateful for financial support fromthe Emmy-Noether program of the DFG. M.A.N. is sup-ported by the John Templeton Foundation, the NSF/NIHjoint program in mathematical biology (NIH GrantR01GM078986), the Bill and Melinda Gates Foundation(Grand Challenges Grant 37874), and J. Epstein.

[1] Nowak MA (2006) Evolutionary Dynamics. Harvard Uni-versity Press, Cambridge, MA.

[2] Eigen M, Schuster P (1977) The hypercycle. A principleof natural self-organization. Part A: Emergence of thehypercycle. Die Naturwiss 64:541–565.

[3] Eigen M, McCaskill J, Schuster P (1989) The molecularquasi-species. Adv Chem Phys 75:149–263.

[4] Wilke C (2005) Quasispecies theory in the context of pop-ulation genetics. BMC Evol. Biol. 5:44.

[5] Nowak MA (1992) What is a quasispecies? Trends inEcology and Evolution 7:118–121.

[6] van Nimwegen E, Crutchfield J (2000) Metastable evo-lutionary dynamics: crossing fitness barriers or escapingvia neutral paths? Bull Math Biol 62:799–848.

[7] Weinreich D, Delaney N, DePristo M, Hartl D (2006)Darwinian evolution can follow only very few mutationalpaths to fitter proteins. Science 312:111–114.

[8] Poelwijk FJ, Kiviet DJ, Weinreich DM, Tans SJ (2007)Empirical fitness landscales reveal accessible evolutionarypaths. Nature 445:383–386.

[9] Traulsen A, Iwasa Y, Nowak MA (2007) The fastest evo-lutionary trajectory. J. Theor. Biol. 249:617–623.

[10] Iwasa Y, Michor F, Nowak MA (2004) Stochastic tunnelsin evolutionary dynamics. Genetics 166:1571–1579.

[11] Frank SA (2007) Evolutionary Dynamics of Cancer.Princeton Univ. Press.

[12] Fisher J (1959) Multiple-mutation theory of carcinogen-esis. Nature 181:651–652.

[13] Nordling C (1953) A new theory on cancer inducingmechanism. Br J Cancer 7:68–72.

[14] Armitage P, Doll R (1954) The age distribution of cancer

and a multi-stage theory of carcinogenesis. Br J Cancer8:1–12.

[15] Knudson AG (1971) Mutation and cancer: Statisticalstudy of retinoblastoma. Proc. Natl. Acad. Sci. U.S.A.68:820–823.

[16] Michor F, Iwasa Y, Nowak MA (2004) Dynamics of can-cer progression. Nat. Rev. Cancer 4:197–205.

[17] Vogelstein B, Kinzler K (2004) Cancer genes and thepathways they control. Nat. Med. 10:789–799.

[18] Sjoblom T, Jones S, Wood L, Parsons D, Lin J, et al.(2006) The consensus coding sequences of human breastand colorectal cancers. Science 314:268–274.

[19] Wood LD, Parsons DW, Jones S, Lin J, Sjoblom T, et al.(2007) The genomic landscapes of human breast and col-orectal cancers. Science 318:1108–1113.

[20] Jones S, Chen WD, Parmigiani G, Diehl F, BeerenwinkelN, et al. (2008) Comparative lesion sequencing providesinsights into tumor evolution. Proc. Natl. Acad. Sci.U.S.A. 105:4283–4288.

[21] Jones S, Zhang X, Parsons DW, Lin JCH, Leary RJ,et al. (2008) Core signaling pathways in human pancre-atic cancers revealed by global genomic analyses. Science321:1801–1806.

[22] Beerenwinkel N, Antal T, Dingli D, Traulsen A, KinzlerKW, et al. (2007) Genetic progression and the waitingtime to cancer. PLoS Comput. Biol. 3:e225.

[23] Gerstung M, Beerenwinkel N (2008) Waiting time modelsof cancer progression. Mathematical Population Studies,to appear. arXiv:0807.3638

[24] Weinreich D (2005) The rank ordering of genotypic fit-ness values predicts genetic constraint on natural se-

10

lection on landscapes lacking sign epistasis. Genetics171:1397–1405.

[25] Moran PAP (1962) The Statistical Processes of Evolu-tionary Theory. Oxford: Clarendon Press, Oxford.

[26] Weinreich DM, Watson R, Chao L (2005) Perspective:sign epistasis and genetic constraint on evolutionary tra-jectories. Evolution 56:1165–1174.

[27] DePristo MA, Hartl DL, Weinreich DM (2007) Muta-tional reversions during adaptive protein evolution. Mol.Biol. Evol. 24:1608–1610.

[28] Karlin S, Taylor HMA (1975) A First Course in Stochas-tic Processes. London: Academic, 2nd edition edition.

[29] Ewens WJ (2004) Mathematical Population Genetics.Springer, New York.

[30] Crow JF, Kimura M (1970) An Introduction to Popula-tion GeneticsTheory. Harper and Row, New York, NY.

[31] Goel N, Richter-Dyn N (1974) Stochastic Models in Bi-ology. Academic Press, New York.

[32] Antal T, Scheuring I (2006) Fixation of strategies for anevolutionary game in finite populations. Bull Math Biol68:1923–1944.

[33] Atwood KC, Schneider KL, Ryan FJ (1951) Periodic se-lection in escherichia coli. Proc. Natl. Acad. Sci. U.S.A.37:146–155.

[34] Gillespie JH (1983) Some properties of finite populationsexperiencing strong selection and weak mutation. AmNat. 121:691–708.

[35] Gillespie JH (2004 (2nd edition)) Population Genetics :A Concise Guide. The Johns Hopkins University Press.

[36] Weinreich DM, Chao L (2005) Rapid evolutionary escapeby large populations from local fitness peaks is likely innature. Evolution 59:1175–1182.

[37] Nowak MA, Michor F, Komarova NL, Iwasa Y (2004)Evolutionary dynamics of tumor suppressor gene inacti-vation. Proc. Natl. Acad. Sci. U.S.A. 101:10635–10638.

[38] Fisher RA (1930) The Genetical Theory of Natural Se-lection. Clarendon Press, Oxford.

[39] Muller HJ (1932) Some genetic aspects of sex. Am Nat66:118–138.

[40] Gerrish PJ, Lenski RE (1998) The fate of competingbeneficial mutations in a asexual population. Genetica102/103:127–144.

[41] Park SC, Krug J (2007) Clonal interference in large pop-ulations. Proc. Natl. Acad. Sci. U.S.A. 104:18135–18140.

[42] Wilke C (2004) The speed of adaptation in large asexualpopulations. Genetics 167:2045–2053.

[43] Jain K, Krug J (2007) Deterministic and stochasticregimes of asexual evolution on rugged fitness landscapes.Genetics 175:1275–1288.

[44] Press WH, Teukolsky SA, Vetterling WT, Flannery BP(2007) Numerical Recipes 3rd Edition: The Art of Scien-tific Computing. New York: Cambridge University Press.

[45] Cooper T, Rozen D, Lenski RE (2003) Parallel changesin gene expression after 20,000 generations of evolution inescherichia coli. Proc. Natl. Acad. Sci. U.S.A. 100:1072–1077.

[46] Lee TH, DSouza LM, Fox GE (1997) Equally parsimo-nious pathways through an RNA sequence space are notequally likely. J. Mol. Evol. 45:278–284.

[47] Chao L, McBroom SM (1985) Evolution of transpos-able elements: an IS10 insertion increases fitness in Es-cherichia coli. Molecular Biology and Evolution 2:359–369.