Fixation probability of a nonmutator in a large population of asexual mutators

TL;DR

This paper calculates the fixation probability of a nonmutator in a large asexual mutator population, revealing complex dependencies on mutation rates and the influence of beneficial mutations, using a multitype branching process.

Contribution

It introduces a model to compute fixation probabilities of nonmutators considering mutation effects and beneficial mutations, providing new insights into mutation rate evolution.

Findings

Fixation probability is nonmonotonic with mutator mutation rate.

Beneficial mutations increase nonmutator fixation probability.

Fixation probability decreases exponentially at high mutation rates.

Abstract

In an adapted population of mutators in which most mutations are deleterious, a nonmutator that lowers the mutation rate is under indirect selection and can sweep to fixation. Using a multitype branching process, we calculate the fixation probability of a rare nonmutator in a large population of asexual mutators. We show that when beneficial mutations are absent, the fixation probability is a nonmonotonic function of the mutation rate of the mutator: it first increases sublinearly and then decreases exponentially. We also find that beneficial mutations can enhance the fixation probability of a nonmutator. Our analysis is relevant to an understanding of recent experiments in which a reduction in the mutation rates has been observed.

Equations

$$W(k)=(1-s{)}^k,$$

$$1-\pi (k,t)=\sum _{n=0}^{\infty }{\psi }_n(k){\left[\sum _{j}M(k\to j)(1-\pi (j,t+1))\right]}^n.$$

$$1-\pi (k,t)=\exp \left[-\frac{W(k)}{\overline{W}(t)}\sum _{j}M(k\to j)\pi (j,t+1)\right],$$

$$1-\pi (k)=\exp \left[-\frac{W(k)}{\overline{W}}\sum _{j}M(k\to j)\pi (j)\right].$$

$$\Pi =\sum _{k=0}^{\infty }p(k)\pi (k).$$

$$M(k\to k+i)=u_d{\delta }_{i,1}+(1-u_d){\delta }_{i,0},$$

$$1-{\pi }_d(k)=\exp \left[-\frac{(1-s{)}^k}{e^{-U_d}}\left(u_d{\pi }_d(k+1)+(1-u_d){\pi }_d(k)\right)\right].$$

$$k_{\max }=\lfloor \left(U_d-u_d\right)/s\rfloor .$$

$${\pi }_d(k)+\frac{{\pi }_d^2(k)}{2}\approx (1+U_d-sk)\left[u_d{\pi }_d(k+1)+(1-u_d){\pi }_d(k)\right]$$

$$\frac{{\pi }_d^2(k)}{2}\approx u_d{\pi }_d(k+1)+(U_d-sk-u_d){\pi }_d(k),0\le k\le k_{\max }.$$

$${\Pi }_d=2(U_d-u_d).$$

$$Q_d^2(k)-\left(1-\frac{u_d}{U_d}-\frac{sk}{U_d}\right)Q_d(k)-\frac{u_d}{U_d}{Q}_d(k+1)=0,$$

$$Q_d(k)=1-\frac{sk}{U_d},$$

$$Q_d(k)=1-\frac{sk}{U_d}+\frac{u_d}{U_d}\Delta (k).$$

$$Q_d(k)=\frac{{\pi }_d(k)}{2U_d}=1-\frac{sk}{U_d}-\frac{u_d}{U_d}\frac{s}{U_d-sk},$$

$${\Pi }_d=\sqrt{\frac{2s{U}_d}{\pi }}\left[1-\frac{u_d}{2U_d}\ln \left(\frac{2s{U}_d}{u_d^2}\right)\right].$$

$$Q_d(k)=\frac{{\pi }_d(k)}{2U_d}\approx 1-\frac{sk}{U_d}-\frac{u_d}{U_d}\left[1-{\left(\frac{u_d}{s}\right)}^{-2^{k+1-k_{\max }}}\right].$$

$${\Pi }_d$$

$$\approx$$

$${\Sigma }_1$$

$${\Sigma }_2$$

$$M(k\to k+i)=\{\begin{array}{l}{u}_d{\delta }_{i,1}+(1-u_d){\delta }_{i,0},k=0\\ u_d{\delta }_{i,1}+u_b{\delta }_{i,-1}+(1-u_b-u_d){\delta }_{i,0},k>0.\end{array}$$

$$1-\pi (k)=\{\begin{array}{l}\exp \left[-\frac{W(0)}{\overline{W}}\left((1-u_d)\pi (0)+u_d\pi (1)\right)\right],k=0\\ \exp \left[-\frac{W(k)}{\overline{W}}\left(u_b\pi (k-1)+(1-u_d-u_b)\pi (k)+u_d\pi (k+1)\right)\right],k>0,\end{array}$$

$$\frac{{\pi }^2(k)}{2}\approx s(\bar{k}-k)\pi (k),k\ge 0.$$

$$\pi (k)=2s(\bar{k}-k),0\le k\le \lfloor \bar{k}\rfloor .$$

$$\Pi =2s\sum _{k=0}^{\lfloor \bar{k}\rfloor }(\bar{k}-k)p(k),$$

$$\frac{\Pi }{{\Pi }_d}=\sqrt{\frac{{\sigma }^2}{U_d/s}}\approx 1-\frac{U_b}{2U_d},$$

$$\frac{dp(0,t)}{dt}$$

$$\frac{dp(k,t)}{dt}$$

$$p_d(k)=e^{-U_d/s}\frac{(U_d/s{)}^k}{k!}.$$

$$p(k)=\frac{(U_d/U_b{)}^{k/2}{J}_{k+\frac{(2-{\zeta }_0)\sqrt{U_b{U}_d}}{s}}\left(\frac{2\sqrt{U_b{U}_d}}{s}\right)}{{\sum }_{m=0}^{\infty }(U_d/U_b{)}^{m/2}{J}_{m+\frac{(2-{\zeta }_0)\sqrt{U_b{U}_d}}{s}}\left(\frac{2\sqrt{U_b{U}_d}}{s}\right)}.$$

$$J_{1+\frac{(2-\zeta )\sqrt{U_b{U}_d}}{s}}\left(\frac{2\sqrt{U_b{U}_d}}{s}\right)-\left(2-\sqrt{\frac{U_b}{U_d}}-\zeta \right)J_{\frac{(2-\zeta )\sqrt{U_b{U}_d}}{s}}\left(\frac{2\sqrt{U_b{U}_d}}{s}\right)=0.$$

Full text

††journal: Journal of Theoretical Biology

1 Introduction

Following the conclusion that mutation rates are subject to the action of evolutionary forces [40], there have been various experimental [5, 45, 39, 10, 30, 26, 41, 29, 46, 37] and theoretical [22, 25, 42, 44, 20, 32, 47, 38, 27, 8, 17, 19, 18, 11] works on the evolution of mutation rates. Many theoretical and empirical studies on adapting populations [36] have shown that the mutator alleles that elevate the mutation rates can reach a high frequency by generating beneficial mutations and hitchhiking with them [28].

However, once the population has adapted to an environment, due to high rate of production of deleterious mutations, the mutators experience a selective disadvantage and the nonmutator allele that lowers the mutation rate is favored due to indirect selection. Indeed, in a long term evolution experiment on E. coli, the frequency of mutators increased in three out of twelve replicate lines while the population was rapidly adapting [39]. But when the rate of fitness increase slowed down considerably, one of the mutator lines experienced a decrease in its mutation rate [46]. Several other experiments have also provided evidence for the rise in frequency of nonmutator allele in an adapted population [45, 30, 29, 37].

In this article, we are interested in a theoretical understanding of the evolution of mutation rates in adapted populations. In particular, using a multitype branching process [14, 33], we study the fixation probability of a nonmutator allele in a large asexual population of mutators that is moderately well adapted. In a recent work by us [19], this question was addressed when the nonmutator arises in the background of strong mutators whose mutation rate is ten to hundred fold higher than the nonmutator [39, 31]. However, as experiments show that the mutation rate decreases merely by a factor two to three in an adapted population [29, 46], here we undertake a more general investigation by allowing the nonmutator’s mutation rate to be comparable to that of the mutator. We also address how beneficial mutations in the mutator that increase its fitness affect the nonmutator fixation. Unlike in [19] where this question was studied in a limited parameter regime, here aided by an exact solution for the population frequency distribution that was obtained recently [16], we explore the parameter space completely.

The article is organised as follows: we define the model, and describe the simulation details and a multitype branching process in Sec. 2. The results for the fixation probability of a nonmutator are obtained in Sec. 3 followed by a discussion in Sec. 4.

2 Models and methods

2.1 Individual-based computer simulations

We consider an asexual population of mutators in which a mutation, irrespective of its location on the genome, changes fitness by a constant factor. Thus the fitness of an individual carrying $k\geq 0$ deleterious mutations (or, in the $k$th fitness class) is given by

[TABLE]

where $0<s<1$ is the selection coefficient. The population is of finite size $N$ and evolves via the standard Wright-Fisher dynamics [9] in which a parent in the fitness class $k$ is selected with a probability equal to $W(k)/N{\overline{W}(t)}$, where ${\overline{W}(t)}$ is the average fitness of the population at generation $t$. The selection step is followed by mutations; we employ a single-step mutation model in which mutations are allowed to occur in the neighboring fitness classes only. In an individual carrying $k>0$ unfavorable mutations, a deleterious mutation occurs at rate $U_{d}$ to fitness class $k+1$ and a beneficial one at rate $U_{b}<U_{d}$ to fitness class $k-1$. In the fittest individual, only deleterious mutations are allowed.

Motivated by a long-term evolution experiment on E. coli in which the nonmutator allele emerged in a mutator population when its fitness had almost saturated [39, 46], we allow the nonmutator to appear after the mutator population has attained a steady state [19, 18]. The invading nonmutator with deleterious and beneficial mutation rates $u_{d}$ and $u_{b}$, respectively, carrying $k$ unfavorable mutations arises in the mutator subpopulation in the $k$th fitness class with a probability equal to the stationary fraction of that subpopulation.

We measured the fixation probability of a single copy of nonmutator in a large population of mutators of strength $\lambda>1$ which is given by the ratio $U_{d}/u_{d}=U_{b}/u_{b}$. As we are interested in adapted populations in which beneficial mutations are rare, we first ignore the beneficial mutations completely (as discussed in Sec. 3.1) and then include beneficial mutations (see Sec. 3.2). In the former case, as the population size $N$ is finite, Muller’s ratchet [12] operates in the mutator population and there is no true steady state. For this reason, we simulated large enough populations in which the Muller’s ratchet clicks very slowly [15] and the mutator population is close to the stationary state of an infinitely large population. The fixation probability of a nonmutator was obtained using $10^{5}$ independent stochastic realizations of the mutator population; the results are shown in Figs. 1 and 2 when beneficial mutations are ignored and in Fig. 5 when they are taken into account.

2.2 Multitype branching process

In a finite population of mutators, a rare nonmutator allele - although beneficial due to indirect selection - can get lost because of stochastic fluctuations. But if it manages to survive random genetic drift, the nonmutator population can reach a frequency comparable to that of the mutators or even substitute them. Then it is interesting to ask: what is the probability that a rare beneficial allele arising in a large resident population does not go extinct? The branching process [14, 33] is tailor-made to answer precisely such questions and here we employ it to obtain an analytical understanding of our simulation results.

Let $1-\pi(k,t)$ denote the extinction probability that a nonmutator arising at generation $t$ in a mutator background with $k$ deleterious mutations is eventually lost. If such a nonmutator gives rise to $n$ offspring in the next generation with probability $\psi_{n}(k)$, then all the $n$ lineages must go extinct in order to contribute to the probability $1-\pi(k,t)$. Furthermore, if mutations are also allowed to occur in the nonmutator from fitness class $k$ to $j$ with a probability $M(k\rightarrow j)$, then summing over the number of offspring produced, we can write [21]

[TABLE]

For the Wright-Fisher process described in the last subsection, the offspring number distribution can be approximated by a Poisson distribution ($\psi_{n}(k)=e^{-\mu_{k}}\mu_{k}^{n}/n!$) with mean $\mu_{k}$ equal to the relative fitness ${W(k)}/{\overline{W}(t)}$. We then arrive at

[TABLE]

where we have used that $\sum_{j}M(k\rightarrow j)=1$. As discussed in Sec. 2.1, we assume that the nonmutator appears only after the mutator population has attained a steady state ($t\to\infty$). In this limit, the fixation probability becomes independent of time and in the following, we drop the time argument to denote the stationary state properties. We are thus required to solve the following nonlinear equation for $\pi(k)$,

[TABLE]

The solutions of the above equation give the fixation probability of a nonmutator in the fitness class $k$. But the probability that the nonmutator appears in the $k$th fitness class is given by the frequency $p(k)$ of the mutator population in that fitness class. Thus the total fixation probability obtained by summing over all the mutator backgrounds can be expressed as [21]

[TABLE]

In (4) and (5) above, we will use the deterministic results for the equilibrium frequency distribution $p(k)$ and the average fitness ${\overline{W}}$ of the mutator as we are working with large mutator populations in which the stochastic fluctuations can be neglected (see Appendix A1 for details). The mutation matrix $M(k\to j)$ is defined in (6) and (22) below. The recurrence equation (4) along with (5) are implemented using the software Wolfram Mathematica. The numerical results thus obtained are compared with those from stochastic simulations in Figs. 1 and 2, and we see a very good agreement. Therefore, in most of the following section, we will discuss results obtained using the multitype branching process.

3 Results

3.1 Only deleterious mutations

Since the beneficial mutation rates are much smaller than their deleterious counterparts [34], as a first approximation, we set the beneficial mutation rates $U_{b}$ and $u_{b}$ equal to zero and denote the quantities of interest with a subscript $d$.

The frequency distribution of the mutator $p_{d}(k)$ is known to be Poisson-distributed with mean $U_{d}/s$ [23, 12] which gives the mean fitness ${\overline{W}}_{d}=\sum_{k=0}^{\infty}W(k)p_{d}(k)=e^{-U_{d}}$ (also see Appendix A1). Furthermore, as for the mutators, we assume that the mutations in the nonmutator occur in the neighboring fitness classes only so that

[TABLE]

with $u_{d}\ll 1$. From (4), we thus obtain

[TABLE]

Note that $\pi_{d}(k)=0$ for all $k\geq 0$ is an exact solution of the above equation. However, as explained below, a nontrivial solution for the fixation probability exists if the number of mutations carried by the nonmutator are small enough.

If a nonmutator arising with $k$ detrimental mutations escapes random genetic drift and displaces the resident mutator population, the steady state fitness of the resulting population is given by $e^{-u_{d}}(1-s)^{k}$ [12]. The additional factor $(1-s)^{k}$ reflects the fact that every individual in such a population carries at least $k$ deleterious mutations. Then the maximum number of mutations that the nonmutator can carry so as to have a selective advantage over the mutator population is determined by the condition $e^{-u_{d}}(1-s)^{k}>e^{-U_{d}}$ [21], or $k<k_{\max}$ where

[TABLE]

Here, $\lfloor x\rfloor$ denotes the maximum integer less than or equal to $x$. Equation (7) along with the boundary condition $\pi_{d}(k_{\max}+1)=0$ can be solved numerically in a straightforward manner for $k\leq k_{\max}$ [21].

However, to obtain an analytical solution, we approximate (7) using the fact that all the variables ($u_{d},U_{d},s,\pi_{d}$) are smaller than one and furthermore, the product $sk<1$ due to (8). Taking logarithm on both sides of (7) and using the expansion $\ln(1-x)$ $\approx$ $-x-x^{2}/2$ for $x\ll 1$, we obtain

[TABLE]

which can be further simplified to yield the following nonlinear recursion equation,

[TABLE]

Some remarks are in order: on dividing both sides of the above equation by $(2s)^{2}$, we first note that the scaled fixation probability $\pi_{d}(k)/(2s)$ is a function of two scaled mutation rates, viz., $u_{d}/s$ and $U_{d}/s$. The three cases considered in the following subsections are classified according to whether these ratios lie below or above one. (Of course, one can also choose to scale the variables by one of the mutation rates.) Second, for a given integer $k_{\max}$, the model parameters lie in the range $k_{\max}<(U_{d}-u_{d})/s<k_{\max}+1$. Thus, when $k_{\max}$ is large, we may ignore the fact that it is an integer and write $k_{\max}\approx(U_{d}-u_{d})/s$. For small integer $k_{\max}$, some simple cases are worked out in Appendix A2.

3.1.1 Case I: $u_{d}<U_{d}<s$

In this case, as the selective effect of a deleterious mutation is large, any nonmutator carrying nonzero deleterious mutations gets eliminated from the population and therefore $\pi_{d}(1)=0$ (this conclusion also follows from (8)). Using this boundary condition in (10), we obtain a quadratic equation for $\pi_{d}(0)$ whose nonzero root is given by $2(U_{d}-u_{d})$. Furthermore, from (A3), since $p_{d}(0)\approx 1$ when $U_{d}<s$, we find that the total fixation probability defined in (5) is given by [27, 19]

[TABLE]

Thus the total fixation probability is simply twice the fitness advantage $e^{-u_{d}}-e^{-U_{d}}\approx U_{d}-u_{d}$ conferred by the nonmutator [13, 27, 19].

3.1.2 Case II: $u_{d}<s<U_{d}$

As $U_{d}$ is the largest variable in this parameter regime, on dividing both sides of (10) by $(2U_{d})^{2}$, we can rewrite it as

[TABLE]

where $Q_{d}(k)=\pi_{d}(k)/2U_{d}$. Since $u_{d}/U_{d}$ is the smallest parameter here, we first ignore the terms containing $u_{d}/U_{d}$ in the above equation and immediately find that the fixation probability decays linearly with the fitness class,

[TABLE]

which shows that a nonmutator has a low chance of fixation if it arises in a mutator background with close to $U_{d}/s$ mutations. However, as the mutator frequency is Poisson-distributed with mean $U_{d}/s$ (see (A3)), the nonmutator is most likely to occur in fitness classes in the neighborhood of $U_{d}/s$ and therefore such low-fitness classes can still contribute to the total fixation probability. Furthermore, due to the form of (13) above, the total fixation probability (5) can be interpreted as the average (positive) deviation from the mean which, by virtue of (A3), is given by $\sqrt{U_{d}/s}$ and therefore we expect $\Pi_{d}\propto\sqrt{sU_{d}}$ (also see (16) below).

The above discussion applies when the mutator is strong, i.e., its mutation rate is much larger than that of the nonmutator ($U_{d}\gg u_{d}$) [19]. However, for weak mutators for which $u_{d}/U_{d}$ is not negligible, corrections to the above behavior can be found by expanding the probability $\pi_{d}(k)$ in a power series in the small parameter $u_{d}/U_{d}$:

[TABLE]

Substituting the above expression in (10) and retaining terms to linear order in $\Delta(k)$, we obtain

[TABLE]

which shows that the nonmutator’s chance of fixation is lowered when the mutator is weaker.

The total fixation probability (5) obtained by summing over the mutator backgrounds is calculated in Appendix A3 and we find that

[TABLE]

When the mutation rate $u_{d}=0$, the second term in the bracket on the RHS vanishes and we recover the result in [19]. The reduction in the fixation probability is, however, not appreciable for moderately strong mutators. For $s=10^{-4},U_{d}=10s$ and $\lambda=50,100,200$, using (4), we find the fixation probability to be $0.96,0.98,0.99$ times the fixation probability when $\lambda\to\infty$, respectively. The corresponding numbers obtained using (16) are $0.939,0.962,0.977$ which overestimate the reduction by $1-2\%$.

3.1.3 Case III: $s<u_{d}<U_{d}$

We again consider (12) as both $s/U_{d}$ and $u_{d}/U_{d}$ are small here. Although $s/U_{d}$ is the smallest parameter, we can not neglect the term $sk/U_{d}$ in (12) as it increases with the fitness class $k$. To tackle this case, as described in Appendix A4, we first obtain an approximate solution for $Q_{d}(k)$ for large fitness classes and then use its properties to arrive at an approximate solution for all the fitness classes which is given by

[TABLE]

In the above expression, when $k\ll k_{\max}$, the last term on the RHS can be neglected and the fixation probability decreases linearly with the fitness class as in the last subsection. But for $k\lesssim k_{\max}$, the decay is faster than an exponential. Equation (17) for $\pi_{d}(k)$ is plotted in the inset of Figs. 3 and  4 for strong and weak mutator, respectively, against the results obtained by solving (4) numerically, and we find a very good agreement.

When the mutator is strong ($U_{d}\gg u_{d}$), the term containing $u_{d}/U_{d}$ on the RHS of (17) can be neglected and we find that the fixation probability decays linearly with the fitness class as indeed supported by the inset in Fig. 3. Then, as shown in Appendix A4, the total fixation probability (5) is given by

[TABLE]

When $U_{d}\ll 2s\lambda^{2}$ or equivalently, $U_{d}\gg u_{d}^{2}/(2s)$, the above equation shows that the fixation probability increases as $\sqrt{U_{d}}$ (as in the last subsection where $u_{d}<s<U_{d}$). However, in the opposite parameter regime, $\Pi_{d}$ decreases exponentially with $U_{d}$. This can be understood as follows: Because of (A3) for the mutator frequency distribution $p_{d}(k)$, the nonmutator is most likely to arise in fitness classes with $U_{d}/s-\sqrt{U_{d}/s}<k<U_{d}/s$ mutations. However, the fixation probability $\pi_{d}(k)$ is zero in this interval if $k_{\max}<(U_{d}/s)-\sqrt{U_{d}/s}$ which, on using $k_{\max}\approx(U_{d}-u_{d})/s$, implies that when $U_{d}>s\lambda^{2}$, the chances of nonmutator fixation are considerably reduced. Figure 3 shows that expression (19) matches well with the exact numerical calculations up to an additive constant as we have neglected the contribution from the $u_{d}$-dependent terms in (17).

When the mutator is weak ($U_{d}\gtrsim u_{d}$), as the inset of Fig. 4 shows, the nonlinear decay of the fixation probability for large fitness classes can not be ignored. In Appendix A4, the total fixation probability is calculated for large $k_{\max}$ and we find that $\Pi_{d}=\Sigma_{1}+\Sigma_{2}$, where

[TABLE]

and the proportionality constants are function of $\lambda$. We thus find that $\Pi_{d}\sim e^{-\frac{U_{d}}{2s\lambda^{2}}}$ for $U_{d}/s\gg 2\lambda^{2}$ in agreement with the data in Fig. 4.

3.2 Both deleterious and beneficial mutations

We now consider the case when the beneficial mutation rates for mutator and nonmutator denoted by $U_{b}$ and $u_{b}$, respectively, are nonzero. It was first observed in [19] that beneficial mutations in the mutator can increase the fixation probability of a nonmutator. This counterintuitive effect was shown in a limited parameter regime where $U_{b}\ll s$ and $u_{b}=0$. Our purpose here is to explore the validity of this result for a broader set of parameters.

3.2.1 When mutation rates $u_{b},u_{d}$ are nonzero

When both deleterious and beneficial mutations occur in the nonmutator, the mutation rate $M(k\to j)$ from fitness class $k$ to $j$ in the nonmutator is given by

[TABLE]

Using this in (4) for the fixation probability, we obtain

[TABLE]

where the average fitness of the mutator population ${\overline{W}}=e^{-s{\bar{k}}}$ with ${\bar{k}}$ being the average number of deleterious mutations. For nonzero $U_{b}$, the mean ${\bar{k}}$ given by (A6) is smaller than $U_{d}/s$, as one would intuitively expect.

When beneficial mutations are ignored, by virtue of $u_{b}=0$, equation (23) for the probability $\pi(k_{\max})$ closes (i.e., it does not involve the fixation probability in any other fitness class) and therefore $\pi(k_{\max})$ can be determined. This result then allows one to numerically calculate the fixation probability in lower fitness classes ($k<k_{\max}$) [21]. However, for nonzero $u_{b}$, equation (23) shows that for $k>0$, the probability $\pi(k)$ is coupled to the fixation probability in both the neighboring fitness classes. Thus if the fixation probability is zero beyond a fitness class $k_{\max}$, a calculation of $\pi(k_{\max})$ requires the knowledge of $\pi(k_{\max}-1)$. For this reason, it is difficult to analyse (23) even numerically and we have not been able to come up with an efficient method to do so.

However, using the simulation method described in Sec. 2.1, we have studied this case for $U_{d}/s>1$ and our results are shown in Fig. 5. We find that for large enough $U_{b}/s$, the total fixation probability $\Pi$ is smaller than $\Pi_{d}$; this behavior is expected as the mutational load carried by the mutators is reduced due to beneficial mutations resulting in a decrease in the fixation probability of nonmutator. However, when $U_{b}/s$ is small, the fixation probability $\Pi>\Pi_{d}$ which is surprising. Figure 5 also shows that the ratio $\Pi/\Pi_{d}$ is larger for weaker mutators and thus a lower bound on the ratio is obtained when $\lambda\to\infty$ (or, equivalently $u_{b}=u_{d}=0$).

3.2.2 When mutation rates $u_{b},u_{d}$ are zero

For $u_{d}=u_{b}=0$, as the surprising effect of beneficial mutations described above survives [19] and the recursion equations for the probability $\pi(k)$ are amenable to analysis (see below), we now study this case using the multitype branching process.

Assuming that the mutation rates and the selection coefficient are small, on proceeding in a manner similar to that in the last subsection, (23) yields

[TABLE]

The nontrivial solution of the above equation is given by

[TABLE]

The total fixation probability is then obtained as

[TABLE]

where $p(k)$ is given by (A4). In [19], the frequency $p(k)$ was obtained by solving (A2) numerically; however, an exact expression for $p(k)$ that was obtained recently [16] allows us to extend our previous results (see (27) below).

When $U_{d}/s\ll 1$, the average number of deleterious mutations is smaller than one ($\lfloor\bar{k}\rfloor=0$) and we have $\Pi\approx 2s{\bar{k}}$ which, on using (A7), yields the relative fixation probability $\Pi/\Pi_{d}\approx 1-(U_{b}/s)$, in agreement with (11) of [19].

When $U_{d}/s\gg 1$, we need to consider the parameter regimes $U_{b}/s\ll 1$ and $\gg 1$ separately. In the former case, the fixation probability $\Pi$ is larger than $\Pi_{d}$ as supported by a perturbation theory in $U_{b}/s$ which yields $\Pi/\Pi_{d}=1+(U_{b}/s)$ [19]. But as Fig. 6 shows, this quantitative dependence agrees with the numerical calculations of (26) in a narrow parameter range ($U_{b}/s\lesssim 0.01$) while the relative fixation probability stays above one for $U_{b}/s<(U_{d}/s)^{-1}$. For larger $U_{b}/s$, however, the relative fixation probability is below unity. For $U_{d}/s\gg 1,U_{b}/s\gg 1$, the mutator frequency distribution is well approximated by a Gaussian function [16] with mean ${\bar{k}}$ and variance $\sigma^{2}$ given by (A6) and (A8), respectively. Using these results in (26), we find that

[TABLE]

which matches well with the numerical data shown in the inset of Fig. 6. The nonmonotonic change in $\Pi/\Pi_{d}$ is due to the discreteness of $k_{\max}$, as explained in Appendix A2.

4 Discussion

In this article, we have studied how the fixation probability of a rare nonmutator that arises in a large adapted population of asexual mutators depends on the mutation rates and the selection coefficient. Motivated by a long-term experiment on E. coli in which two nonmutator lineages arose in the mutator population when the fitness growth had slowed down considerably [46], we have calculated the fixation probability in the stationary state of the mutator population.

Figure 1 shows that the fixation probability $\Pi_{d}$ (that ignores beneficial mutations) increases with selection coefficient $s$. This trend is consistent with the intuitive expectation that when the selective cost of a deleterious mutation is high, the mutation rate should be low. However, for sufficiently large $s$ where all mutations can be treated as lethal [27], the probability $\Pi_{d}$ saturates to a maximum value (see Sec. 3.1.1). In this parameter regime, only one type of nonmutator - the one without any deleterious mutations - has a nonzero chance of fixation and the classic single-locus theory applies [13, 27]. But for smaller selection coefficients, a nonmutator can carry many deleterious mutations and a multilocus analysis is required as has been done here using a multitype branching process [33]. It is known that the fixation of a beneficial mutant is impeded due to interference from deleterious mutations when the loci are tightly linked [21]. Viewing the nonmutator allele as a beneficial mutant under indirect selection, the reduction in the fixation probability in the multiloci scenario follows from this result. The above discussion is also qualitatively consistent with the expectation that sexual populations in which loci are weakly linked are more likely to have low mutation rates [20, 43, 35].

As seen in Fig. 2, the fixation probability of the nonmutator also increases with the mutator strength $\lambda=U_{d}/u_{d}$ where $U_{d}(u_{d})$ is the deleterious mutation rate of the mutator (nonmutator); this is because a population of strong mutators would carry more deleterious mutations and hence more likely to get lost. To put it differently, a nonmutator with higher mutation rate has a lower chance of fixation because it also accrues deleterious mutations which weaken its advantage over the mutators. These expectations are indeed borne out by our analyses in Secs. 3.1.2 and 3.1.3.

The above discussion suggests that a nonmutator with a mutation rate comparable to that of the mutator [29, 46] is unlikely to fix. However, from our detailed analysis in Sec. 3.1.3, we arrive at a novel and important conclusion that a nonmutator is most likely to fix when the mutator population has a mutation rate $U_{d}\sim s\lambda^{2}$ (also, see Figs. 3 and 4). In a long term evolution experiment [46], the mutation rate of the ${\it E.coli}$ population carrying mutator allele decreased merely by a factor two at large times. As this event occurred in two independent lineages, using the data for $U_{d}$ and $\lambda$ given in Table 2 of [46] in the above criteria for most probable mutation rate, we find the selection coefficient $s\sim 10^{-2}$ so that the selection to deleterious mutations is moderate in this experiment. We also note that for a wide range of selection coefficients $s=10^{-4}-10^{-2}$ and $U_{d}=0.03$ [46], the mutation rate is most likely to decrease by a factor $20$ or less; the above discussion thus provides an explanation for a small decrease in the mutation rates. The nonmonotonic behavior of the fixation probability with mutator’s mutation rate, shown in Figs. 3 and 4, may be understood as follows: if $U_{d}$ is too small, the nonmutator does not offer significant advantage and therefore has a small chance of fixation. On the other hand, if $U_{d}$ is too large, most individuals in the mutator population carry a large number of deleterious mutations. In this case, a nonmutator is favored only if it arises in mutator backgrounds with few deleterious mutations. But such mutator subpopulations are rare when $U_{d}$ is large thereby rendering the nonmutator fixation unlikely (also, see the discussion below (19)).

A factor studied here that can decrease the chances of nonmutator fixation is the occurrence of beneficial mutations in the mutator population as a result of which it carries lesser mutational load. However, we point out that when the mutator is weak, the beneficial mutations can enhance the fixation probability, see Fig. 5. As explained in [19], this counterintuitive effect arises because of the competition between two opposing factors: while beneficial mutations improve the fitness of the mutator population and hence adversely affect the fixation probability of the nonmutator, they also present a fitter mutator background for the nonmutator to arise thus augmenting its chance of fixation. This factor may also contribute to the fixation of the nonmutator in the experiment of [46] but the relevant beneficial mutation rates do not seem to be available.

Throughout this discussion, we have assumed that the mutator population is large enough so that random genetic drift may be ignored. A negative correlation between effective population size and mutation rates has been found in recent studies [41, 26] and rationalised using a simple argument [27, 19]; however, a more detailed and general analysis is needed for a better understanding of this relationship. For example, all the studies mentioned here except [18] have examined the process of mutation rate reduction on a single-peaked, multiplicative fitness landscape; more complex fitness landscapes [7] could be used to understand the evolution of mutation rates. Moreover, the physiological costs [22, 24, 6, 20, 3] associated with the decline in mutation rates can be included in a future study to gain an insight into its effect on the mutation rate evolution.

Acknowledgements

The authors thank Sona John for discussions during the early stages of this work. A. James thanks CSIR for the funding.

Appendix A1 Deterministic model for frequency distribution of mutators

For small selection coefficients, we can study the deterministic evolution of the mutator frequency distribution $p(k,t)$ in continuous time [4]. The fitness $W(k)$ defined in the discrete time model (see (1)) and the fitness $w(k)$ in the continuous time model (discussed below) are related through $W(k)=e^{w(k)}$ which gives $w(k)\approx-sk$. Then the mutator frequency distribution obeys the following deterministic equations [16],

[TABLE]

where $\bar{k}=\sum_{k=0}^{\infty}kp(k,t)$ is the average number of deleterious mutations in the mutator population. The last term on the RHS of the above equations is the selection term while the other terms correspond to mutations in the neighboring fitness classes. In the stationary state where the LHS is zero, we will denote the frequency by $p(k)$.

In the absence of beneficial mutations ($U_{b}=0$), the above equations simplify to yield the Poisson-distributed mutator frequency with mean $U_{d}/s$ [23, 12]:

[TABLE]

Thus the average fitness ${\overline{W}}_{d}=e^{-s\sum_{k=0}^{\infty}kp_{d}(k)}=e^{-U_{d}}$.

When the beneficial mutations also occur, the stationary state population fraction is given by [16]

[TABLE]

Here, $J_{n}(z)$ is the Bessel function of first kind with order $n$ and argument $z$ [1] and $\zeta_{0}$ is the smallest root of the following equation:

[TABLE]

For $U_{b}=0$, the result (A4) above reduces to (A3) [16].

For completeness, below we quote the results obtained in [16] that are pertinent to the discussion here. When beneficial mutations are allowed, as one would expect, the average number of deleterious mutations is smaller than $U_{d}/s$ and given by

[TABLE]

For nonzero $U_{b}$, the variance of the number of deleterious mutations is larger than the mean and given by

[TABLE]

Appendix A2 When the integer $k_{\max}$ is small

Here we consider some simple cases where the fixation probability can be found exactly using the recursion equation (10).

$\underline{k_{\max}=0}$: In this case, the difference in the mutation rates lie in the range $0<U_{d}-u_{d}<s$. As $\pi_{d}(1)=0$, we immediately have $\pi_{d}(0)=2(U_{d}-u_{d})$ and therefore

[TABLE]

where, for brevity, we have defined

[TABLE]

From the above expression, we find that for a given $\lambda$, the total fixation probability first increases and then decreases as a function of $\mu$ when $0<\mu<\lambda/(\lambda-1)$.

$\underline{k_{\max}=1}$: Using $\pi_{d}(2)=0$ in (10), we find the nonnegative solutions for the fixation probability to be

[TABLE]

As a result, the total fixation probability is given by

[TABLE]

and is a nonmonotonic function of the ratio $\mu$ when $\lambda/(\lambda-1)<\mu<2\lambda/(\lambda-1)$.

From the above simple examples and the inset of Fig. 2, we conclude that when $k_{\max}$ is small, the total fixation probability increases and decreases several times as $U_{d}/s$ is varied. However, as Fig. 2 also shows, this effect is unimportant for large $k_{\max}$.

Appendix A3 Fixation probability when $u_{d}<s<U_{d}$

Using (A3) and (15) in the sum (5), we find that the total fixation probability is given by

[TABLE]

We first note that for $u_{d}=0$, the sum in the above equation can be done exactly [19] and is given by

[TABLE]

We verify that the results (A11) and (A15) for $k_{\max}=0,1$, respectively, are reproduced from the above equation when $u_{d}=0$.

To estimate the fixation probability $\Pi_{d}$, we first approximate the Poisson distribution $p_{d}(k)$ with mean and variance $U_{d}/s$ by a Gaussian distribution with these properties [19]. On approximating the sum in (A16) by an integral, we get

[TABLE]

where $\mu=U_{d}/s$ and $\nu=u_{d}/s$ as defined in (A12). For small $\nu$ and large $\mu$, on carrying out the above integrals, we obtain

[TABLE]

The second term in the bracket on the RHS is obtained on using that the exponential integral $E_{1}(z)=\int_{z}^{\infty}dt~{}e^{-t}/t\approx-\ln z$ for small $z$ (see 5.1.1 and 5.1.11, [1]). We remark that although (A16) is linear in the small parameter $\nu$, the correction term is nonlinear due to the logarithmically-diverging second integral in (A19).

Appendix A4 Fixation probability when $s<u_{d}<U_{d}$

We first study the behavior of the fixation probability for fitness classes close to $k_{\max}$. Using $Q_{d}(k_{\max}+1)=0$ and $k_{\max}\approx(U_{d}-u_{d})/s$ for large $k_{\max}$ (see the discussion after (10)) in equation (12), we get $Q_{d}(k_{\max})=0$ which further yields $Q_{d}(k_{\max}-1)=s/U_{d}$. To obtain simple expressions for other fitness classes, we note that for fitness classes close to the boundary $k_{\max}$, the coefficient of $Q_{d}(k)$ in (12) can be neglected yielding a simpler recursion equation,

[TABLE]

which can be easily solved subject to the boundary condition at $k=k_{\max}-1$ and we get

[TABLE]

Note that the above fixation probability is of a double exponential form ($e^{-e^{-x}}$) and therefore decays faster than an exponential as the fitness class $k$ increases towards $k_{\max}$. Furthermore, for large $k_{\max}-k$, it saturates to $\lambda^{-1}$.

To find the behavior for small $k$, we now write the complete solution $Q_{d}(k)$ as

[TABLE]

Substituting this in (12), we get

[TABLE]

The expression (A26) is exact for $Q^{<}_{d}(k)$. But as it is nonlinear, we find an approximate solution for $Q^{<}_{d}(k)$ using the fact that $Q^{>}_{d}(k)$ saturates quickly to $\lambda^{-1}$ to arrive at (A27). If we now set (A27) to zero, we get

[TABLE]

The above solution shows that the error committed in ignoring (A27) is of the order $su_{d}/U_{d}^{2}$ which is negligible in the parameter regime under consideration. Using (A23) and (A28) in (A24), we finally obtain the fixation probability in (17).

To find the total fixation probability, we need to perform the following sums,

[TABLE]

where $p_{d}(k)$ is the Poisson distribution given by (A3) and $\mu,\nu$ are defined in (A12).

Strong mutator ($\lambda\gg 1$): Neglecting the mutation rate of the nonmutator in (A29) when $U_{d}\gg u_{d}$, we find that

[TABLE]

Alternatively, in (A30), we can approximate the Poisson distribution by a Gaussian and the sum by an integral (as in Appendix A3) to get the above result.

Weak mutator ($\lambda\gtrsim 1$): The first sum on the RHS of the above equation can be expressed as ((26.4.21), [1])

[TABLE]

where we have used that when $\nu\gg 1$, the integrand in (A35) is appreciable for $t\ll 1$. Solving the integral in the last equation, we obtain

[TABLE]

where $\textrm{erfc}(x)=(2/\sqrt{\pi})\int_{x}^{\infty}dt~{}e^{-t^{2}}$ is the complementary error function and we have used the asymptotic expansion of error function to obtain the last expression ((7.1.23), [1]).

The second sum on the RHS of (A29) can be estimated via a saddle-point method [2]. Denoting the logarithm of the summand by $f(k)$ and expanding it about its turning point at $k_{*}$ where $f^{\prime}(k_{*})=0$, we have

[TABLE]

where prime denotes a derivative with respect to $k$. In the above equation, $k_{*}$ is a solution of the following equation,

[TABLE]

and the second derivative is given by

[TABLE]

Taking logarithms on both sides of (A43), we find that

[TABLE]

which shows that $k_{*}$ is close to $k_{\max}$. This result furthermore yields $f^{\prime\prime}(k_{*})\sim\lambda^{-1}$. Since the validity of the saddle-point method requires $f^{\prime\prime}(k_{*})$ to be large [2], our approximation is good when $\lambda$ is not too large. We thus have

[TABLE]

Note that in the last equation, the factors in the bracket are a function of $\lambda$ and therefore the sum in (A46) is proportional to $\sqrt{\mu}e^{-\frac{\mu}{2\lambda^{2}}}$ for fixed $\lambda$.