The pace of evolution across fitness valleys

Evolutionary dynamics is based on natural selection, mutation and genetic drift nowak:2006bo . It can be illustrated as the dynamics of a population in an abstract, typically high-dimensional fitness landscape. Since individuals with higher fitness produce more offspring, the average density of individuals is highest close to the fitness maxima. Many such features as the stationary population density in the fitness landscape or the mutation rate under which a population can still be concentrated around a fitness maximum have been addressed eigen:1977aa ; eigen:1989aa ; wilke:2005aa ; nowak:1992aa . An important question is how a population evolves towards a fitness peak via several intermediate states. If the intermediate states have the same fitness as the initial state, then evolution to higher fitness states is neutral at first and thus poses no significant problems nimwegen:2000 . If the intermediate states have lower fitness than the initial state, then a fitness valley has to be overcome and it is more difficult to reach the fitness peak weinreich:2006aa ; poelwijk:2007aa . In this case, population stuck on a local peak cannot escape by natural selection alone, because there is no evolutionary trajectory with successively advantageous mutations. Instead, neutral genetic drift becomes important.

Here, we consider the dynamics of these systems from a different perspective. We address the average time a population needs to transfer from one peak to another one. For small mutation rates and finite populations, we calculate this average time analytically. When mutation rates are high, we can describe the system by a set of differential equations and obtain the relevant times from a numerical integration of the differential equations. In this framework, the relevant question is how fast a population evolves traulsen:2007ee .

In particular, we can address the question whether a population evolves faster from one peak to another via $d$ mutations if

• (i)

  mutations have to occur in a certain order, i.e. only a 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 are identical in both the cases and equal to that of the initial state. Thus the only difference remaining is the number of available paths. When the order of mutations is not fixed then multiple paths are available and the evolutionary dynamics will be faster when compared to a single path. We can then ask the question: Does a population evolve faster on a narrow ridge or a broad valley? This implies that we move away from the simplest case mentioned above and decrease the fitness in the intermediate states of the multiple paths compared to the fitness in the intermediate states of the single path. We show how the pace of evolution depends on the depth of the valley, the number of intermediate states and the size of the population.

In general, evolutionary dynamics depends crucially on the size of the population. In a small population a single mutation will typically reach fixation or extinction before another mutation can arise. The population moves as a whole step by step on the fitness landscape. For large populations, even for small mutation rates usually multiple types arise at the same time. This results in a non-zero population density in many states at the same time. For intermediate mutation rates, the population can either move stepwise across the fitness landscape or move several stpdf without getting concentrated in one of the intermediate states. This phenomenon has been termed stochastic tunneling iwasa:2004aa . If the mutation rates are too small, tunneling does not occur because it is unlikely that a second mutation arises before the first one has reached fixation or has gone extinct. If the mutation rates are high, tunneling occurs trivially, because the system can be approximated by differential equations for the densities in the different states. These different scenarios including the limiting cases of stepwise evolution (typical for small populations) and continuous evolution (typical for large populations) can also be observed when the population size is kept constant, but the mutation rates are increased. For computer simulations increasing the mutation rate is more convenient than simulating huge populations for moderate mutation rates.

One important example for an evolutionary process in which the timescales are of crucial importance is the somatic evolution of cancer frank:2007aa . Cancer progression has been investigated mathematically since the 1950s fisher:1959 ; nordling:1953 ; armitage:1954 . Of special interest are the tumor suppressor genes knudson:1971aa ; michor:2004aa . In a normal cell, there are two alleles of the tumor suppressor gene. The mutation in the first allele is neutral if the second wild-type allele can sufficiently perform the function. Inactivation of both the alleles confers a selective advantage to the cell and can lead to cancer progression. This is an example in which the order of mutations does not matter. Many cancers also require certain particular mutations that initiate cancer growth and pave the way for the accumulation of further mutations vogelstein:2004aa . Recently, it has been shown that after cancer initiation, a large number of different mutations may be involved in cancer progression sjoeblom:2006sj ; wood:2007ri ; jones:2008ty ; jones:2008ec . So far, it is unclear if the mutations have to occur in a specific order or if there is more variation in the order beerenwinkel:2007aa ; gerstung:2008aa .

For simplicity, we consider only very simple fitness landscapes here in which the fitness in all the intermediate states is identical. In natural systems, these fitness values will differ and also the mutation rate may not be constant. In addition, sometimes the order of mutations will matter and sometimes, it will not. Thus, sometimes a particular mutation will be a prerequisite to obtain a new function, but sometimes new mutations do not require any prerequisites. For example, this is the case in the evolution of resistance to $\beta$ lactam antibiotics studied by weinreich:2005aa and weinreich:2006aa . However, here we focus on a very simple model to highlight the general aspects of the dynamics by analytical and numerical considerations.

This paper is organized as follows. We begin with the description of the two ways to order the mutations, the single path and the hypercube. We then derive analytical approximations of the fixation times for small mutation rates and discuss the effect of the different parameters on the fixation times. Next, we address the dynamics for intermediate and high mutation rates. Finally, we explore biological examples which can be modeled using this approach.

To model evolutionary dynamics in a haploid population of size $N$, we use the Moran process nowak:2006bo ; moran:1962ef . In each time step, one individual is selected at random, but proportional to fitness. It produces one offspring, which replaces a randomly chosen individual. In one generation, each individual reproduces on average once. During reproduction, mutations occur with probability $\mu$. We are interested in the time it takes until $d$ mutations reach fixation in the population, starting from a homogeneous population in the initial state without any mutants. Moreover, we aim to explore the dynamical features of this process. We restrict ourselves to two different cases that allow the derivation of some analytical results.

The analytical approach is only valid as long as the mutation rate is small, $\mu\ll N^{-2}$. For higher mutation rates, the population does not have to consist of at most two different types at any time. Instead, $d$ mutations can be fixed in the population without sequentially fixing one after the other. This process has been termed stochastic tunneling and is of great importance in the context of cancer initiation iwasa:2004aa ; michor:2004aa ; nowak:2004bb ; beerenwinkel:2007aa . Tunneling across fitness valleys is more likely than tunneling across a flat fitness landscape (see Fig. 4). Even for $d=2$, the evolutionary dynamics is characterized by a doubly stochastic process, which makes analytical approaches tedious iwasa:2004aa . As discussed above, for $\mu N^{2}\ll 1$ the population usually contains at most two different types. In this case, the probability of stochastic tunneling will be very small. On the other hand, for $\mu N>1$, at least one mutant is produced per generation. Thus, the probability of stochastic tunneling approaches $1$. For $N^{-2}<\mu<N^{-1}$, the mutations are sometimes fixed sequentially and sometimes via stochastic tunneling. Fig. 5 shows how the tunneling probability increases from $0$ to $1$ in this interval.

For intermediate mutation rates, it is likely that the population contains more than two different types. The types with beneficial mutations will compete for fixation. This process is termed clonal interference crow:1970ck ; fisher:1930fi ; muller:1932aa ; gerrish:1998aa ; park:2007aa . Clonal interference has been considered to slow down adaptation, but recently it has been shown that it can have a positive influence on a rugged fitness landscape gerrish:1998aa ; wilke:2004aa ; jain:2007aa .

The states in a single path can be characterized by the number of mutations. In the hypercube, the states are characterized by the types of mutations that have already occurred. Thus, there are many different types that have undergone a specific number of mutations. However, all types that have already accumulated $k$ mutations can be pooled into the error class $k$. The number of different types in the error class $k$ is given by $\binom{d}{k}=\frac{d!}{k!(d-k)!}$.

In a single path, a population is said to tunnel across a state if it passes through a state without ever reaching fixation in that state. Analogously, in a hypercube a population said to tunnel across an error class if it passes through that error class without ever reaching fixation in it. Within the error classes, tunneling can occur across individual states, but also across several states at once. This means that the whole population passes only across that particular state and not across any other, without ever reaching fixation in that particular state. Tunneling across an error class can also occur in a second way: the population can use all of the available states in the error class, but the total number of individuals in the error class never reaches $N$. Thus, the probability of tunneling via the individual states is always lower than the probability of tunneling across the error classes. Fig. 6 shows the relation between the different probabilities of tunneling in the hypercube with respect to the rate of mutation $\mu$ for the special case $d=2$.

Due to higher effective rates of mutation, the probability of tunneling across a hypercube is expected to be greater than or equal to the probability of tunneling across a single path. However, numerical simulations reveal that for $d=2$ the probability of tunneling in a single path is higher than in the hypercube. This is a special case: for $d=2$ in a hypercube, the number of states into which the initial state can mutate into is $2$. The effective rate of mutations is thus twice as much as in the single path. The number of states which can be mutated into next is one, both in the single path and the hypercube. Thus the rate at which the individuals are pushed into the first state is higher in hypercube than in the single path while the rate of individuals being pushed out is the same. Thus there is a higher probability of reaching fixation in the first error class in a hypercube (see Fig. 6). We only observe this effect for $d=2$, for $d>2$, the probability of tunneling is higher in a hypercube than in a single path, as expected (see Fig.  4).

For $\mu N>1$, the stochastic features of the dynamics become less important. In this case, the system can be described by a set of $d+1$ deterministic differential equations for the fraction $x_{k}(t)$ of the population that has $k$ mutations jain:2007aa . Obviously, we have $\sum_{k=0}^{d}x_{k}(t)=1$. Transitions out of state $0$ occur with probability $T_{0\to}=(1-\frac{x_{0}}{\phi}u_{0\to 0})x_{0}$, where $\phi=x_{0}+(1-x_{0}-x_{d})\,s+x_{d}\,r$ is the average fitness of the population. This includes all the reproductive events except for the one where a type $0$ is produced. Transitions into state $0$ occur with probability $T_{\to 0}=\frac{x_{0}}{\phi}u_{0\to 0}(1-x_{0})$. Thus, the fraction of individuals in the initial state follows the differential equation

The probability that an offspring is of type $k$ is given by $\lambda_{k}=\sum_{j=0}^{k}\frac{x_{j}r_{j}}{\phi}u_{j\to k}$. The difference between the hypercube and the single path only occurs in the quantity $u_{j\to k}$, which is given above for both cases. The sum in $\lambda_{k}$ is over all individuals with $k$ or less mutations and $r_{j}$ is the fitness of individuals with $j$ mutations. This leads to the differential equation for the fraction of individuals with $k$ mutations,

where $k=0,\ldots,d$. Of course, the special case $k=0$ recovers Eq. (15). This set of $d+1$ differential equations describes how the system moves from state $k=0$ to state $k=d$. In general, only a numerical solution of this system of equations is feasible. While this allows us to infer details of the dynamics, our main interest is the time required for fixation of $d$ number of mutations. Thus, we solve the differential equation numerically using a standard Runge-Kutta algorithm press:2007aa . To find an equivalent to the fixation time in a stochastic simulation, we average between fixation ($x_{d}=1$) and the time when there are on average less than $1$ individuals outside the final state ($x_{d}=1-\frac{1}{N}$). Thus, the fixation time is the time when the solution of the differential equation crosses $x_{d}=1-\frac{1}{2N}$.

Fig. 7 shows an overview of the fixation times, covering the full range of mutation rates. For small mutation rates, we have sequential fixation of mutations and the time can be well approximated by Eqs. (8) and (10). For high mutation rates, the numerical solution of Eq. (16) leads to a good approximation for the fixation times.

We have determined the average time during which a population moves from a certain initial state to a final state of higher fitness. The initial and the final states are separated by a fixed number of mutations $d$. The mutations jointly confer a fitness advantage to the final mutant which can be represented by a peak in the fitness landscape. If the intermediate mutations need to occur in a specific order for the evolution of the final mutant then it corresponds to the single path. Otherwise, evolution occurs on a hypercube and there are $d!$ ways of reaching the final state.

We have explored the simplest system in which the fitness in all intermediate states is the same. As expected, the fixation times on a hypercube are shorter than on a single path, due to the presence of multiple paths available in a hypercube. This observation leads to the question: for which parameters does the hypercube show shorter fixation times than the single path, even with an added disadvantage? The fitness in the intermediate states was then set to lower values than the ones in the single path. Up to a certain threshold value of the fitness of the intermediate states, the hypercube shows shorter fixation times than in the single path. The value of the threshold depends on the population size, total number of required mutations and the fitness in the final state.

The fixation times for large populations largely depend on the fitness function and are qualitatively independent of the order of mutations. Let us first focus on a flat landscape: when the intermediate states have a fitness equal to the fitness of the initial wild-type, then for small mutation rates large populations have shorter fixation times than small populations. This is because the neutral rate of evolution does not depend on the population size. But the waiting time for fixation of the last mutation becomes shorter with larger population size. For intermediate mutation rates, tunneling starts earlier in larger populations. This leads to a marked decrease in the fixation time with larger population size. For high mutation rates, the time to fixation is no longer dominated by the time for the first mutant to reach the final state, but by the time until all individuals are in that state. Due to this, for high mutation rates the time required for fixation can be shorter in smaller population as compared to larger populations. Next, we focus on fitness valley: If the fitness landscape consists of a valley with reduced fitness of the intermediate states, small populations have an advantage for small mutation rates, as they can easily leave the initial state and enter the valley. But for high mutation rates, large populations reach fixation faster, because they can explore states within and beyond the fitness valley more easily.

Our numerical simulations reveal that tunneling can be 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 even in relatively large populations. Concrete values for fixation times are collected in Table 1. They reveal that even in long-term studies of experimental evolution, it is difficult to observe the consecutive fixation of neutral mutants cooper:2003aa . Consecutive fixation of advantageous mutants, however, is significantly faster. For example, while Table 1 reveals a fixation time of $\sim 10^{11}$ generations on a single path for $d=10$, $s=1$ and $N=10^{6}$, an optimal choice of the intermediate fitness values traulsen:2007ee would lead to a fixation time of $\sim 10^{7}$ generations.

While we have focused on the simplest possible system which allows analytical approximations, experimental studies reveal of course a much higher complexity. weinreich:2006aa studied experimentally the point mutations in the $\beta$-lactamase gene of bacteria. $\beta$ lactam antibiotics are commonly used, but the bacteria can develop resistance to the drugs. Five point mutations in a particular allele of the $\beta$-lactamase gene increases the resistance of the bacteria to cefotaxime by a factor of $\sim 100,000$. Theoretically the mutations leading from the wild-type allele to the resistant allele can occur in $5!=120$ ways. These can be represented by a hypercube of $d=5$. But in only $18$ of the $120$ trajectories, the intermediate mutations are either neutral with respect to the initial state or beneficial. Weinreich and colleagues have shown that these have the highest probability of realization. For all beneficial intermediates the fastest way to reach the final state would be when the relative fitness increase between any two consecutive mutations is constant traulsen:2007ee , but usually in nature several different mutations are available and the population first evolves to states that provide the highest selective advantage.

In another experimental study the sequence space of the 5s rRNA of a marine bacterium, Vibrio proteolyticus was explored lee:1997 . The sequences from Vibrio proteolyticus and Vibrio alginolyticus differ in only four positions. All the possible intermediates were constructed by the authors and the fitness of each was calculated chao:1985aa . Two of the valid intermediates have a fitness lower than the initial wild-type. We have shown how such fitness valleys can be crossed by exploring the phenomenon of tunneling or multiple mutations (for high mutation rates). Thus, the population does not need to move in a Wrightian fashion (the whole population moving as a whole across the valley).

The theory discussed herein deals with basic evolutionary concepts which are important to the kind of biological examples described above. More complex properties of the experimental studies like more general cases of epistasis and compensatory mutations can easily be incorporated, 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 $2^{d}!$ distinct epistatic patterns. Thus, one should rather focus on concrete systems instead. For example, one could simulate the dynamics in a system with experimentally derived fitness values and mutation rates. Not all the paths of a hypercube might be accessible for selection, but still some of them might prove to be significant depending upon the particular values of the parameters, such as fitness values and population size. Our goal here was to characterize the simplest features of the dynamics of a population crossing a fitness valley. This approach can be helpful when more realistic scenarios are addressed.

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