The rate of multi-step evolution in Moran and Wright-Fisher populations
TThe rate of mutli-step evolution in Moran and Wright-Fisherpopulations
Stephen R. Proulx a, ∗ a Ecology, Evolution and Marine Biology Department, UC Santa Barbara, Santa Barbara, CA 93106,USA
Abstract
Several groups have recently modeled evolutionary transitions from an ancestral allele toa beneficial allele separated by one or more intervening mutants. The beneficial allele canbecome fixed if a succession of intermediate mutants are fixed or alternatively if successivemutants arise while the previous intermediate mutant is still segregating. This latterprocess has been termed stochastic tunneling. Previous work has focused on the Moranmodel of population genetics. I use elementary methods of analyzing stochastic processesto derive the probability of tunneling in the limit of large population size for both Moranand Wright-Fisher populations. I also show how to efficiently obtain numerical resultsfor finite populations. These results show that the probability of stochastic tunneling istwice as large under the Wright-Fisher model as it is under the Moran model.
Keywords:
Population Genetics, Stochastic Process, Stochastic Tunneling, FixationProbability
1. Introduction
Evolutionary biologists have long understood that transitions between adaptive setsof traits may involve multiple substitutions separated by neutral or maladaptive inter-mediate states (Wright, 1932). There has been a resurgence of interest in these ideas, inpart because of advances in methods to measure epistatic interactions (e.g. Tong et al.,2001, 2004) and ability to observe evolutionary trajectories (Weinreich and Chao, 2005).Several researchers have modeled evolutionary processes when epistatic interactions allowfor multiple genotypes to have the same direct effect on fitness but experience differentevolutionary dynamics because of differences in their genetic robustness(van Nimwegenet al., 1999; de Visser et al., 2003; Proulx and Phillips, 2005; Draghi et al., 2010) or thelocal mutational landscape (Wilke et al., 2001; O’Fallon et al., 2007). These scenarioscan be called circum-neutral because alternative genotypes differ in their long-term evo-lutionary dynamics only because of the genomic circumstances in which they are found(Proulx and Adler, 2010).Several groups have extended the theory to describe the rates and probability oftransition along a multi-step evolutionary pathway. Weinreich and Chao (2005) took the ∗ Corresponding Author
Preprint submitted to Theoretical Population Biology June 1, 2018 a r X i v : . [ q - b i o . P E ] J u l pproach of calculating the total waiting time along various pathways and comparing therelative waiting times to reach a final state. Hermisson and Pennings (2005) considered ascenario where previously accumulated genetic variation may become adaptive followingan environmental shift. In this scenario the population genetic dynamics of standingvariation plays an important role in determining how evolution proceeds at the next stepin the process (see also Kopp and Hermisson (2009)). Iwasa et al . (2004) derived approx-imate results on the waiting time and probability of a two-step sequence of mutationaltransitions using the Moran model, while Iwasa et al. (2003) derived results in a Wright-Fisher model for a scenario where multiple mutations are required to escape the immuneresponse. These results have been utilized by several other groups to study the rate ofmulti-step evolutionary processes (Durrett and Schmidt, 2008; Lynch, 2010; Lynch andAbegg, 2010). Several other works have explored the probability and timing of multi-step processes, as well as exploring the validity of approximations (Schweinsberg, 2008;Weissman et al., 2009; Durrett et al., 2009). Both Schweinsberg (2008) and Weissmanet al. (2009) have presented branching process approximations for large populations thatare equivalent to the large population size limit results for the Moran model presentedhere.The goal of this paper is to show how the finite population processes for both theMoran model and the Wright-Fisher model can be written and solved using the methodof first step analysis. This helps to clarify some of the terms described by Iwasa etal . (2004) and gives an algorithm for efficiently solving the finite population Moran model.The Moran tunneling probabilities have previously been applied to Wright-Fisher pop-ulations without verifying that these results still hold. I show that the Wright-Fishertunneling probabilities differ from the Moran probabilities by a factor of 2. This correc-tion will allow stochastic tunneling results to be applied to a wider range of scenarios.I also compare the large population size approximations for the rate of tunneling withsimulations and exact calculations for small population size. By considering the population level evolutionary process as a series of transitionsbetween populations fixed for a single genotype we can calculate the waiting time forthe population to become fixed for secondary mutations. So long as
N µ << r ≤
1, while the second mutational step (secondary mutant) is assumed to have fitness a > r is exactly one, the first mutationalstep has no direct effect on fitness and the primary mutants can be considered circum-neutral (Proulx and Adler, 2010). Such circum-neutral substitutions do not directly affectreproductive fitness but do alter the long-term evolutionary trajectory of the population.The ancestral population can evolve to be fixed for the secondary mutant either through2 sequential mutational pathway or because a lineage of primary mutants destined forextinction produces a secondary mutant which is destined for fixation, a process termedstochastic tunneling by (Komarova et al., 2003).The waiting time until a secondary mutation becomes fixed can be expressed in termsof the waiting times for the sequential and tunneling paths. I define the per generationprobability of successful sequential substitutions S and S and the per generation prob-ability of the opening of a successful tunnel as T . The waiting time for the transitionbetween population states is well described by an exponential waiting time so long aspopulation size is not too small (Iwasa et al., 2005). This means that the process ischaracterized by a race between waiting for a primary mutation to arise and fix and thestart of a tunneling pathway. The expectation of the total waiting time until a secondarymutation is given by E [ t ] = T ( T + S ) + S ( S + S + T ) S ( T + S ) , (1)where the first term represents the contribution to the expected waiting time from tunnel-ing pathways and the second term represents the contribution from sequential pathways.If T = 0 this is simply the sum of the waiting times for primary and secondary mutationsto sequentially fix. This approximation ignores the time that it takes for beneficial mu-tations to spread through the population and the amount of time that primary mutantsare segregating before a secondary mutation arises. The time required for alleles destinedto fix to spread to fixation is typically much smaller than the waiting times for them toarise, and in any case it can be simply added to the total waiting time (see Lynch andAbegg, 2010).The per generation probabilities of sequential fixation are S = N µ U ( r ) , (2) S = N µ U ( a/r ) , (3)where N is the haploid population size(for simplicity I assume this is approximately theeffective population size as well), µ is the probability that an ancestral allele will mutateinto a primary mutant, µ is the probability that the primary mutant will mutate intoa secondary mutant, and U ( x ) is the fixation probability of a mutation with relativefitness x when initially present as a single copy. Because this follows sequential fixationof mutants, the secondary mutant is invading into a population fixed for the primarymutant, giving it a relative fitness of a/r .Following Iwasa et al . (2004), the probability of tunneling can be written as T = (1 − U ( r ))(1 − E [no secondary substitution | extinction]) , (4)where E [no secondary substitution | extinction] represents the probability that no success-ful secondary mutations arise while the primary mutant is segregating conditioned onthe eventual extinction of the lineage of primary mutants. This can be related to theunconditional expectation by E [no secondary substitution | extinction] = E [no secondary substitution](1 − U ( r )) (5)(Iwasa et al., 2004). This provides a simple relationship between calculations made usingthe conditioned trajectory of primary mutations and the unconditioned trajectory ofprimary mutations. 3 . Moran Model The Moran model (Moran, 1962) follows a population of size N described by a vector x where x i indicates the number of individuals of genotype i . In the Moran process, eachunit of time either 2 elements of x change by one unit each in opposite directions or x remains constant. This model is often conceptually described as one where an individualis chosen to reproduce at random weighted by their reproductive output. Population sizeis kept constant by choosing one of the original population members to die (it may be theone that reproduced). Mutation causes offspring to differ from their parent’s genotypewith a probability defined by the mutation rate. The scale of time in this model ispopulation size dependent; a generation is measured in terms of N time steps.Because the number of individuals of each genotype can change by at most one uniteach time step, the Moran model can be expressed as a Markov chain whose matrixdefinition is tridiagonal. This is in sharp contrast to the Wright-Fisher model whereany population state can move to any other state in one generation (albeit with lowprobabilities). This simple matrix structure allows many features of the stochastic processto be expressed algebraically.The evolution of the population can be described as a series of transitions betweenstates described by the complement of segregating alleles. Assuming that primary mu-tations occur rarely enough so that multiple primary mutants do not typically arisetogether, the transition probabilities between states can be based on the introduction ofa single individual mutant. As this approximation breaks down more error will be in-troduced into the transition probabilities. The ancestral population is monomorphic forthe ancestral genotype. Following the introduction of a primary mutant, the populationwill evolve with two genotypes for some time until either the lineage of primary mutantsgoes extinct or a secondary mutant lineage arises and does not go extinct.Following the introduction of a primary mutant, the population is composed of i primary mutants and N − i ancestral alleles. In the absence of mutation, the Markovtransition probabilities are P r ( i → i −
1) = i ( N − i ) N ( ir + ( N − i )) P r ( i → i ) = i r + ( N − i ) N ( ir + ( N − i )) P r ( i → i + 1) = ir ( N − i ) N ( ir + ( N − i )) , where 0 < i < N and r is the relative fitness of the primary mutant. Note that thismodel ignores the change in the number of primary mutants due to their mutation intosecondary mutants. Following the introduction of the primary mutant, it may produce alineage that eventually goes extinct or eventually becomes fixed in the population. It isuseful to describe the population process conditioned on the eventual extinction of theprimary mutation in order to consider these two scenarios separately. The conditionalprocess can simply be described byPr( i → j | extinction) = Pr( i → j ) π j π i π i represents the the probability that the lineage goes extinct given that there arecurrently i mutants in the population (Ewens, 1973).I use first step analysis (Taylor and Karlin, 1984) to in order to find the total proba-bility that no successful secondary mutant is spawned from a lineage of primary mutantsdestined to eventual extinction. Let v i be the probability that no successful secondarymutants are spawned from a lineage beginning with i mutants. v i can be implicitlydefined as the probability that no successful mutants arise in the current time step mul-tiplied by the probability that no successful mutants are produced in the future. For theprocess conditioned on eventual extinction of the primary mutant lineage we have v i = (1 − ωr iN ) (cid:18) P r ( i → i − π i − π i v i − + P r ( i → i ) π i π i v i + P r ( i → i + 1) π i +1 π i v i +1 (cid:19) , (6)where the composite parameter ω = µ U ( a ). Note that v = 1 and π N = 0. This then isa system of N − N − T = N µ (1 − v ) π . For finite populations the system of equations can be represented as a tridiagonalmatrix. This system can be solved numerically using a mathematical computing pack-age. Many results are known for the matrix inverse and eigenvectors of tridiagonalmatrices (Usmani, 1994; da Fonseca, 2007). These results can be used to numericallycalculate the eigenvectors for even large population size because they involve less than4 N multiplications and additions. Because they use recursive calculations there is noconstraint imposed by memory levels and computation time basically increases linearly.On a MacPro with a 3.2 GHz Xeon processor running Mathematica 7 the calculation forpopulation size of 10 takes about 50 seconds.Given our system of equations (6) we can write a matrix equation of the form Av = x , where v is the vector of v i and A = a b c a b c a b c . . . . . .. . . . . . b N − c N − a N − x = x . The elements of the matrix can be found from equation (6) and I provide formulas for a i , b i , c i and x below. Note that this system has N − θ i = a i θ i − − b i − c i − θ i − ,θ = 1 ,θ = a φ i = a i φ i +1 − b i c i φ i +2 ,φ N = 1 ,φ N − = a N − . Each element of A − can be expressed as an algebraic expression of θ and φ . Because v is the first element of A − x and because only the first element of x is non-zero we needonly calculate A − , . A − , = θ φ θ N − . (7)Given the system of equations (6) we have x = ( N − N − ωr ) N ( N − r )(1 − π ) a i = ( N − i ) + i rN ( N − i + ir ) N − irωN − b i = ( N − i ) ir (1 − π i +1 ) N ( N − i + ir )(1 − π i ) N − irωNc i = ( N − i − i + 1)(1 − π i ) N ( N − i + ir )(1 − π i +1 ) N − ( i + 1) rωN Finally, the total probability that no successful secondary mutations are produced whilea lineage descending from a single primary mutant is extant is v = φ θ N − (cid:18) − (cid:16) − rωN (cid:17) N − N ( N − r )(1 − π ) (cid:19) . (8)This method can be used to numerically solve for v and is reasonably quick even in largepopulations.
3. Wright-Fisher Populations
In the Wright-Fisher formulation the population at generation t + 1 is found bysampling gametes produced in generation t . So long as the number of gametes producedis reasonably large, the probability distribution for adults in generation t + 1 is binomial6uch that P r (0 →
0) = 1 ,P r ( i → j ) = (cid:18) Nj (cid:19) (cid:18) irir + ( N − i ) (cid:19) j (cid:18) − irir + ( N − i ) (cid:19) N − j ,P r ( N → N ) = 1 , represents the probability that the population goes from i mutants to j mutants inone generation. Again I use a first step analysis to calculate the probability that nosecondary mutations arise beginning from a single mutant. Exact calculations of thefixation probabilities for the finite Wright-Fisher model are not available, so this systemcannot be converted to the process conditioned on eventual extinction of the primarymutant. This means that the probability of tunneling will have to be back-calculatedfrom equation (5). If the primary mutant becomes fixed then the probability that asuccessful secondary mutation is spawned is 1.I define ˜ v i as the probability that no successful secondary mutations are spawnedstarting from the state where i primary mutants are present. For each possible state i , the probability that no successful secondary mutants are produced is the probabilitythat none of the i primary mutants immediately produce a successful secondary mutant((1 − ω ) i ) multiplied by the sum of the probabilities that the next generation contains j primary mutants multiplied by the probability that a lineage starting with j mutantsnever produces a successful secondary mutant. This is slightly different from the Moranmodel where only one secondary mutant can possibly arise at each time point. This gives˜ v = 1 , ˜ v i = (1 − ω ) i N (cid:88) j =0 P r ( i → j )˜ v j , ˜ v N = 0 . The equations for 1 ≤ i ≤ N − v j terms as follows0 = (1 − ω ) i N (cid:88) j =0 P r ( i → j )˜ v j − ˜ v i = (1 − ω ) i i − (cid:88) j =0 P r ( i → j )˜ v j + P r ( i → i )(˜ v i −
1) + N (cid:88) j = i +1 P r ( i → j )˜ v j . This is a system of linear equations in ˜ v j and can be written in matrix form as A˜v = x where A = . . . . . . − ω ) P r (1 →
0) (1 − ω ) P r (1 → − − ω ) P r (1 → . . . . . . (1 − ω ) P r (2 →
0) (1 − ω ) P r (2 →
1) (1 − ω ) P r (2 → − . . . . . . ... ... ... . . . ...0 0 . . . . . . v = ˜ v ˜ v ...˜ v N x = . The solution to ˜v = A − x can be found numerically using standard mathematical pack-ages. Compared to the solutions for the Moran model, the Wright-Fisher model requiresmany more computational operations for the same population size.Recall that ˜ v is the unconditioned probability that no successful secondary mutationsare produced from a lineage of primary mutants descending from a single primary mutant.We can calculate the probability that no successful secondary mutations are producedconditioned on the eventual extinction of a lineage of primary mutants descending from asingle initial mutant using the approximate fixation probability for a single copy mutant(Crow and Kimura, 1970) and equation (5) as v = ˜ v (cid:18) − − e − r − − e − N ( r − (cid:19) .
4. Large Population Size Approximations
In very large populations, the dynamics of newly introduced mutants can be modeledas a branching process. In this limit, segregating mutants do not interact and mean fitnessis not altered by their spread in the population. The probability that a newly introducedprimary mutant leads to the eventual maintenance of a secondary mutant can then becalculated based on this branching process, while fixation probabilities of beneficial allelesmay still be modeled using finite population size results. In other words, population sizeenters the calculations in two ways; the finite population size causes frequency dependentinteractions between segregating mutant alleles and causes the fixation probability ofbeneficial mutants to increase. The large population size approximation takes this firstpopulation size to be large while leaving this second population size finite.The main feature of the branching process approximation that allows this to becalculated is that the probability that no secondary mutations persist as a function ofthe number of initial mutants scales as v i = α i , where α represents the probability thatno secondary mutations persist when a single primary mutant is introduced. Once α hasbeen solved for the per generation probability of tunneling is described by T = N µ (1 − α )(1 − U ( r )) . (9) The calculation of the expected probability of a successful secondary mutation in thelarge population limit begins with the same first step analysis equations as in a finitepopulation. Starting with equation (6) setting i = 1 and taking the limit as N → ∞ wehave v = v (1 + r (1 + ω )) − r . (10)8nserting v i = α i M into equation (10) gives α M = r ( ω +1) − √ (1 − r ) +(2+ r ( ω +2)) rω r . (11)Note that this is slightly different from the formula Iwasa et al. (2004) give be-cause theirs is an expression in terms of the unconditional process. Recalling that E [ P | extinction] = E [ P ](1 − U ( r )) their expression can be rewritten in terms of theconditional process as α IMN = 2 − (cid:112) (1 − r ) + 2(1 + r ) rω r (1 − U ( r )) . (12)In general, these equations are extremely similar and will often be so close as to benumerically indistinguishable for any empirical applications. Notable exceptions arewhen the population is small enough that the Iwasa et al. (2004) breaks down and α IMN > r = 1) the two expressions become α M = 1 + 12 (cid:16) ω − (cid:112) ω (4 + ω ) (cid:17) (13) α IMN = (cid:0) − √ ω (cid:1) N − N . (14)Recalling that ω = µ − /a − /a N and taking the limit as N → ∞ we getlim N →∞ α M = 1 − (cid:114) (1 + a − a µ ) − a − a µ (15)lim N →∞ α IMN = 1 − (cid:114) a − a µ (16)(See (Schweinsberg, 2008) for an alternate derivation). These expressions are most dif-ferent when a is large and µ is large. However, even for an unrealistically large value µ = 10 − the expressions are different by less than 1%. Using the branching process approximation we have v j = α j WF . For 1 ≤ i ≤ N − v i = (1 − ω ) i N (cid:88) j =0 (cid:18) Nj (cid:19) (cid:18) irir + ( N − i ) (cid:19) j (cid:18) − irir + ( N − i ) (cid:19) N − j α j WF . (17)For any specific i , v i can be approximated by noting that the binomial probabilitiesapproach a Poisson distribution as N → ∞ (Ross, 1988; Iwasa et al., 2003). Rewritingequation (17) in the limit of large N using the Poisson approximation gives v i = (1 − ω ) i ∞ (cid:88) j =0 e − ir ( ir ) j j ! α j WF = (1 − ω ) i (cid:16) e − r (1 − α WF ) (cid:17) i . (18)9 similar result was obtained by Iwasa et al. (2003) for multi-step pathways involvingmultiple routes to a higher fitness mutant. For a primary mutant lineage starting with1 individual the probability of tunneling is the solution of α WF = e − r (1 − α WF ) (1 − ω ) . (19)This implicit definition of α WF can be easily solved numerically for specific values of r and ω (figure B.1).The implicit equation for the probability that no successful secondary mutants areproduced has a convenient interpretation. In large Wright-Fisher populations the dis-tribution of offspring produced by a single parent is Poisson. A newly arising primarymutant has probability ω of immediately mutating and spawning a lineage of secondarymutants that is destined to fix. If it does not (with probability 1 − ω ) then it producesa Poisson distributed number of offspring mutants with mean r , each of which has aneffectively independent probability 1 − α WF of spawning a secondary mutant lineage itselfdestined to fix. The probability that none of these primary mutant offspring spawn asuccessful secondary mutant lineage is simply the 0 term from the Poisson distributionwith parameter r (1 − α WF ). Thus, α WF is equal to the product of the probability that thelone primary mutant does not immediately produce a successful secondary mutant withthe probability that none of the primary mutant offspring produce a successful secondarymutant lineage. Define T M and T WF as the probability that a successful secondary mutation arisesfrom a lineage of primary mutants founded by a single primary mutant in the Moran andWright-Fisher models, respectively.For the Moran model where r = 1, T M = (cid:112) ω (1 + ω/ − ω (i.e. 1 − α M fromequation (13)). For small ω this can be approximated using a series expansion (seeappendix Appendix A). Likewise the Wright-Fisher model can be approximated usingequation (19) and again expanding around small ω (see appendix Appendix A). T M ≈ √ ω − ω O [ ω ] / , (20) T WF ≈ √ ω − ω O [ ω ] / . (21)Noting that the fixation probability approaches U ( a ) = ( a −
1) in large Moran populationsand U ( a ) = 2( a −
1) in large Wright-Fisher populations and substituting ω = 2 µ U ( a )gives T M ≈ (cid:112) µ ( a − − µ ( a − , (22) T WF ≈ (cid:112) µ ( a − − µ ( a − . (23)These calculations show that for the same level of µ and a the Wright-Fisher tunnelingprobability is a factor of 2 larger than for the Moran model.In both models the probability of tunneling depends on the distribution of the num-ber of primary mutants spawned by a lineage founded by one primary mutant and the10urvature of the function relating the probability a successful secondary mutant arises tothe total number of primary mutants in a given lineage. In the limit of large populationsize, where a branching process approximation is valid, this distribution can be exactlydetermined (Dwass, 1969). Examining the relationship between the expectation of theprobability no successful secondary mutant arises and the distribution of the primarymutant lineage size can help to understand this effect. Once population size gets to bereasonably large, the distribution of primary mutant lineage size remains almost constantexcept at the tail. In the Wright-Fisher case, when r = 1 we have a critical branchingprocess. Conditioned on eventual extinction of the primary mutant lineage, we can calcu-late the expected number of primary mutants spawned in a lineage arising from a singleprimary mutant (see Appendix B). Alternatively, we could calculate the distribution ofprimary mutant lineage sizes and then calculate the total probability that no successfulsecondary mutations are spawned. Figure B.3 shows how the curvature of the functionfor the probability of successful secondary mutations as a function of primary mutationlineage size would lead to a reduced total probability that a secondary mutation becomesfixed. Because of Jensen’s inequality, the expectation of the function of the random vari-able is larger than the function of the expectation. Because the probability cannot bebelow 0, no matter how large the primary mutant lineage, the contribution of very largeprimary mutation lineages is negligible. Thus, T is considerably smaller than we wouldpredict based only on the expected number of primary mutants spawned by a lineagedestined to eventually become extinct.For cases where r < ω is more straightforward (seeappendix Appendix A). Approximating and substituting in the values for ω gives T M ≈ µ ( a − r − r (24) T WF ≈ µ ( a −
1) 11 − r . (25)For r near one this approximation fails breaks down and either the exact equation or theneutral approximation should be used (figure B.1). Again, the probability of tunnelingis larger for Wright-Fisher populations, but now they are different by a factor of 2 r .The factor of 2 is again because of higher fixation probabilities in the Wright Fishermodel while the factor of r is due to the Moran model assumption that mutations onlyoccur during reproduction. If, alternatively, mutations were assumed to happen at aconstant rate per generation (where generations encompass N elementary steps of theMoran process), then the results would only differ because of their fixation probabilities(to the first order approximation).These results agree with the Moran population results of Iwasa et al . (2004). Iwasa et al . (2004) characterize the regime where r < − r descendants (see Appendix B). An important caveat is that equations (24) and (25) only11pply when r < ω is small. In these cases equations (11) and (19) should be usedinstead.
5. Multiple intermediate steps
This approach readily extends to the case where multiple intermediate mutationalsteps stand between the ancestral state and any mutants that have improved fitness.Such multi-step probabilities have been calculated for the Moran model by Iwasa etal . (2004) , Schweinsberg (2008), Lynch and Abegg (2010); Lynch (2010), and Weissmanet al. (2009). Schweinsberg (2008) considered the neutral case and found that the ap-proximation is valid when µ is small relative to the inverse of the population size squared.Lynch (2010) considered both the multistep tunneling probability for both neutral anddeleterious intermediates but did so by assuming that each intermediate deleterious mu-tation first rose to mutation selection balance. Additionally, Iwasa et al. (2004) consid-ered a immune response escape mutants in a multi-step Wright-Fisher model. Weissmanet al. (2009) examined the probability of multi-step tunneling for an arbitrary numberof intermediate mutations and rigorously defined the regions of parameter space wherethe stochastic tunneling approximation is valid. They note that the stochastic tunnelingapproximation, as used here, is only valid when population size is large. Specifically, therequired population size goes up if there are many intermediate mutations and they areonly weakly selected against (Weissman et al., 2009).The probability that a multistep tunnel is opened can be calculated recursively. Usingequation (11), but writing ω generically as µU , I define a function g ( µ, U ) = 1 − r ( µU + 1) − (cid:112) (1 − r ) + (2 + r ( µU + 2)) rµU r . (26)For a 2-step process, the probability of a tunnel is simply g ( µ , U ( a )). Further recursionsgive the total probability for longer tunnels, so that a three-step tunnel opens withprobability g ( µ , g ( µ , U ( a ))).In the case where the sequence of mutations all have the same mutation rate wethink of the probability a k -step tunnel is opened by suppressing the variable µ in g andrecursively applying g to U ( a ). The solution can be found graphically by cobwebbingthe graph of g (figure B.2) (Adler, 1998). Longer tunnels have lower probabilities, andthe decrease in probability depends on the shape of g . If g (cid:48) (0) < g has no fixed point), butif g (cid:48) (0) > g (cid:48) (0) < r < µ .For Wright-Fisher populations the probability of tunneling was described implicitlyas a solution to a transcendental equation. However, the same recursive approach canbe taken to find successive tunneling probabilities. Substituting ω = µU into equation(19) and noting that at a fixed point of the recursion α = 1 − U gives1 − U = e − r ( U ) (1 − µU ) , (27)where solving for U gives the fixed point of the recursion. Decreasing r decreases thefixed point. Approximating around U = 0 and solving gives U ∞ = 2( µ + r − r (2 µ + r ) , (28)12here U ∞ represents the fixed point of the recursion. If r < − µ then there is no positivefixed point. This fixed point represents the infinite recursion for the probabilities andrelies on the assumption that the time-scales of fixation and mutation can be treatedseparately. However, Weissman et al. (2009) found that these conditions narrow as thethe length of the pathway considered increases. As the length of the pathway increases,the probability of a series of sequential fixation events increases. However, when equation(28) has no fixed point we can still infer that tunneling across long valleys will not occur.
6. Simulations of finite population
The approach taken by Iwasa et al . (2004) involved first approximating the Moranprocess by a small time-step approximation and then using special functions and heuristicarguments to arrive at an approximation for the rate of tunneling. This method implic-itly assumes large population size and explicitly ignores some higher-order terms. Myapproximation explicitly assumes large population size to derive a result for the limit aspopulation size goes to infinity. Numerical solutions can be used to evaluate the abilityof these large-population size approximations to predict the rate of tunneling in finitepopulations. I simulated the discrete time Moran and Wright-Fisher models in order toassess the accuracy of the different approximations. My simulation draws waiting timesfor mutations and then tracks individuals in populations while multiple mutations aresegregating. Once the secondary mutation has reached a significant size its fixation isvirtually guaranteed, and the simulation is stopped.For the Moran model, the Iwasa et al . (2004) solution and my approximation aredisplaced from the numerical solution in opposite directions (figure B.4). The Iwasa etal . (2004) approximation overestimates the waiting time to a secondary mutation. Bothapproximations do extremely well once population size is larger than about 100 (for theparameter values in figure B.4). For smaller population sizes the variance in the wait-ing time is so large that, in practice, it would be hard to distinguish the alternativeapproximations. It is interesting to note that, in terms of displacement from the ac-tual waiting time, the Iwasa et al . (2004) approximation performs better, even thoughit intentionally ignores some terms. This is apparently because the large-population ap-proximation underestimates the waiting time (because mean fitness is not altered by thespread of mutants in the large-population limit), and by chance the terms that Iwasa etal . (2004) exclude happen to push the approximation further off.For the Wright-Fisher model, I compared my solution with the one presented byLynch and Abegg Lynch and Abegg (2010) (figure B.5). Lynch and Abegg applied theIwasa et al . (2004) approximation to Wright-Fisher populations but modified the calcu-lations to adjust for the possibility of multiple mutational hits in a single generation. Myapproximation for the rate of tunneling in Wright-Fisher populations includes a factorof √ et al . (2004) Moran approximationto Wright-Fisher populations. At small population sizes, the sequential fixation pathwaydominates and the two approaches yield similar predictions. At intermediate populationsize the predictions are most different, but still have qualitatively similar patterns. Forthe parameters used in figure B.5, the two predictions always differ by less than about0 .
3% and are within about 8 million generations of each other. For intermediate popula-tion sizes, the finite population matrix solution does capture the behavior of the system.13s population size grows larger the simulated results move towards my approximation,just as the matrix solution converges towards my approximation.
7. Conclusions
Multi-step evolution is becoming more widely recognized as an important componentof the evolutionary process (Weinreich and Chao, 2005; Hermisson and Pennings, 2005;Kopp and Hermisson, 2009; Durrett and Schmidt, 2008; Lynch, 2010; Lynch and Abegg,2010). Previous analyses were derived for the Moran model (Iwasa et al., 2004; Schweins-berg, 2008; Weissman et al., 2009) and for specific instances of the Wright-Fisher model(Iwasa et al., 2003) . The Moran model results have been used in models of Wright-Fisherpopulations (Lynch and Abegg, 2010) and the distinction between the two models has notreceived much attention. The results I present here are derived using elementary meth-ods from stochastic processes theory. I present exact expressions for the limiting case oflarge population size in both Moran and Wright-Fisher populations. Some methods forefficiently numerically evaluating the finite population size solutions are presented.The Wright-Fisher and Moran models differ in two different but related ways. First,the branching process calculation of the probability of tunneling depends on a sum overthe number of mutant offspring produced by a single mutant (compare equations (11) and(17)). For the same mean selection against primary mutants, the variance in the numberof offspring under the Wright-Fisher model is 1/2 as large as under the Moran model.The tunneling probability also depends on ω , the product of the secondary mutationrate and the probability of fixation of the secondary mutation. Second, the fixationprobability also depends on the distribution of the number of offspring and again differsby a factor of two between the models. For the same increase in fitness, the probabilityof fixation is twice as large under the Wright-Fisher model as under the Moran model.Because these are both introduced inside a square-root function, the total effect is thatthe probability of tunneling is twice as large under the Wright-Fisher model as comparedto the Moran model. The general prediction is that when the offspring distribution haslower variance then the probability of tunneling will increase.For Wright-Fisher populations, the probability of tunneling is the solution to anexponential equation which has a simple intuitive explanation. The total probabilitythat no successful secondary mutants are produced is the product of the probabilitiesthat the initial primary mutant does not immediately produce a successful secondarymutant, 1 − ω and the probability that none of its primary mutant progeny produce asuccessful secondary mutant.My analysis shows that tunneling in the Wright-Fisher model is more likely than inthe Moran model. Stochastic simulations show that the Wright-Fisher approximationdoes indeed capture the mean behavior of the evolutionary process once populationsize is relatively large. The improvement over applying the Moran approximation tothe Wright-Fisher scenario is quite minor, but there is no added difficulty in using thiscorrect approximation in the future. Acknowledgements
The work was improved by numerous discussions including those with Ricardo Azevedo,Michael Lynch, Alexey Yanchukov, Leah Johnson, Daniel Weissman, and the Theory14unch group. The comments of three reviewers greatly improved the manuscript. Thiswork was supported by NSF grant EF-0742582.
Appendix A. Approximations around small ω The probability of tunneling in a large Moran population is T M = (cid:112) ω (1 + ω/ − ω . (A.1)The second term ( − ω ) is already linear and does not need further approximation. Wewould like to approximate this function for small ω as a power series of terms ω i/ . Define f (ˆ ω ) = (1 + ˆ ω ) and write the radical as (cid:112) ωf (ˆ ω ), where ˆ ω will later be set equal to ω .We will construct a Taylor’s series approximation around ˆ ω = 0. Note that f (0) = 1 and f (cid:48) (0) = 1 /
4, and all higher derivatives of f are 0. The series is as follows (cid:112) ωf (ˆ ω ) ≈ √ ω + ˆ ω ωf (cid:48) (0)( ωf (0)) − / − ˆ ω
2! ( ωf (cid:48) (0)) ( ωf (0)) − / + . . . (A.2)Taking just the zero order terms gives T M ≈ √ ω − ω O [ˆ ω ] . (A.3)After setting ˆ ω = ω , the general expression for the i th term in the Taylor’s expansionis s i = ( − ( i +1) ω ((2 i +1) / i ! (1 / i (1 / i − j =1 (2 j − . (A.4)This series converges so that (cid:80) ∞ i =0 s i = √ ω (cid:112) ω/
4. So a good approximation is T M ≈ √ ω (cid:112) ω/ − ω . (A.5)For the Wright-Fisher model, when r = 1, α is represented implicitly by α = e − α (1 − ω ). Recall that T WF = 1 − α and express it in terms of Lambert’s W (the solution of W ( z ) e W ( z ) = z ) gives T WF = 1 + W ( ω − e ) . (A.6)Many methods of approximating W are known, we use the power series expansion pre-sented by Corless et al. (1996), W ( z ) ≈ − p − / p + · · · , (A.7)where p = (cid:112) ez + 1). Substituting this back into equation (A.6) gives T WF ≈ √ ω − ω + O [ ω ] / . (A.8)For situations where r < T M = 1 − r ( ω + 1) − (cid:112) (1 − r ) + (2 + r ( ω + 2)) rω r . (A.9)15way from r = 1 this can be approximated using a Taylor’s series expansion to give T M ≈ rω − r − ( rω ) (1 − r ) + O [ ω ] . (A.10)For the Wright-Fisher model with r < α as a function of ω giving α ( ω ) = e − r (1 − α ( ω )) (1 − ω ) , (A.11)and note that α (0) = 1. Differentiating the equation with respect to ω gives α (cid:48) ( ω ) = rα (cid:48) ( ω ) e − r (1 − α ( ω )) (1 − ω ) − e − r (1 − α ( ω )) . (A.12)Setting ω = 0 gives α (cid:48) (0) = rα (cid:48) ( ω ) − ⇒ α (cid:48) (0) = − − r . (A.13)Thus, T WF = 1 − α ( ω ) ≈ ω − r = 2 µ ( a −
1) 11 − r . (A.14) Appendix B. Mean number of primary mutants spawned by a lineage des-tined to become extinct
Appendix B.1. Moran populaitons
First step analysis can be used to calculate the mean number of mutants spawnedby a lineage descending from a single initial mutant. Again I condition on the eventualextinction of the mutant linage. Note that in the Moran model, time is scaled by thepopulation size and this must be kept in mind when using results based on the numberof mutants produced to estimate the probability of a secondary mutation. The first stepanalysis yields the system of equations D = P r (1 → π π (1 + D ) + P r (1 → D ) + P r (1 → π π (1 + D ) (B.1) D i = P r ( i → i − π i − π i ( i + D i − ) + P r ( i → i )( i + D i ) + P r ( i → i + 1) π i +1 π i ( i + D i +1 )(B.2) D N − = P r ( N − → N − π N − π N − ( N − D N − ) + P r ( N − → N − N − D N − ) , (B.3)where D i is defined as the expected cumulative weight of the number of descendantsfrom a mutant lineage starting with i mutants, conditioned on eventual extinction ofthe lineage (Weissman et al., 2009). This weight represents the mutational opportunityfor the lineage of primary mutants and is scaled by the population size. The boundary16onditions are that D = 0 while D N does not need to be defined since no transition tostate N is possible. In the neutral case where r = 1 this reduces to D = 1 N (1 + D ) + (cid:18) − N − N (cid:19) (1 + D ) + N − N (1 + D ) (B.4) D i = i ( N − i + 1) N ( i + D i − ) + (cid:18) − i ( N − i ) N (cid:19) ( i + D i ) + i ( N − i − N ( i + D i +1 ) D N − = 2( N − N ( N − D N − ) + (cid:18) − N − N (cid:19) ( N − D N − ) . The solution D i = iN / N time steps in the Moran model we can say that on average a single mutant producesan effective number of N/ N/ r < i mutants must produce i timesas many descendants as a lineage starting with 1 mutant. I define D ∗ i as the number ofdescendants produced scaled to the population size such that D i = D ∗ i N . Solving for D ∗ i +1 and taking the limit as N → ∞ gives D ∗ i +1 = ( r + 1) D ∗ i − D ∗ i − − r . The additional conditions that D ∗ i i = D ∗ i − i − = D ∗ i +1 i +1 implies that D ∗ i = i − r . (B.5)Thus, for the Moran model, the scaled number of mutants produced by a single mutantapproaches 1 / (1 − r ). For finite populations, the value is within 1% of this limit forpopulations larger than 100 when r < . Appendix B.2. Wright-Fisher populations
For Wright-Fisher populations I calculate the number of descendants when the pri-mary mutation is neutral. Consider a population with N haploid genotypes. The firststep equations are D i = N (cid:88) j =0 P r ( i → j ) π j π i ( i + D j ) ,P r ( i → j ) = (cid:18) Nj (cid:19) (cid:18) iN (cid:19) j (cid:18) N − iN (cid:19) N − j ,π i = N − iN . Because the mutant allele is neutral there is no effect on mean fitness as the numberof primary mutants changes in the population. This means that, conditioned on non-fixation of the primary mutation, the number of descendants left by i primary mutants17ust just be i times the number left by 1 primary mutant. Thus D i = D i . Insertingthis relationship back into the system of equations gives D i = iD = 1 N − i N (cid:88) j =0 (cid:18) Nj (cid:19) (cid:18) iN (cid:19) j (cid:18) N − iN (cid:19) N − j ( N − j )( i + jD ) , which is can be written as the sum of terms involving the first and second (non-central)moments of the binomial distribution with parameters N and i/N . Simplification gives D i = i (cid:18) D ( N − N (cid:19) . Solving for D gives D = N . Thus D i = N i. (B.6)So for a haploid Wright-Fisher population the mean number of neutral mutants spawnedby a lineage destined to extinction is equal to the number of haploid genomes present inthe population.For deleterious mutations it is not possible to write the conditional process for finitepopulations because no closed-from solution for the fixation probabilities is available.Instead, I calculate the average number of descendants in a large population using thePoisson approximation. The first step equations are D i = N (cid:88) j =0 P r ( i → j )( i + D j ) ,P r ( i → j ) = e − ir ( ir ) j j ! . In very large populations the branching approximation applies and D i = iD , yielding D i = D i = i + D ir. Solving for D shows that D i = i − r . (B.7)18 probability that a wild-type allele mutates to produce a primary mutation µ probability that a primary mutant allele mutates to produce a secondary mutation N Number of haploid genomes in the population U ( x ) probability an allele with relative fitness x becomes fixed when initially present as a single copy U ∞ probability that a tunnel of infinitely many steps will open. r fitness of primary mutants relative to the wild-type a fitness of secondary mutants relative to the wild-type S probability that a primary mutation destined to become fixed arises in a given generation S probability that a secondary mutation destined to become fixed arises in a populationcomposed entirely of primary mutants T probability, in a population of wild-type alleles, of a primary mutant destined to (beforethe primary mutant becomes fixed) give rise to a secondary mutant that then becomes fixed π i probability of eventual extinction of a lineage descending from i primary mutants ω composite parameter equal to µ U ( a ) v i probability that no successful secondary mutations are produced from a lineage descendingfrom i primary mutants, conditional on the non-fixation of the primary mutation v vector of the probability that no successful secondary mutations are produced˜ v i unconditional probability that no successful secondary mutations are produced from a lineagedescending from i primary mutants α probability that no successful secondary mutants are spawned from a lineage descending froma single primary mutant α IMN approximate α derived by Iwasa et al . (2004) T M for the Moran model, the probability that a single primary mutant will produce a lineagethat produces a successful secondary mutant. T WF for the Wright-Fisher model, the probability that a single primary mutant will produce a lineagethat produces a successful secondary mutant.19 igures .000 0.005 0.0101 r 5. Ω A r P r ob a b ilit yo f t unn e li ng B Figure B.1: The solution for the probability of tunneling in large Wright-Fisher populations. PanelA shows the probability of tunneling as a function of ω and 1 − r (note that the selection coefficient s = 1 − r ). When 1 − r ≈ − r increasesthe probability of tunneling approaches that given by equation (25). Panel B shows the decrease of theprobability of tunneling as 1 − r goes up in blue ( ω = 0 . √ ω relative to 1 − r from equation (25). For this value of ω , the two curves are within5% once 1 − r > . (cid:45) (cid:45) (cid:45) (cid:45) U i U i (cid:43) Figure B.2: The probability of tunneling in multistep pathways under the Moran model can be foundby recursively applying the formula for the probability of tunneling. The diagonal line is shown in blackalong with the recursion formula from equation (26) for three different values of r . In all cases, µ = 10 − .The green curve is for r = 1, where each intermediate mutation has the same fitness as the ancestralallele. The graph can be cobwebbed as shown by the black arrows. The initial value is the probabilitythat the beneficial mutant at the end of the pathway becomes fixed, starting from a single individual.Regardless of the initial condition, the probability of tunneling through a long pathway converges on thevalue where the green curve crosses the black line. The blue curve is for slightly deleterious intermediatemutants with r = 0 . r = 0 . umber of primary mutants P r ob a b ilit y no s ec ond a r y m u t a n t s s u ccee d minimum maximum true total probabilityprobability given mean lineage size Figure B.3: A schematic illustration of the relationship between the probability of tunneling and themean number of primary mutants. The probability that no successful secondary mutants are producedis a decreasing function of the primary mutant lineage size with positive second derivative. Assumethat the primary mutant lineage takes on two possible values with equal probability (the minimum andmaximum). The red lines map from the primary mutant lineage size and have a mean given by the reddashed line. The probability that no successful secondary mutant is produced from events having themean lineage size is shown where the red dashed line meets the vertical axis. The true total probabilitythat no successful secondary mutants are produced is found by averaging the probabilities in the twodifferent lineage sizes which can be visualized by finding the midpoint between the blue horizontal lines.Thus, the true probability that no secondary mutants are produced is higher than that found from themean population size. This means that the probability of tunneling is smaller than would be expectedbased on mean population size alone.
16 64 256 1024 4096020 00040 00060 00080 000100 000 Population Size W a iti ng T i m e Figure B.4: Predicted and observed waiting times for tunneling under the Moran model. The black pointsshow the waiting time in units of generations for a successful secondary mutation to arise along with the95% confidence intervals. Each parameter set was simulated 10 times. The yellow curve represents theexpected waiting time when the tunneling pathway is ignored. The orange curve represents the Iwasa et al . (2004) approximation (their equation 9), while the red curve represents my new approximation.The green curve was constructed by numerically solving the system for finite population size. For largepopulation size, the two approximations agree. The parameter values were chosen to have extremelybeneficial secondary mutants so that the tunneling effect is large, µ = 10 − , µ = 10 − , r = 1, and a = 100. The simulations were stopped when the advantageous secondary mutation became fixed orthere were more than 200 secondary mutants (the probability of loss after that point is approximately10 − ). Similar results are obtained for r < (cid:64) N (cid:68) L og W a iti ng T i m e B (cid:64) N (cid:68) L og W a iti ng T i m e Figure B.5: Predicted and observed waiting times for tunneling under the Wright-Fisher model. Theblack points show the waiting time for a successful secondary mutation to arise along with the 95%confidence intervals. Because the Wright-Fisher simulations take much more time than the Moranmodel simulations, each parameter set was simulated only 1000 times, except for N = 10 which wasonly simulated 100 times. The yellow curve represents the expected waiting time when the tunnelingpathway is ignored. The orange curve represents the formula used by Lynch and Abegg based on theIwasa et al . (2004) approximation, while the red curve is based on my Wright-Fisher approximation.The green curve was constructed by numerically solving the system for finite population size. Becausethis requires solving a non-sparse matrix equation that is the size of the population, I was only able tonumerically solve the finite population size model for population size less than 10 . Panel B shows thegraph for these population sizes in more detail. For small population sizes, the simulated data match thenumerical prediction quite well. For population sizes larger than about 10 the observed values agreewith my Wright-Fisher approximation. Parameter values are µ = 10 − , µ = 10 − , r = 1, a = 1 . ITERATURE CITED
Adler, F.R., 1998. Modeling the Dynamics of Life. Brooks Cole, Pacific Grove, CA.Corless, R., Gonnet, G., Hare, D., Jeffrey, D., Knuth, D., 1996. On the Lambert W function. Advancesin Computational Mathematics 5, 329–359.Crow, J.F., Kimura, M., 1970. An introduction to population genetics theory. Harper & Row, NewYork.Draghi, J.A., Parsons, T.L., Wagner, G.P., Plotkin, J.B., 2010. Mutational robustness can facilitateadaptation. Nature 463, 353–355.Durrett, R., Schmidt, D., 2008. Waiting for two mutations: With applications to regulatory sequenceevolution and the limits of darwinian evolution. Genetics 180, 1501–1509.Durrett, R., Schmidt, D., Schweinsberg, J., 2009. A waiting time problem arising from the study ofmulti-stage carcinogen esis. Annals of Applied Probability 19, 676–718.Dwass, M., 1969. Total progeny in a branching process and a related random walk. J Appl Probab 6,682–686.Ewens, W., 1973. Conditional diffusion processes in population genetics. Theoretical population biology4, 21–30.da Fonseca, C.M., 2007. On the eigenvalues of some tridiagonal matrices. Journal of Computational andApplied Mathematics 200, 283–286.Hermisson, J., Pennings, P., 2005. Soft sweeps: Molecular population genetics of adaptation fromstanding genetic variation. Genetics 169, 2335–2352.Iwasa, Y., Michor, F., Komarova, N.L., Nowak, M.A., 2005. Population genetics of tumor suppressorgenes. J Theor Biol 233, 15–23.Iwasa, Y., Michor, F., Nowak, M.A., 2003. Evolutionary dynamics of escape from biomedical interven-tion. Proceedings of the Royal Society of London. Series B: Biological Sciences 270, 2573.Iwasa, Y., Michor, F., Nowak, M.A., 2004. Stochastic tunnels in evolutionary dynamics. Genetics 166,1571–9.Komarova, N., Sengupta, A., Nowak, M., 2003. Mutation-selection networks of cancer initiation: tumorsuppressor genes and chromosomal instability. J Theor Biol 223, 433–450.Kopp, M., Hermisson, J., 2009. The genetic basis of phenotypic adaptation i: Fixation of beneficialmutations in the moving optimum model. Genetics 182, 233–249.Lynch, M., 2010. Scaling expectations for the time to establishment of complex adaptations. P NatlAcad Sci Usa 107, 16577–16582.Lynch, M., Abegg, A., 2010. The rate of establishment of complex adaptations. Mol Biol Evol 27,1404–1414.Moran, P., 1962. The Statistical Processes of Evolutionary Theory. Clarendon Press, Oxford.van Nimwegen, E., Crutchfield, J.P., Huynen, M., 1999. Neutral evolution of mutational robustness. PNatl Acad Sci Usa 96, 9716–20.O’Fallon, B.D., Adler, F.R., Proulx, S.R., 2007. Quasi-species evolution in subdivided populationsfavours maximally deleterious mutations. PNAS 274, 3159–3164.Proulx, S.R., Adler, F.R., 2010. The standard of neutrality: still flapping in the breeze? Journal ofEvolutionary Biology 23, 1339–1350.Proulx, S.R., Phillips, P.C., 2005. The opportunity for canalization and the evolution of genetic networks.Am Nat 165, 147–62.Ross, S., 1988. A first course in probability. Macmillan, New York.Schweinsberg, J., 2008. The waiting time for m mutations. Electron. J. Probab. 13, no. 52, 1442–1478.Taylor, H.M., Karlin, S., 1984. An introduction to stochastic modeling. Academic Press, Orlando.Tong, A., Evangelista, M., Parsons, A., Xu, H., Bader, G., Page, N., Robinson, M., Raghibizadeh, S.,Hogue, C., Bussey, H., Andrews, B., Tyers, M., Boone, C., 2001. Systematic genetic analysis withordered arrays of yeast deletion mutants. Science 294, 2364–2368.Tong, A.H.Y., Lesage, G., et al., 2004. Global mapping of the yeast genetic interaction network. Science303, 808–813.Usmani, R.A., 1994. Inversion of Jacobi’s tridiagonal matrix. Computers & Mathematics with Applica-tions 27, 59–66.de Visser, J.A.G.M., Hermisson, J., Wagner, G.P., Meyers, L.A., Bagheri, H.C., Blanchard, J.L., Chao,L., Cheverud, J.M., Elena, S.F., Fontana, W., Gibson, G., Hansen, T.F., Krakauer, D., Lewontin,R.C., Ofria, C., Rice, S.H., von Dassow, G., Wagner, A., Whitlock, M.C., 2003. Perspective: Evolutionand detection of genetic robustness. Evolution 57, 1959–1972. einreich, D., Chao, L., 2005. Rapid evolutionary escape by large populations from local fitness peaksis likely in nature. Evolution 59, 1175–1182.Weissman, D.B., Desai, M.M., Fisher, D.S., Feldman, M.W., 2009. The rate at which asexual populationscross fitness valleys. Theoretical population biology 75, 286–300.Wilke, C.O., Wang, J.L., Ofria, C., Lenski, R.E., Adami, C., 2001. Evolution of digital organisms athigh mutation rates leads to survival of the flattest. Nature 412, 331–333.Wright, S., 1932. The roles of mutation, inbreeding, crossbreeding and selection in evolution. The SixthInternational Congress of Genetics, Edited by D.F. Jones. Brooklyn Botanic Garden, Menasha, WI. .einreich, D., Chao, L., 2005. Rapid evolutionary escape by large populations from local fitness peaksis likely in nature. Evolution 59, 1175–1182.Weissman, D.B., Desai, M.M., Fisher, D.S., Feldman, M.W., 2009. The rate at which asexual populationscross fitness valleys. Theoretical population biology 75, 286–300.Wilke, C.O., Wang, J.L., Ofria, C., Lenski, R.E., Adami, C., 2001. Evolution of digital organisms athigh mutation rates leads to survival of the flattest. Nature 412, 331–333.Wright, S., 1932. The roles of mutation, inbreeding, crossbreeding and selection in evolution. The SixthInternational Congress of Genetics, Edited by D.F. Jones. Brooklyn Botanic Garden, Menasha, WI. .