Multimutation Leaps in the Equilibrium Regime.

Let us assume (binary) genomes of length L (in the context of this analysis, L should be construed as the number of evolutionarily relevant sites, such as codons in protein-coding genes, rather than the total number of sites), the probability of single mutation μ << 1 per site per round of replication (generation), and constant effective population size N_e >> 1. Then, the transition probability from sequence i to sequence j is (equation 3.11 in ref. 15)qij=μhij(1−μ)L−hij,[1]where hij is the Hamming distance (number of different sites between the 2 sequences). The number of sequences separated by the distance h is equal to the number of ways h sites can be selected from L, that is,Nh=L!(L−h)!h!≈Lhh!,[2]where the last, approximate expression is valid under the assumption that L >> 1 and L >> h (h can be of the order of 1).

Assuming also µ << 1, we obtain a typical combinatorial probability of leaps over the distance h:Q(h)∼Nhq(h)=Ph(Lμ)≡(Lμ)he−Lμh!,[3]which is a Poisson distribution with the expectation Lμ.

In steady state, the probability of fixation of the state i is proportional to exp(−νxi) whereν=Ne−1(Moran process)[4]ν=2(Ne−1) (haploid Wright-Fisher process)ν=2Ne−1 (diploid Wright-Fisher process)

and xi=−ln⁡fi where fi is the fitness of the genotype i [xi is analogous to energy in the Boltzmann distribution within the analogy between population genetics and statistical physics (16)]. For other demographic structures and assumptions on the mutation process, the relationship between fixation probability can quantitatively differ while retaining the same form. In particular, for a population that produces offspring by binary division (fission), ν≈4Ne (17, 18).

Then, the rate of the occurrence and fixation of the transition i→j is (15)Wij=qijν(xi−xj)exp[ν(xi−xj)]−1.[5]The distribution function of the fitness differential Δij=xi−xj has to be specified (hereafter, we refer to x as fitness, omitting logarithm for brevity). We analyze first the case without epistasis, that is, with additive fitness effects of individual mutations:Δ(h)=y1+y2+…+yh,[6]where yi are independent random variables with the distribution functions Gj(yj). Then, the distribution function of the fitness difference isρh(Δ)=∏j∫dyjGj(yj)δ(∑jyj−Δ)=∫−∞∞dk2πe−ikΔ∏j∫dyjeikyjGj(yj),[7]which is obtained by using the standard Fourier transformation of the delta function.

Now, let us specify the distribution of the fitness effects of mutations Gi(yi), assuming an exponential dependency of the probability of a mutation on its fitness effect, separately for beneficial and deleterious mutations:Pi(yi)={Die−ϵiyi,yi>0Dirieϵiyi,yi<0,[8]where Di is the normalization factor, ri is the ratio of the probabilities of beneficial and deleterious mutations, and ϵi is the inverse of the characteristic fitness difference for a single mutation (discussed below). For simplicity, we assume here the same decay rates for the probability density of the fitness effects of beneficial and deleterious mutations. Empirical data on the distributions of fitness effects of mutations (19, 20) clearly indicate that ri≪1. From the normalization condition,Di=ϵi1+ri≈ϵi.[9]Note that the mean of the fitness difference (selection coefficient) when the distribution of the fitness effects is given by 8 is|si|=∫dyiGi(yi)yi≈1ϵi.[10]For simplicity, we start with an assumption that the values of Di and ri are independent of i. For the model 8∫−∞∞dyG(y)eiky=iD(1k+iϵ−rk−iϵ).[11]Then, from Eq. 5, the fixation rate of an h-mutation leap is equal toφ(h)=∫−∞∞dΔνΔeνΔ−1ρh(Δ).[12]Substituting 11 into 7, we obtainρh(Δ)=−ih+1ϵh(h−1)!(1+r)hdh−1dkh−1[(1−rk+iϵk−iϵ)he−ikΔ]k=−iϵ, Δ> 0ρh(Δ)=ih+1ϵh(h−1)!(1+r)hdh−1dkh−1[(k−iϵk+iϵ−r)he−ikΔ]k=iϵ, Δ < 0.[13]Consider first the case r = 0 (all mutations are deleterious). Then, ρh(Δ<0)=0. For Δ > 0, that is, decrease of the fitness, we haveρh(Δ)=Δh−1ϵh(h−1)!e−ϵΔ.[14]Then, the fixation rate 12 of an h-mutation leap is equal toφ(h)=zh(h−1)!∫0∞dtthe−ztet−1=hzhζ(h+1,z+1),[15]where z=ϵ/ν and ζ(x,y) is the Hurwitz zeta function ζ(x,y)=∑k=0∞1(k+y)x. Therefore, the rate of fixation for leaps of the length h is equal to W(h)=Ph(Lμ)φ(h).

In one extreme, if z≫1 (ν|s|≪1, neutral landscape), φ(h)≈1 and mutations are fixed at the rate they occur. In the opposite extreme case of strong negative selection (z≪1,ν|s|≫1), φ(h)≈hzhζ(h+1) where ζ(x) is the Riemann zeta function. For a rough estimate, ζ(h+1) can be replaced by 1, and then, W(h)≈Lμe−LμPh−1(Lμz). In this case, the maximum of W(h) is reached at h=Lμz≅Lμ/ν|s|, which gives a nonnegligible fraction of multimutation leaps (h>1) among the fixed mutations only for Lμ≥ν|s|. However, in this case, the value of W(h) at this maximum is exponentially small because e−Lμ<e−ν|s|. Therefore, in the regime of strong selection against deleterious mutations and at high mutations rates (Lμ≥ν|s|), multiple mutations actually dominate the mutational landscape, but their fixation rate is extremely low. Qualitatively, this conclusion seems obvious, but we now obtain the quantitative criteria for what constitutes “strong selection.” We find that, even for ν|s|∼10, the rate of multimutation leaps (h=4) can be nonnegligible (>10⁻⁴ per generation; Fig. 2A) at the optimal Lμ values, whereas for ν|s|∼100, any leaps with h>1 are unfeasible (Fig. 2B).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Rates of leaps on a landscape dominated by deleterious mutations. Rates of transitions are plotted against the per-genome mutation rate (Lμ) and the leap length for different strengths of selection (A: ν|s| = 10 and B: ν|s| = 100). Contour lines indicate orders of magnitude and start from the rate of 10⁻⁵ leaps per generation.

Under a more realistic model, all values of ϵi (the inverse of the fitness effect of a mutation) are different. For Δ>0 and r=0 (no beneficial mutations), using Eq. 13, we getρh(Δ)=∫−∞∞dk2πe−ikΔ∏jiϵjiϵj+k.[16]For example, in Kimura’s neutral evolution model (21), ϵi is a binary random variable that takes a value of ∞ (|si|=0, neutral mutation), with the probability f, and a value of 0 (|si|=∞, lethal mutation), with the probability 1−f. Then, ρh(Δ)=fhδ(Δ), φ(h)=fh and Lμ is replaced with Lfμ in Eq. 3, a trivial replacement of the total genome length L with the length of the part of the genome where mutations are allowed, Lf. Accordingly, W(h)=Ph(Lμf), and multimutation leaps become relevant for Lμf≥1.

Let us now estimate the probability of leaps with beneficial mutations (Δ<0). Assuming r≪1 (rare beneficial mutations), Eq. 13 takes the formρh(Δ)≈hrϵ2h−1e−ϵ|Δ|[17]and the fixation rate of a leap including beneficial mutations isφ(h)=hr2h−1zζ(2,z).[18]If z≫1 (weak positive selection), zζ(2,z)≈1, so that the role of beneficial mutations is negligible. If z≪1 (strong positive selection),φ(h)=hr2h−11z.[19]Comparing Eq. 19 with the result for Δ>0 (Eq. 18), one can see that, in this case, beneficial mutations are predominant among the fixed mutations ifr>(2ϵ/ν)h+1.[20]In this regime, multimutation leaps (h>4) occur at nonnegligible rates under sufficiently high (but not excessive) mutation rates (Fig. 3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Rates of leaps on a landscape combining beneficial and deleterious mutations. Rates of leaps are plotted against the per-genome mutation rate (Lμ) and the leap length for different strengths of selection (A and C: ν|s| = 10; B and D: ν|s| = 100) and for different frequencies of beneficial mutations (A and B: r = 10⁻⁴; C and D: r = 10⁻³). Contour lines indicate orders of magnitude and start from the rate of 10⁻⁵ leaps per generation.

The model considered above assumes independent effects of different mutations (no epistasis, “ideal gas of mutations” model). Now, let us take into account epistasis. In the case of strong epistasis, effects of combinations of different mutations are increasingly strong, diverse, and, effectively, unpredictable, resulting in a rugged fitness landscape (22). In the limit of epistasis strength and unpredictability, epistasis creates numerous highly beneficial combinations that, once they occur, are highly likely to be fixed, and a far greater number of highly deleterious combinations that are immediately lethal. Due to the effective randomness of genetic interactions, we consider the resulting landscape as essentially random for h>1, with the frequency of the beneficial combinations r independent of h. In this case, the effective number of fixed leaps with h>1 is simplyW(h)=Q(h)f∼Ph(Lμ)r.[21]If all single mutations (h=1) are deleterious (ν|s|≫1), their rate of fixation (Eq. 15) can be approximated by W(1)=P1(Lμ)1ν|s|, whereas for all leaps of the length h>1, effective number of fixed leaps is W(h>1)=(1−P0(Lμ)−P1(Lμ))f. Therefore, the condition for W(h>1)>W(1) isr>Lμe−Lμ(1−(Lμ+1)e−Lμ)1ν|s|.[22]In the high-mutation regime (Lμ≫1), multiple mutations occur orders of magnitude more frequently than single mutations, overwhelming the difference of scale between r and 1/ν|s|, and making multimutation leaps much more likely. In the low-mutation regime (Lμ≪1), the balance between single and multiple mutations tends to 2(1Lμ−1)≈2Lμ and the condition for the dominance of multimutation leaps becomes r>2Lμ1ν|s|. Around the Eigen threshold (Lμ≈|s|) (23), the condition corresponds to r>2ν|s|2, that is, the frequency of beneficial multimutation combinations should be unrealistically high to sustain evolution by multimutation leaps.

A Nonequilibrium Model of Stress-Induced Mutagenesis.

The analysis presented above suggests that the necessary condition for fixation of multimutational leaps is the high-mutation regime. At low mutation rates (Lμ≪1), multimutation (h>1) events occur too rarely to be fixed in realistic settings even if the frequency of beneficial combinations among them is reasonably high. However, in the high-mutation regime (Lμ≫1), the above analysis is problematic for 2 reasons. First, the expression for the fixation rate (Eq. 5) is technically valid only for the case when the new mutation is either fixed or lost before the emergence of the next one, which implies Lμ<1/ν≪1. Second, under any realistic model of the fitness landscape, most mutations should be deleterious. Thus, Lμ>1 implies that most of the progeny carries 1 or more mutations, and therefore suffers from these deleterious effects. Under these conditions, the assumption of constant N_e is unrealistic, because the size of such a population will decrease under the mutational load, down to an eventual crash.

The complete analysis of the behavior of a variable-size population under the high-mutation regime and strong mutational effects is currently beyond the state of the art. Therefore, here we analyze a simplified model of the short-term behavior of a (microbial) population after the onset of stress-induced mutagenesis (Lμ≫1).

Consider a microbial population consisting of N0 individuals. Under typical conditions, the population is in an equilibrium, so that approximately N0/2 individuals survive the average generation span and produce N1≈N0 progeny by division (here we consider simple asexual division as the progeny-generating process whereby each surviving individual produces 2 offspring; other demographic models can be accommodated without loss of generality). The typical mutation rate is low (Lμ0≈1/Ne≈1/N0, according to refs. 24 and 25), so the population can be considered homogeneous. Upon the onset of unfavorable conditions, the survival rate of the wild-type individuals drops to fw≪1/2 and the mutation rate in the stressed individuals increases such that Lμ>1.

If fw is not too small (fwN0≫1), the immediate wild-type survivors produce 2fwN0 first-generation progeny. With the expected number of mutations per descendant being Lμ, the distribution of the number of mutations in the progeny is given by the Poisson distribution with the expected number of mutants with h mutations of 2fwPh(Lμ)N0.

Let us consider a mutation landscape that is dominated by deleterious mutations with strong sign epistasis. All single mutations are deleterious, so the survival of their carriers over the generation time is f1≪fw. An overwhelming majority of multimutation combinations have even stronger negative effects, so for h>1,fh≪f1. Some small fraction rh of these combinations, however, is strongly beneficial in the new conditions, conferring to their carriers the survival rate of ∼1/2.

What would the rh(h) function look like? Intuitively, rh(h) should decay to 0 at large h, or at least not grow, as it is overwhelmingly likely that a sufficiently large set of mutations would contain a subset that it unconditionally lethal. Here, for simplicity, we consider a general form of rh(h) that equals 0 for h=1 and monotonically decays with h from r2 at an arbitrary rate.

If the deleterious effect of mutations is strong enough (0≈fh≈f1≪fw), then the only plausible source of beneficial mutants is the population of wild-type individuals (neither single mutants nor multiple mutants that do not carry the beneficial combinations survive to the next generation). The population of the wild-type individuals decays exponentially through both the diminished survival and through mutations, reaching fw(2fwP0(Lμ))k−1N0 at the k-th generation after the onset of the unfavorable conditions. Ignoring stochastic fluctuations, the total number of wild-type individuals that survive until the population collapse can be estimated asNw∞≈fwN0/(1−2fwP0(Lμ)),[23]which is approximately equal to fwN0 if fwP0(Lμ)≪1/2.

Over the combined lifetimes of the surviving wild-type individuals, the expected number of beneficial mutants isE(NB∞)≈2Nw∞∑h=2∞(rhPh(Lμ))=2Nw∞∑h=2∞(rhe−Lμ(Lμ)hh!),[24]which depends on the genome-wide mutation rate Lμ and the shape of the rh(h) function.

Let us first consider the 2 extreme cases of rh(h). In the limit of a completely flat function (rh=r2 for all h>2), Eq. 23 gives E(NB∞)≈2Nw∞r2(1−P0(Lμ)−P1(Lμ)). This function asymptotically reaches the value of 2Nw∞r2 with Lμ→∞. In the other extreme of a rapidly decaying rh(h), that is, rh=0 for h>2, Eq. 23 gives E(NB∞)≈2Nw∞r2P2(Lμ). This function reaches its maximum at Lμ=2 with E(NB∞)≈4Nw∞r2e−2.

It can be shown that the estimates for all other monotonically decaying rh(h) functions reach their maxima at finite values of Lμ with4Nw∞r2e−2≤E(NB∞)≤2Nw∞r2.[25]Indeed, let us consider first the simplest model rh(h)=r2ξh−2 (0<ξ<1). Then, Eq. 24 takes the formE(NB∞)≈2Nw∞r2ξ−2e−Lμ(eLμξ−Lμξ−1).[26]As a function of Lμ, the quantityφ(Lμ)=e−Lμ(eLμξ−Lμξ−1)[27]reaches the maximum at Lμ∗=x/ξ, where x is the solution of the equationex−1x=11−ξ[28]with the value at the maximumφ(Lμ∗)=(1−ξ)1−ξξxξ(1−ξ+x)−1/ξ.[29]For ξ≪1 (rapid decay of rh), Lμ∗=2 and E(NB∞)≈4Nw∞r2e−2. In the opposite limit of slowly decaying rh, ξ→1, Lμ∗≈ln11−ξ and φ(Lμ∗)≈1.

For a general slowly decaying function rh(h), one can find thatLμ∗≈ln|rh|drhdh|h≈Lμ∗|.[30]Importantly, even in this case, the optimal mutation rate Lμ∗ increases only logarithmically with the decay rate; furthermore, the optimum value is notably robust to changes in rh(h) (Fig. 4).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Abundance of beneficial multimutation combinations depending on the mutation rate. Abundance of beneficial multimutation combinations, φ(Lμ), given by Eq. 24, relative to r2. (A) rh(h)=r2(h−1)−α with α=0 (blue), α=0.5 (orange), α=1 (green), and α=2 (red). (B) rh(h)=r2ξh−2 de with ξ=1 (blue), ξ=0.9 (orange), ξ=0.5 (green), and ξ=0.25 (red).

The approximate condition for population survival, E(NB∞)>2, can be derived from Eqs. 23 and 25 and is bounded from below byr2>e22N0fw[31]at the optimal value of Lμ.