The competition between simple and complex evolutionary trajectories in asexual populations

On rugged fitness landscapes where sign epistasis is common, adaptation can often involve either individually beneficial “uphill” mutations or more complex mutational trajectories involving fitness valleys or plateaus. The dynamics of the evolutionary process determine the probability that evolution will take any specific path among a variety of competing possible trajectories. Understanding this evolutionary choice is essential if we are to understand the outcomes and predictability of adaptation on rugged landscapes. We present a simple model to analyze the probability that evolution will eschew immediately uphill paths in favor of crossing fitness valleys or plateaus that lead to higher fitness but less accessible genotypes. We calculate how this probability depends on the population size, mutation rates, and relevant selection pressures, and compare our analytical results to Wright-Fisher simulations. We find that the probability of valley crossing depends nonmonotonically on population size: intermediate size populations are most likely to follow a “greedy” strategy of acquiring immediately beneficial mutations even if they lead to evolutionary dead ends, while larger and smaller populations are more likely to cross fitness valleys to reach distant advantageous genotypes. We explicitly identify the boundaries between these different regimes in terms of the relevant evolutionary parameters. Above a certain threshold population size, we show that the probability that the population finds the more distant peak depends only on a single simple combination of the relevant parameters.


Background
In an adapting population, evolution often has the potential to follow many distinct mutational trajectories. In order to predict how the population will adapt, we must understand how evolution chooses among these possibilities. Many experimental and theoretical studies have analyzed this question, focusing primarily on the simple case where epistasis is absent, so that each mutation has some fixed fitness effect [1][2][3][4][5][6]. This work can explain the probability that a given mutation will fix as a population adapts, as a function of its fitness effect, the population size, mutation rate, distribution of fitness effects of However, the fitness effect of a mutation often depends on the genetic background in which it occurs. A particularly interesting form of this phenomenon, sign epistasis, occurs when several mutations are individually neutral or deleterious but their combination is beneficial [7]. Sign epistasis has been observed repeatedly in experiments [8][9][10][11][12][13], and plays a central role in the evolution of complex phenotypes that involve multiple interacting components. When sign epistasis is present, adaptation can involve passing through genotypes of lower fitness -i.e. a population may have to cross a fitness valley or plateau. Thus the fate of a mutation depends not only on its fitness, but also on its adaptive potential [14].
Several recent theoretical studies have analyzed the evolutionary dynamics of fitness valley crossing [15][16][17][18][19][20]. This work has focused on calculating the rate at which adapting populations cross a valley or plateau, in the absence of any other possible mutational trajectories. However, individually beneficial mutations may often compete with more complex evolutionary trajectories. We must then ask how likely evolution is to eschew the immediately uphill paths, and instead cross valleys or plateaus to reach better but less accessible genotypes. In other words, when the fitness landscape is rugged, we wish to understand whether evolution will take the more "farsighted" path to reach distant advantageous genotypes, rather than a "greedy" trajectory that fixes immediately beneficial mutations regardless of whether these may lead to evolutionary dead ends.
In this article, we analyze this evolutionary choice between immediately beneficial mutations and more complex mutational trajectories that ultimately lead to higher fitness. We calculate the probability that an adapting population will follow each type of competing trajectory, as a function of the population size, mutation rates, and selection pressures. We focus on asexual populations, where the only way for a population to acquire a complex adaptation is for a single lineage to acquire each mutation in turn. Our analysis is similar in spirit to earlier work which also considered the tradeoff between short-term and longterm fitness advantages [21][22][23][24]. However, these earlier studies dealt with competition between different strictly uphill or neutral paths, and considered the case where the less beneficial initial mutation led to better longterm evolutionary opportunities. In contrast, our analysis describes the competition between uphill mutations and more complex trajectories. While these two cases can be qualitatively similar in very small populations, they lead to very different dynamics in larger populations where the sign of the effect of the intermediate mutation can play a crucial role.
Our results show that population size has a crucial impact on how "farsighted" evolution can be. This dependence is not monotonic: evolution at intermediate population sizes is most "greedy", while both larger and smaller populations are more likely to eschew uphill paths in favor of complex trajectories. In large populations, our results show that a single parameter reliably predicts the extent of this evolutionary "foresight" across a wide range of parameters. Finally, we describe how our analysis can be generalized to predict how evolution will choose among even more complex trajectories, such as broad fitness valleys with multiple intermediate genotypes, and we discuss evolution in genotype spaces with many possible evolutionary paths.

Methods
We are interested in how a population makes an evolutionary choice when confronted with multiple possible mutational trajectories. Specifically, we focus on the extent to which adaptation proceeds by crossing fitness valleys rather than acquiring immediately beneficial (uphill) mutations. Of course, the relative frequency of valley crossing will depend on the number of available fitness valleys, their depth, and the fitness advantage of the multiple-mutants, as well as the distribution of fitness effects (DFE) of the uphill mutants. Our goal is to understand how the prevalence of valley crossing depends on these factors.

Model
Throughout most of this article, we consider the simplest context in which we can address this question: the choice between a single uphill path and a single fitness valley. Specifically, we consider a haploid asexual population of constant size N which can either acquire an uphill mutation (u) that confers an immediate fitness advantage s u , or alternatively acquire a deleterious fitness valley intermediate (i) with fitness deficit δ i on which background a double-mutant (v) with fitness s v > s u can arise. This scenario is illustrated in Figure 1. We also consider the case of a fitness plateau, where δ i = 0.
Because we are interested in the evolutionary choice between competing mutational trajectories, we assume that these two trajectories are mutually exclusive (i.e. the mixed genotypes ui and uv are strongly deleterious), so that only one genotype (either u or v) can eventually fix in the population. As a measure of evolutionary foresight, we analyze the probability that the double-mutant v fixes as a function of the relevant mutation rates, selection coefficients, and population size. In some situations, we could imagine that after either genotype u or v fixes, another set of competing potential trajectories become available. In this case, our analysis predicts the long-term relative ratio of fixed uphill versus valley-crossing mutations. In the Discussion, we consider how this model can be extended to the situation where there are many different competing uphill paths and valleys, and to broader fitness valleys involving multiple intermediate genotypes.

Simulations
In addition, we compare our analytical predictions for valley crossing probability to Wright-Fisher simulations. Each simulated population was evolved until either the uphill genotype or valley-crossing genotype fixed. Valley crossing probabilities were then inferred from the number of trials in which the valley-crossing genotype fixed, out of 1000 trials per parameter set.

Results
In the absence of the uphill genotype, fitness valley crossing can be modeled as a homogeneous Poisson process with rates as calculated by [17]. In small populations, the primary role of the uphill genotype is to introduce an effective time limit on this process: once an uphill mutation destined to survive drift first occurs, it very quickly fixes, leading to the extinction of the wild-type. The probability of valley-crossing can thus be calculated as the probability that the intermediate i fixes before the uphill genotype u. An example of this is shown in Figure 1b.
In larger populations, the dynamics are more complex, as illustrated in Figure 1c. Rather than leading to a single cutoff time for valley-crossing to occur, the single-mutant occurs and gradually increases in frequency. This leads to a decline in the size of the wild-type background on which intermediate and valley-crossing mutants can arise, and a corresponding increase in the mean fitness of the population ( Figure 1c). These effects gradually reduce the rate at which intermediates are produced, and make these intermediates effectively more deleterious relative to the mean fitness. These factors reduce the rate of the valley-crossing process. Thus valley-crossing becomes an inhomogenous Poisson process, with a rate that depends on the random appearance time T u of the uphill mutant.
In general, these effects of interference and tunneling are complex. However, the analysis becomes simpler in two specific regimes. When the expected drift time of the intermediate genotype is short, we can neglect the changing background fitness due to the uphill mutant during this drift time ( Figure 1d). Alternatively, for very large populations (Nμ > 1), the Poisson process approximation breaks down and both uphill and intermediate mutations can be treated deterministically (Figure 1e), and only the valley-crossing genotype must be treated stochastically.
These various regimes are illustrated in Figure 2. We now analyze each in turn, assuming weak selection (s j < 1) for all genotypes throughout. Taken together, this provides a complete picture of the probability that evolution will eschew the immediately uphill path in favor of the more complex adaptation.

Small populations
When the population size is small enough that the probability of stochastic tunneling is very low, the population is generally clonal or nearly clonal, and moves in Markovian jumps between neighboring genotypes. The transition between genotypes i and j occurs at rate where π ij is the probability that a single j mutant will give rise to a lineage that fixes, given by the standard formula, We refer to this as the sequential fixation regime. Because we are considering neutral and weakly deleterious intermediates, we account for the possibility of back-mutation to the wild type if the intermediate fixes. Therefore, the process can be modeled as an absorbing states Markov chain, where the wild type and intermediates act as transient states, and the uphill and double-mutant genotypes act as absorbing states. From elementary Markov chain theory, we find As the population size increases, π wu → 2s u and π wi → 0, so π wu π wi → ∞, and P cross → 0. Thus we find that within the sequential fixation regime, larger population sizes are less likely to cross fitness valleys.

Stochastic tunneling
For large populations, the probability that deleterious intermediates will fix declines drastically, and successful double mutants will instead arise on the unfixed singlemutant background in a process known as stochastic tunneling [16]. This transition occurs when Here N wt is the wild-type population size, and p i (t d ) is the cumulative probability that a single-mutant lineage will give rise to a successful double-mutant lineage by time t d after it appears. This probability is given by [17]: where and δ i,eff and s v,eff are the fitnesses of the intermediate and valley-crossing genotypes relative to the (time-dependent) mean population fitness. These effective fitnesses and N wt depend on the background at time t d in a way we must now consider. The background in turn is determined by the frequency f u (T v + t d ) of the uphill genotype at time t d after the appearance of the first single-mutant intermediate lineage destined for success at time T v (Figure 1c). Thus the uphill genotype frequency sets the fitness background on which valley-crossing probabilities are determined.
To calculate the relevant effective parameters, we note that the appearance of uphill lineages destined for success can be modeled as a Poisson process. Moreover, because the valley-crossing genotypes make up a tiny fraction of the population (unless the double-mutant has already established), we can treat the genetic background on which these uphill lineages appear as essentially fixed. Therefore, the first uphill lineage destined to survive genetic drift will appear at time T u , distributed exponentially with rate λ u = Nμ u π wu ≈ Nμ u (2s u ).
Once a successful uphill lineage appears, we assume it establishes in time τ est = γ e /(2s u ), where γ e ≈ .577 is the Euler-Mascheroni constant [25], and then sweeps deterministically according to wheret ≡ t − T u − τ est is the time after establishment of the uphill mutant.
Conditioning on the appearance time T u , we can thus work out our effective parameters (12) These effective parameters encompass the two main effects of the sweeping uphill mutation on the valley crossing probability: the first represents the declining wild-type background on which new mutations can arise, and the remaining two represent the decreasing relative fitness of the valley-crossing lineage.
We are interested in the probability that a doublemutant lineage destined for success appears before the uphill genotype fixes. Integrating over all possible appearance times T u , this is given by: This integral is a complete solution for the probability of valley-crossing in the stochastic tunneling regime, provided that the population size is small enough that the Poisson process approximation above holds. Although it does not have a simple closed-form solution, we can easily evaluate the integral numerically. Alternatively, there is a simple and relevant parameter regime in which background fitnesses change slowly. We now consider this case, and show that it allows us to evaluate our expression for the valley-crossing probability explicitly. In a later section below, we turn to the alternative case where the Poisson process approximation breaks down, and we can instead treat all single-mutants deterministically; the valley-crossing probability also simplifies considerably in this very large population regime.

Slowly changing background fitness
One of the main complications of Equation (13) is the integral over possible drift times t d , which reflects the increasing effective deleteriousness of the intermediate genotype as the uphill genotype sweeps to fixation and increases the mean fitness of the background population. However, when s u is small or δ i is large, this background fitness changes slowly compared to the intermediate drift time. In this case, we can treat the background during intermediate drift as effectively constant (Figure 1d). This eliminates the need for an integral over t d , since the timedependent probability of success of a single-mutant at time t is fully determined by We can further simplify the analysis if we treat the probability of crossing the valley as a function of two probabilities: the probability P v,1 that the first successful valley-crossing lineage appears before the first successful uphill lineage establishes, and the probability P v,2 that a successful valley-crossing lineage appears after the uphill mutant establishes: The calculation of P v,1 takes place on a purely wild-type background, so we can use where p v , the probability that the intermediate lineage survives drift long enough to give rise to an ultimately successful double mutant lineage, is given by [17]: Meanwhile, λ u is unchanged from the original analysis. P v,1 is determined by a race between these two exponential random variables. Using basic properties of the exponential, the probability that the double-mutant appears before the uphill genotype establishes is therefore where τ represents the difference between the mean drift time τ drift of the valley-crossing lineage and the mean establishment time τ est of the uphill lineage, Here τ est is as given above, and from [17], we can approximate the drift time as If the successful uphill lineage establishes, the first successful valley-crossing lineage still has a chance to appear and outcompete it, albeit on a declining wild-type background. Thus for P v,2 we get a similar integral as in the original analysis. However, since we are approximating p i as constant, the rate of successful lineage generation simplifies tô (20) Integrating and assuming mutation rates are small compared to selection pressures, we find: Combining these results, we find We expect this result to be valid provided that δ i is effectively constant over the expected drift time. This will hold when

Semi-deterministic tunneling
We now consider the case where Nμ i > 1 and Nμ u > 1, and hence the Poisson process approximation used to derive Equation (13) breaks down. Fortunately, in this regime the number of single-mutant intermediates and uphill mutants in the population are well approximated by their deterministic expectation (Figure 1e). Thus the only random variable is the appearance time of the first successful double-mutant lineage, which occurs with rate Because intermediates never make up a large portion of the population, f u is unaffected by f i , and hence is still given by Equation (9). We can then approximate the frequency f i of the single-mutant intermediates using mutation-selection balance with a declining wild-type population: where f i * gives the independent deterministic dynamics (mutation-selection balance) of single mutants on the wild-type background, and (1 − f u ) is the size of the wildtype background. It is useful to transform this into an integral in the frequency domain: We note from equation 9 that We further approximate Combining these expressions and assuming μ u s u , we find where This gives where we have defined the useful quantity Thus we see that in very large populations, the probability of valley-crossing depends in a simple way on the single composite parameter . The form of this composite parameter depends crucially on whether the fitness cost of the intermediate genotype is large or small, as defined in Equation (31).

Discussion
Our results have shown that the size of a population strongly influences how "farsighted" it can be. In small populations, genetic drift is strong relative to selection, so the evolutionary dynamics proceeds by sequential fixation. Since fixation of a deleterious intermediate becomes less likely in larger populations, this means that increasing the population size initially decreases the relative influence of fitness-valley crossing. However, as the population size increases further, beneficial mutations take longer to fix, maintaining diversity in the population and allowing double mutants to stochastically tunnel on the declining wild-type background [15][16][17]. Together, these effects lead to a non-monotonic relationship between population size and the probability that evolution will favor the complex adaptation over the directly uphill path, as illustrated in Figure 3a. This nonmonotonic dependence on population size is similar in spirit to the results of earlier work analyzing evolution on epistatic landscapes in the absence of fitness valleys [21][22][23].
It is interesting to note that P cross does not immediately begin to rise with the onset of tunneling. Instead, the dependence is more complex, as a consequence of the tradeoff between increasing mutation rates and fixation times. Nevertheless, for populations in which the transition to valley-crossing behavior occurs in the semideterministic regime, we can derive a simple expression for the threshold size at which the population will tend to cross valleys with probability P cross . A straightforward inversion of (32) gives valid for large population sizes in the semi-deterministic regime. Thus in this regime the threshold size above which a population exhibits a given degree of foresight (i.e. has a particular P cross ) depends only on . To illustrate this, in Figure 3b we show P cross as a function of N for a variety of simulations across different values of μ u , μ i , μ v , s u , δ i and s v . It is clear that even across a wide parameter range, N is a reliable predictor of valley crossing probability in the semi-deterministic limit.

Multiple intermediates and evolutionary predictability
Recently, Szendro et al. [26] simulated evolution across a wide range of population sizes on an experimentallyderived epistatic fitness landscape, finding that "evolutionary entropy" (i.e. unpredictability over all possible outcomes, given by S = − p j log p j ) varied nonmonotonically with population size. Specifically, these authors found that entropy initially decreased with N above a characteristic population size N ∝ 1/μ, before increasing again above a second characteristic size N ∝ 1/μ 2 ; they argued that these points were related to the supply rate of single and double mutants respectively. Our analysis is consistent with these results. For example, the increase in entropy at N ∝ 1/μ 2 found by [26] corresponds to our result that valley-crossing begins to significantly influence evolution when N ∼ 1/ : this is approximately proportional to 1/μ 2 , albeit with an additional log dependence on μ from the γ factor that would be harder to observe experimentally.
We find related behavior if we extend our analysis to valleys with more intermediates. As a simple example, we consider a fitness landscape where we add deleterious intermediates u 0 and v 0 (each occurring at rate μ 0 and leading to fitness cost δ i ) to the uphill and valleycrossing branches respectively, so that we now have competition between a single-intermediate valley and a two-intermediate valley. In a large enough population, mutation-selection balance will ensure that If we assume that these sub-populations are large enough that double-mutants behave deterministically, then we find the crossing probability obeys Thus our analysis of valleys with multiple intermediates suggests that the Nμ 2 entropy peak is not unique: as the population grows larger, there should be entropy peaks corresponding to foresight across valleys with increasing numbers of intermediates. The emergence of such second peaks has been observed in simulations [23], and our model offers a quantitative outline of where such peaks should occur given the relevant evolutionary parameters. In the example above, for instance, there would be an entropy peak approximately proportional to Nμ 3 , and in general, we could expect entropy peaks at points approximately proportional to Nμ n for (n − 1)-intermediate valleys. In practice, however, the semideterministic approximation will break down for any sizable number of intermediates, unless the population size is unrealistically large.

Many paths
Throughout this paper, we have assumed the presence of a single uphill mutation and fitness valley. We now consider how our analysis can be extended to predict how evolution chooses among many such possible mutational trajectories.
In small populations that are in the sequential fixation regime, we simply add additional transient transition matrix elements representing different mutations, with the uphill mutations transitioning to the uphill absorbing state, and similarly for the valley-crossing mutations. When stochastic tunneling is important, we must instead add the rates of single valley-crossing mutants to get a total rate of and similarly find the total rate of uphill mutants that are destined to survive drift, where in the last equality we have replaced a discrete collection of uphill mutations with a continuous distribution of uphill fitness effects ρ(s), and a discrete collection of mutation rates μ u with a total beneficial mutation rate U b . This is valid as long as u 1. Once an uphill mutation destined to survive drift occurs, the probability that it has fitness s is given by the ratio between its partial rate and the total rate; formally, the probability density is given by: Using these expressions, we can integrate our results from the analysis over all possible trajectories. However, we note that if there are a large number of weakly beneficial mutations, it is possible the first successful lineage to appear will be outcompeted by a stronger uphill mutation that arises later but fixes first. Our analysis applies provided that we consider only uphill mutations that reach a significant portion k of the population before a new, more fit uphill mutant is expected to be produced: i.e. u τ k = ∞ s cutoff Nμ s π(s)ρ(s)ds (τ k ) < 1, where τ k is the expected time for a single-mutant destined for success to make up frequency k of the population. This is consistent with our intuition that as the population size grows larger, we increasingly expect the mutations of largest effect to dominate the dynamics.