Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2711507

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Model
- 3. Small mutation rates
- 4. Intermediate mutation rates
- 5. High mutation rates
- 6. Discussion
- References

Authors

Related links

J Theor Biol. Author manuscript; available in PMC 2010 August 7.

Published in final edited form as:

Published online 2009 April 24. doi: 10.1016/j.jtbi.2009.04.011

PMCID: PMC2711507

NIHMSID: NIHMS112932

The publisher's final edited version of this article is available at J Theor Biol

See other articles in PMC that cite the published article.

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

Evolutionary dynamics is based on natural selection, mutation and genetic drift (Nowak, 2006). 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 and Schuster, 1977; Eigen et al., 1989; Wilke, 2005; Nowak, 1992). 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 (van Nimwegen and Crutchfield, 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 et al., 2006; Poelwijk et al., 2007). 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 et al., 2007).

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

- mutations have to occur in a certain order, i.e. only a single evolutionary trajectory is available or
- 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 steps without getting concentrated in one of the intermediate states. This phenomenon has been termed stochastic tunneling (Iwasa et al., 2004). 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 continuos 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, 2007). Cancer progression has been investigated mathematically since the 1950s (Fisher, 1959; Nordling, 1953; Armitage and Doll, 1954). Of special interest are the tumor suppressor genes (Knudson, 1971; Michor et al., 2004). 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 and Kinzler, 2004). Recently, it has been shown that after cancer initiation, a large number of different mutations may be involved in cancer progression (Sjöblom et al., 2006; Wood et al., 2007; Jones et al., 2008a,b). 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 et al., 2007; Gerstung and Beerenwinkel, 2008).

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 *β* lactam antibiotics studied by Weinreich (2005) and Weinreich et al. (2006). 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, 2006; Moran, 1962). 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 *μ*. 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.

If the mutations can occur only in a particular order, we have a single evolutionary path, see Fig. 1 for an illustration. Individuals in the initial state have fitness *r*_{0} = 1 and individuals in the final state have fitness *r _{d}* > 1. It is instructive to characterize an individual by a string of

If the order of mutations does not matter, evolutionary dynamics takes place on a hypercube in *d* dimensions cf. Fig. 1. Thus, there are 2^{d} different types of individuals. In the initial state, we have *d* possible mutations. In the next step, *d* − 1 mutations are available. Consequently, we have *d*! possible paths to fixation. Again, we assume *r*_{0} = 1 and *r _{d}* > 1. Further, all individuals with some, but not all mutations have fitness

If the mutation probability is small, we do not need to make specific assumptions on the mutation process. But when the mutation probability increases, we can no longer be certain that only a single mutation occurs during reproduction. For simplicity, we do not consider the possibility of backward mutations. Although back mutations are often relevant, especially to escape from evolutionary dead ends (DePristo et al., 2007), it is not straightforward to define the speed of evolution in a system with backward mutations. This is due to the fact that for sufficiently high mutation rates, fixation in the final state might never occur. Other definitions of the end state of the system become arbitrary to a certain extend. The probability *u _{m→m+k}* that the offspring of an individual with

$${u}_{m\to m+k}=\left(\begin{array}{c}d-m\\ k\end{array}\right)\phantom{\rule{0.2em}{0ex}}{\mu}^{k}{(1-\mu )}^{d-m-k}.$$

(1)

This equation is valid for the hypercube, where the order of mutations does not matter. Here,
$\left(\begin{array}{c}d-m\\ k\end{array}\right)$ is the number of different types of mutants with *k* additional mutations, *μ ^{k}* is the probability that mutations occur at

For small mutation probabilities, double mutations can be neglected. Since mutations occur rarely, we can calculate the average time until *d* mutations are fixed in the population analytically. Let us first address the evolutionary dynamics when mutants with fitness *r _{m}* are already present in a resident population of fitness

$${T}_{j}^{+}=\frac{{r}_{m}j}{{r}_{m}j+{r}_{w}(N-j)}\frac{N-j}{N}.$$

(2)

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

$${T}_{j}^{-}=\frac{{r}_{w}(N-j)}{{r}_{m}j+{r}_{w}(N-j)}\frac{j}{N}.$$

(3)

The probability that *k* mutants take over the entire population is given by (Nowak, 2006; Karlin and Taylor, 1975; Ewens, 2004; Crow and Kimura, 1970)

$${\phi}_{k}\left(\frac{{r}_{m}}{{r}_{w}}\right)=\frac{1+{\sum}_{i=1}^{k-1}{\prod}_{j=1}^{i}\frac{{T}_{j}^{-}}{{T}_{j}^{+}}}{1+{\sum}_{i=1}^{N-1}{\prod}_{j=1}^{i}\frac{{T}_{j}^{-}}{{T}_{j}^{+}}}=\frac{1-{\left(\frac{{r}_{w}}{{r}_{m}}\right)}^{k}}{1-{\left(\frac{{r}_{w}}{{r}_{m}}\right)}^{N}}.$$

(4)

If a mutant reaches fixation, the average number of generations this process takes is given by (Goel and Richter-Dyn, 1974; Antal and Scheuring, 2006)

$${\tau}_{\mathrm{fix}}\left(\frac{{r}_{m}}{{r}_{w}}\right)=\frac{1}{N}\sum _{k=1}^{N-1}\sum _{l=1}^{k}\frac{{\phi}_{l}}{{T}_{l}^{+}}\prod _{m=l+1}^{k}\frac{{T}_{m}^{-}}{{T}_{m}^{+}}.$$

(5)

For a neutral process with *r _{m}* =

The total time *τ* until a mutation reaches fixation in a population is the sum of the waiting time until a successful mutant occurs and the fixation time of the mutant *τ* = *τ*_{wait} + *τ*_{fix}. The waiting time is the inverse of the mutation rate divided by the probability that a particular mutant is successful,

$${\tau}_{\mathrm{wait}}\left(\frac{{r}_{m}}{{r}_{w}}\right)=\frac{1}{\mu N}\frac{1}{{\phi}_{1}\left(\frac{{r}_{m}}{{r}_{w}}\right)}.$$

(6)

For *μ* → 0, we have *τ*_{wait} → ∞, but *τ*_{fix} remains approximately constant. Thus, *τ* ≈ *τ*_{wait} for small mutation rates. In principle, we could calculate *τ*_{fix} in the presence of mutations. But since our approximation is only valid for small mutation rates, this will be a minor correction.

For *μ* *N*^{−2}, the population is homogeneous most of the time. Only occasionally, a mutant arises and reaches fixation or goes to extinction. The total time until *d* mutations are fixed in the population is the sum of the waiting times for the successful mutants plus the time of the *d* fixation events. For a single path with initial fitness 1, intermediate fitness *s* and final fitness *r*, we find for the total time *τ ^{S}*

$$\begin{array}{ll}{\tau}^{S}& ={\tau}_{\mathrm{wait}}(s)+(d-2){\tau}_{\mathrm{wait}}(1)+{\tau}_{\mathrm{wait}}(r/s)\\ & +{\tau}_{\mathrm{fix}}(s)+(d-2){\tau}_{\mathrm{fix}}(1)+{\tau}_{\mathrm{fix}}(r/s).\end{array}$$

(7)

For small *μ*, we have *τ _{fix}*

$${\tau}^{S}=\frac{1}{\mu}\left[\frac{1}{N{\phi}_{1}(s)}+d-2+\frac{1}{N{\phi}_{1}(r/s)}\right]$$

(8)

Consider now a “fitness valley”, in which the intermediate states have fitness *s* < 1, but the final state has fitness *r* > 1. To move from the fitness peak in the initial state to the fitness state in the final state, first a disadvantageous mutation has to be fixed in the population. Since
${\phi}_{1}(s<1)\ll \frac{1}{N}$, the waiting time of such an event is very long. The waiting time for the neutral mutations,
${\tau}_{\mathrm{wait}}(1)=\frac{1}{\mu}$ and the waiting time for a successful mutation, *τ*_{wait} (*r / s*) are significantly shorter. Thus, *τ ^{S}* is dominated by
$\frac{1}{\mu N{\phi}_{1}(s)}$ for

Fixation time for a single path (squares) and a hypercube (circles) with small mutation rates (*μ* *N*^{−2}) for different intermediate fitness values. Evolutionary dynamics is always faster in the hypercube. Solid lines show the analytical **...**

If the order of mutations is arbitrary, evolutionary dynamics occurs on a hypercube. In this case, the whole process will be faster, as we have *d*! possible paths instead of a single one. Now, the waiting times in the different states depend on the number of mutations that are still available. For the total time *τ ^{H}*, we obtain,

$$\begin{array}{ll}{\tau}^{H}& =\frac{1}{d}{\tau}_{\mathrm{wait}}(s)+\sum _{k=1}^{d-2}\frac{1}{d-k}{\tau}_{\mathrm{wait}}(1)+{\tau}_{\mathrm{wait}}\left(\frac{r}{s}\right)\\ & +{\tau}_{\mathrm{fix}}(s)+(d-2){\tau}_{\mathrm{fix}}(1)+{\tau}_{\mathrm{fix}}\left(\frac{r}{s}\right)\end{array}$$

(9)

Note that the time of the fixations alone is identical for the hypercube and the single path. Neglecting these fixation times for small *μ* (as *τ _{fix}*

$${\tau}^{H}=\frac{1}{\mu}\left[\frac{1}{\mathit{dN}{\phi}_{1}(s)}+\sum _{k=1}^{d-2}\frac{1}{d-k}+\frac{1}{N{\phi}_{1}(\frac{r}{s})}\right].$$

(10)

Since
$\frac{1}{dN{\phi}_{1}(s)}<\frac{1}{N{\phi}_{1}(s)}$ and
${\sum}_{k=1}^{d-2}\frac{1}{d-k}<{\sum}_{k=1}^{d-2}1=d-2$ it is obvious that *τ ^{H}* <

Next, we derive expressions for some interesting thresholds of the waiting times in the limit of small mutation rates. Since evolutionary dynamics is always faster if many paths are available, we now compare a fitness valley in which many paths are available to a single path in which the order of mutations is important, but fitness does not decrease in the course of evolution. The basic question we address here is, whether it is faster to cross a broad valley or a narrow ridge in fitness space. In other words, we compare *τ ^{S}* (

$$d-1+\frac{1}{N}\frac{1-\frac{1}{{r}^{N}}}{1-\frac{1}{r}}=\frac{1}{\mathit{dN}}\frac{1-\frac{1}{{s}^{N}}}{1-\frac{1}{s}}+\sum _{k=1}^{d-2}\frac{1}{d-k}+\frac{1}{N}\frac{1-{(\frac{s}{r})}^{N}}{1-\frac{s}{r}}$$

(11)

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

$$d(d-1)-\sum _{k=1}^{d-2}\frac{d}{d-k}=\frac{1}{N}\frac{1-\frac{1}{{s}^{N}}}{1-\frac{1}{s}}\approx \frac{{e}^{N(1-s)-1}}{N(1-s)},$$

(12)

where we used (1 − *x/N*)^{−N} → *e ^{x}* for large

The figure shows the threshold values for which evolution on a hypercube with fitness *s* < 1 in the intermediate states proceeds as fast as on a single neutral path with *s* = 1. Full lines show *N*(1 − *s*) based on the numerical solution of **...**

Next, we address the effect of the intermediate fitness *s*, which has an important influence on the cumulative waiting time *τ*. Fig. 2 shows how the waiting time decreases with increasing fitness in the intermediate states *s*. If *s* comes very close to the fitness in the final state, the waiting time increases. This increase is seen both in the single path and the hypercube. An increase in the intermediate state fitness will not always lead to a reduction in waiting times. Instead, the fixation times reach a minimum when the fitness growth is constant between any two consecutive states (Weinreich and Chao, 2005; Traulsen et al., 2007). For the hypercube, the fastest trajectory will be steeper than on a single path: At first, many mutations are available and a big fitness increase is not necessary. Later, fewer mutations are available and thus, the fitness should increase faster. The precise form of the trajectory will in this case depend on the number of mutations *d* and the population size *N*. We note that a similar reasoning can be applied to construct a fitness landscape that allows to cross a fitness valley fastest: Because the fastest trajectory has the same form regardless if a fitness peak is approached (*r* > 1) or a fitness minimum is approached (*r* < 1). Thus, the fastest way to cross a fitness valley is to descend to the minimum with exponentially decreasing fitness and to increase from the minimum again with exponentially increasing fitness.

Now, we turn to the effect of the intermediate fitness *s* on the individual waiting times. Eq. (8) and Eq. (10) both consist of three terms each. The first term denotes the time required to leave the initial state. The second term is the time spent in moving through all the intermediate states. This second term is independent of *s*, because the transitions are neutral. The last term denotes the time required to reach the ultimate state from the penultimate state. For small values of *s*, the probability to fixate the disadvantageous mutation is very small. Thus, the total time is dominated by the first term. When *s* is increased to a threshold value *s*_{1}, then the time for leaving the first state is identical to the waiting time in the intermediate states. For the hypercube, *s*_{1} is given by
$\frac{1}{d}{\tau}_{\mathrm{wait}}({s}_{1})={\sum}_{k=1}^{d-2}\frac{1}{d-k}{\tau}_{\mathrm{wait}}$ (1), which reduces to

$$\frac{1-\frac{1}{{s}_{1}^{N}}}{1-\frac{1}{{s}_{1}}}=\mathit{dN}\sum _{k=1}^{d-2}\frac{1}{d-k}.$$

(13)

This equation can be solved numerically for specific values of *N* and *d*. For the single path, the right hand side of this equation has to be replaced by *N*(*d* − 2). For *s* > *s*_{1}, the time to cross the intermediate states is larger than the waiting time in the first state. On the hypercube, we can define a second threshold for which the waiting time in the first state is the same as the time required to reach the final state from the penultimate state. This arises because the effective mutation rate in state 0 is *d* times larger than the effective mutation rate in state *d* − 1. The threshold *s*_{2} is given by
$\frac{1}{d}{\tau}_{\mathrm{wait}}({s}_{2})={\tau}_{\mathrm{wait}}(\frac{r}{{s}_{2}})$ or

$$\frac{1}{d}\frac{1-\frac{1}{{s}_{2}^{N}}}{1-\frac{1}{{s}_{2}^{1}}}=\frac{1-{(\frac{{s}_{2}}{r})}^{N}}{1-\frac{{s}_{2}}{r}}$$

(14)

Again,*s*_{2} has to be determined numerically. For a single path, the factor *d*^{−1} in Eq. (14) has to be dropped. Thus, the threshold *s*_{2} occurs for *s* > 1 and is simply given by
${s}_{2}=\sqrt{r}$.

The fixation time is also strongly influenced by the number of mutations *d*. A larger *d* increases the length of the path and usually also the fixation times. For the single path, this increase results only from the increase in the time required to cross the intermediate states, because the time for leaving the initial state and the time to reach the final state from the penultimate state are independent of *d*. The time required to reach the ultimate state from the penultimate state is also independent of *d* for the hypercube, but the time required to leave the initial state decreases with increasing *d*. This is because as *d* increases, there are more states available in the first error class and thus the effective rate of mutation out of the initial state increases. As for the single path, the time to cross the intermediate states increases with *d* in the hypercube. For the hypercube, this interplay can lead to a non-monotonic dependence of the fixation time on *d*. For example, for *N* = 100 and *s* = 0.95, the fixation time *τ ^{H}* decreases with

Increasing the fitness of the final state *r* increases the advantage of the final state over the intermediate states. This will result in the decrease in the time required for the population to make the last move. Increasing *r* has no effect on the time required to cross the intermediate states or the time required to move away from the initial state. As a result, those two times remain constant even as *r* increases, both in the single path and the hypercube.

The analytical approach is only valid as long as the mutation rate is small, *μ* *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 et al., 2004; Michor et al., 2004; Nowak et al., 2004; Beerenwinkel et al., 2007). 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 et al., 2004). As discussed above, for *μN*^{2} 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 *μN* > 1, at least one mutant is produced per generation. Thus, the probability of stochastic tunneling approaches 1. For *N*^{−2} < *μ* < *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.

The probability of tunneling across the hypercube (circles, blue) is larger than in the single path (squares, red) due to the higher effective mutation rate. The tunneling across the valley denoted by the filled symbols (*s* = 0.9) is always larger than **...**

The probability of tunneling across a neutral hypercube (circles) is always higher than the probability of tunneling across a neutral single path (squares). Here, the probability that the system tunnels across at least one state or one error class is **...**

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 and Kimura, 1970; Fisher, 1930; Muller, 1932; Gerrish and Lenski, 1998; Park and Krug, 2007). 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 and Lenski, 1998; Wilke, 2004; Jain and Krug, 2007).

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
$\left(\begin{array}{c}d\\ k\end{array}\right)=\frac{d!}{k!\left(d-k\right)!}$.

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 *μ* for the special case *d* = 2.

For *d* = 2, the probability of tunneling across the intermediate state is slightly higher in the single path (squares) than in the hypercube (open circles), shown for *N* = 100 here. This is because the effective mutation rate into the intermediate state **...**

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 *μ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}*(

$${\dot{x}}_{0}(t)=\frac{1}{N}\left[\frac{{x}_{0}}{\phi}{u}_{0\to 0}(1-{x}_{0})-(1-\frac{{x}_{0}}{\phi}{u}_{0\to 0}){x}_{0}\right].$$

(15)

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→k}, which is given above for both cases. The sum in *λ _{k}* is over all individuals with

$${\dot{x}}_{k}(t)=\frac{1}{N}[{\lambda}_{k}(1-{x}_{k})-(1-{\lambda}_{k}){x}_{k}],$$

(16)

where *k* = 0, …, *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 et al., 2007). 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 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 (Cooper et al., 2003), it is difficult to observe the consecutive fixation of neutral mutants. Consecutive fixation of advantageous mutants, however, is significantly faster. For example, while Table 1 reveals a fixation time of ~ 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 et al., 2007) would lead to a fixation time of ~ 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 et al. (2006) studied experimentally the point mutations in the *β*-lactamase gene of bacteria. *β* lactam antibiotics are commonly used, but the bacteria can develop resistance to the drugs. Five point mutations in a particular allele of the *β*-lactamase gene increases the resistance of the bacteria to cefotaxime by a factor of ~ 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 et al., 2007), 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 et al., 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 and McBroom, 1985). 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, 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.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.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Antal T, Scheuring I. Fixation of strategies for an evolutionary game in finite populations. Bull Math Biol. 2006;68:1923–1944. [PubMed]
- Armitage P, Doll R. The age distribution of cancer and a multi-stage theory of carcinogenesis. Br J Cancer. 1954;8:1–12. [PMC free article] [PubMed]
- Atwood KC, Schneider KL, Ryan FJ. Periodic selection in escherichia coli. Proc Natl Acad Sci USA. 1951;37:146–155. [PubMed]
- Beerenwinkel N, Antal T, Dingli D, Traulsen A, Kinzler KW, Velculescu VE, Vogelstein B, Nowak MA. Genetic progression and the waiting time to cancer. PLoS Comput Biol. 2007;3:e225. [PubMed]
- Chao L, McBroom SM. Evolution of transposable elements: an IS10 insertion increases fitness in Escherichia coli. Molecular Biology and Evolution. 1985;2:359–369. [PubMed]
- Cooper T, Rozen D, Lenski RE. Parallel changes in gene expression after 20,000 generations of evolution in escherichia coli. Proc Natl Acad Sci USA. 2003 Feb;100(3):1072–1077. [PubMed]
- Crow JF, Kimura M. Harper and Row. New York, NY: 1970. An introduction to population genetics theory.
- DePristo MA, Hartl DL, Weinreich DM. Mutational reversions during adaptive protein evolution. Mol Biol Evol. 2007;24:1608–1610. [PubMed]
- Eigen M, McCaskill J, Schuster P. The molecular quasi-species. Adv Chem Phys. 1989;75:149–263.
- Eigen M, Schuster P. The hypercycle. a.principle of natural self-organization. part a: Emergence of the hypercycle. Die Naturwissenschaften. 1977;64:541–565. [PubMed]
- Ewens WJ. Mathematical Population Genetics. Springer; New York: 2004.
- Fisher J. Multiple-mutation theory of carcinogenesis. Nature. 1959;181:651–652. [PubMed]
- Fisher RA. The Genetical Theory of Natural Selection. Clarendon 1930
- Frank SA. Evolutionary Dynamics of Cancer. Princeton Univ Press; 2007.
- Gerrish PJ, Lenski RE. The fate of competing beneficial mutations in a asexual population. Genetica. 1998;102/103:127–144. [PubMed]
- Gerstung M, Beerenwinkel N. Waiting time models of cancer progression. arXiv:0807.3638v1 2008
- Gillespie JH. Some properties of finite populations experiencing strong selection and weak mutation. The American Naturalist. 1983;121(5):691–708.
- Gillespie JH. Population Genetics : A Concise Guide. 2. The Johns Hopkins University Press; 2004.
- Goel N, Richter-Dyn N. Stochastic Models in Biology. Academic Press; New York: 1974.
- Iwasa Y, Michor F, Nowak MA. Stochastic tunnels in evolutionary dynamics. Genetics. 2004;166:1571–1579. [PubMed]
- Jain K, Krug J. Deterministic and stochastic regimes of asexual evolution on rugged fitness landscapes. Genetics. 2007;175:1275–1288. [PubMed]
- Jones S, Chen W-D, Parmigiani G, Diehl F, Beerenwinkel N, Antal T, Traulsen A, Nowak MA, Siegel C, Velculescu VE, Kinzler KW, Vogelstein B, Willis J, Markowitz SD. Comparative lesion sequencing provides insights into tumor evolution. Proc Natl Acad Sci USA. 2008a;105(11):4283–4288. [PubMed]
- Jones S, Zhang X, Parsons DW, Lin JC-H, Leary RJ, Angenendt P, Mankoo P, Carter H, Kamiyama H, Jimeno A, Hong S-M, Fu B, Lin M-T, Calhoun ES, Kamiyama M, Walter K, Nikolskaya T, Nikolsky Y, Hartigan J, Smith DR, Hidalgo M, Leach SD, Klein AP, Jaffee EM, Goggins M, Maitra A, Iacobuzio-Donahue C, Eshleman JR, Kern SE, Hruban RH, Karchin R, Papadopoulos N, Parmigiani G, Vogelstein B, Velculescu VE, Kinzler KW. Core signaling pathways in human pancreatic cancers revealed by global genomic analyses. Science. 2008b;321(5897):1801–1806. [PMC free article] [PubMed]
- Karlin S, Taylor HMA. A first course in stochastic processes. 2. Academic; London: 1975.
- Knudson AG. Mutation and cancer: Statistical study of retinoblastoma. Proc Natl Acad Sci USA. 1971;68:820–823. [PubMed]
- Lee TH, DSouza LM, Fox GE. Equally parsimonious pathways through an rna sequence space are not equally likely. J Mol Evol. 1997;45:278–284. [PubMed]
- Michor F, Iwasa Y, Nowak MA. Dynamics of cancer progression. Nature Reviews Cancer. 2004;4:197–205. [PubMed]
- Moran PAP. The statistical processes of evolutionary theory. Clarendon; Oxford: 1962.
- Muller HJ. Some genetic aspects of sex. The American Naturalist. 1932 Mar-Apr;66(703):118–138.
- Nordling C. A new theory on cancer inducing mechanism. Br J Cancer. 1953;7:68–72. [PMC free article] [PubMed]
- Nowak MA. What is a quasispecies? Trends in Ecology and Evolution. 1992;7:118–121. [PubMed]
- Nowak MA. Evolutionary Dynamics. Harvard University Press; Cambridge, MA: 2006.
- Nowak MA, Michor F, Komarova NL, Iwasa Y. Evolutionary dynamics of tumor suppressor gene inactivation. Proc Natl Acad Sci USA. 2004;101:10635–10638. [PubMed]
- Park S-C, Krug J. Clonal interference in large populations. Proc Natl Acad Sci USA. 2007;104(46):18135–18140. [PubMed]
- Poelwijk FJ, Kiviet DJ, Weinreich DM, Tans SJ. Empirical fitness landscales reveal accessible evolutionary paths. Nature. 2007;445:383–386. [PubMed]
- Press WH, Teukolsky SA, Vetterling WT, Flannery BP. The Art of Scientific Computing. 3. Cambridge University Press; New York: 2007. Numerical Recipes.
- Sjöblom T, Jones S, Wood L, Parsons D, Lin J, Barber T, Mandelker D, Leary R, Ptak J, Silliman N, Szabo S, Buckhaults P, Farrell C, Meeh P, Markowitz S, Willis J, Dawson D, Willson J, Gazdar A, Hartigan J, Wu L, Liu C, Parmigiani G, Park B, Bachman K, Papadopoulos N, Vogelstein B, Kinzler K, Velculescu V. The consensus coding sequences of human breast and colorectal cancers. Science. 2006;314:268–274. [PubMed]
- Traulsen A, Iwasa Y, Nowak MA. The fastest evolutionary trajectory. J Theor Biol. 2007;249(3):617–623. [PMC free article] [PubMed]
- van Nimwegen E, Crutchfield J. Metastable evolutionary dynamics: crossing fitness barriers or escaping via neutral paths? Bull Math Biol. 2000;62:799–848. [PubMed]
- Vogelstein B, Kinzler K. Cancer genes and the pathways they control. Nature Medicine. 2004;10:789–799. [PubMed]
- Weinreich D. The rank ordering of genotypic fitness values predicts genetic constraint on natural selection on landscapes lacking sign epistasis. Genetics. 2005;171:1397–1405. [PubMed]
- Weinreich D, Delaney N, DePristo M, Hartl D. Darwinian evolution can follow only very few mutational paths to fitter proteins. Science. 2006;312:111–114. [PubMed]
- Weinreich DM, Chao L. Rapid evolutionary escape by large populations from local fitness peaks is likely in nature. Evolution. 2005;59:1175–1182. [PubMed]
- Weinreich DM, Watson R, Chao L. Perspective: Sign epistasis and genetic constraint on evolutionary trajectories. Evolution. 2005;56(6):1165–1174. [PubMed]
- Wilke C. The speed of adaptation in large asexual populations. Genetics. 2004;167:2045–2053. [PubMed]
- Wilke C. Quasispecies theory in the context of population genetics. BMC Evolutionary Biology. 2005;5:44. [PMC free article] [PubMed]
- Wood LD, Parsons DW, Jones S, Lin J, Sjoblom T, Leary RJ, Shen D, Boca SM, Barber T, Ptak J, Silliman N, Szabo S, Dezso Z, Ustyanksky V, Nikolskaya T, Nikolsky Y, Karchin R, Wilson PA, Kaminker JS, Zhang Z, Croshaw R, Willis J, Dawson D, Shipitsin M, Willson JKV, Sukumar S, Polyak K, Park BH, Pethiyagoda CL, Pant PVK, Ballinger DG, Sparks AB, Hartigan J, Smith DR, Suh E, Papadopoulos N, Buckhaults P, Markowitz SD, Parmigiani G, Kinzler KW, Velculescu VE, Vogelstein B. The genomic landscapes of human breast and colorectal cancers. Science. 2007;318(5853):1108–1113. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |