The dynamics of alternative pathways to compensatory substitution

Abstract

The role of epistatic interactions among loci is a central question in evolutionary biology and is increasingly relevant in the genomic age. While the population genetics of compensatory substitution have received considerable attention, most studies have focused on the case when natural selection is very strong against deleterious intermediates. In the biologically-plausible scenario of weak to moderate selection there exist two alternate pathways for compensatory substitution. In one pathway, a deleterious mutation becomes fixed prior to occurrence of the compensatory mutation. In the other, the two loci are simultaneously polymorphic. The rates of compensatory substitution along these two pathways and their relative probabilities are functions of the population size, selection strength, mutation rate, and recombination rate. In this paper these rates and path probabilities are derived analytically and verified using population genetic simulations. The expected time durations of these two paths are similar when selection is moderate, but not when selection is weak. The effect of recombination on the dynamics of the substitution process are explored using simulation. Using the derived rates, a phylogenetic substitution model of the compensatory evolution process is presented that could be used for inference of population genetic parameters from interspecific data.

Background

Complex interactions among loci play an important role in many evolutionary processes. They have arisen in the context of adaptation as Sewall Wright's shifting balance theory [1, 2] and in speciation as Dobzhansky-Muller incompatibilities [3, 4]. Such interactions have been shown to play a role in the intramolecular evolution within proteins [5, 6] and RNAs [7, 8], as well as the intermolecular co-evolution between different proteins [9] and the regulation of gene expression [10]. An understanding of the the dynamics of such dependent evolutionary processes is therefore applicable for a wide range of fundamental questions in evolutionary biology.

Often these problems are cast in the framework of a fitness landscape, with peaks corresponding to genotypes of high fitness and valleys to genotypes of low fitness [11, 12]. When the landscape is very flat and all genotypes of similar fitness, different loci may appear to evolve independently. But when the landscape is more rugged and certain combinations of alleles are deleterious, the method and rate at which an entire population can cross from one peak to another become the central issues.

If a haploid genome possesses a combination of alleles at different loci that are associated with high fitness, a mutation at any locus is likely to result in a suboptimal combination of alleles and consequently a haplotype of low fitness. However, a second mutation at another locus may restore the combination of alleles to one of high fitness. An entire population moving from one fitness peak to another in this manner is called compensatory substitution.

Theoretical population geneticists have used simple but powerful models, typically consisting of two loci with two alleles at each locus, to discern the rate at which such shifts can occur [13–21]. The expected time for these shifts is a function of the population size, the mutation rate, the strength of selection against intermediates, and the probability of recombination between loci [13]. With only four allelic combinations, the frequency of these haplotypes in a population becomes a difficult three-dimensional diffusion problem [17, 18]. The approach has uniformly been to simplify the dimensionality of the problem in order to obtain analytical results, with simulations frequently used to verify results and explore problems for which analytic solutions could not be found [13–21].

Kimura demonstrated that population genetic theory could allow for the occurrence of compensatory substitutions, using the simplest possible model of compensatory evolution with two linked loci and two alleles at each locus [13]. Subsequent efforts have extended Kimura's results to models of greater complexity by allowing for: different fitnesses for the two "peak" haplotypes [14, 15, 22], different fitnesses for the two deleterious intermediate haplotypes [14, 15, 17, 19], different mutation rates [14, 17, 19], reversible mutations [18, 19, 21], diploid populations [20, 21], and expanding to a two locus four allele model [18]. These mutation/selection models are depicted in additional file 1.

Each of these studies made simplifying assumptions about the strength of selection and mutation, and in order to compare their results it is useful to know in what ranges of parameter space their results may be applicable. The figures from these papers and the simulations used by the authors are indicators of the range of parameter space in which these authors consider their approaches reasonable. Figure 1 summarizes the approximate ranges of parameter space explored by each of these papers, based on their simulations and figures. This focuses on mutation and selection, both scaled by population size, although most of these papers explored recombination as well.

Figure 1

The parameter ranges explored by previous work on compensatory evolution. θ = 4Nµ where N is the diploid population size, µ is the per-locus mutation rate and s is the selection coefficient against deleterious intermediates. All of these studies also explored recombination. In black: Kimura [13]. In grey: Iizuka and Takefu [14]. In transparent red box: Stephan [17]. In pink: Michalakis and Slatkin [15]. In purple: Higgs [18]. In light blue: Innan and Stephan [19].

While most of this theoretical work assumed that selection (and often mutation) was very strong (but see [19]), there are many cases in which selection against intermediates is not expected to be strong. A high frequency of deleterious intermediate states was observed among homologous stem pairs of RNA between species [23], taken as evidence that selection on RNA intermediates could not be as strong as previously supposed. The role effective population size has played in the evolution of tRNAs was examined [24] and it was concluded that natural selection on paired sites has been weak to moderate. It has also been observed in RNA intramolecular interactions that certain non-canonical intermediate configurations (GU or UG) are more stable than others, and that compensatory substitutions tend to favor paths through less-deleterious intermediates [12]. These studies underscore the need to understand the dynamics of the compensatory substitution process when selection is moderate to weak.

There are two pathways to compensatory substitution. When selection is very weak there is a substantial chance (particularly if the population size is small) that a deleterious mutation can become fixed in the population before a second (compensatory) mutation arises and goes to fixation. This will be referred to as a "Type 1" event (Figure 2A). However, when selection is very strong the deleterious haplotype cannot become fixed in the population. In this case, compensatory substitution can only occur if the second mutation arises before the first mutation is lost due to selection, and the new high-fitness double mutant may then drift to fixation. This kind of two-at-a-time substitution will be called a "Type 2" event (Figure 2B). When selection is of intermediate strength both of these substitution pathways are plausible (Figure 2C).

Figure 2

Illustrations of compensatory substitution pathways. The frequency of the two high fitness (blues) and low fitness (reds) haplotypes in a population are shown over time. (a) Type 1. Under weak selection the deleterious intermediate can readily fix in the population. (b) Type 2. Under strong selection the deleterious intermediates cannot rise above low frequency, and compensatory substitution must await the second mutation arising while the first is still polymorphic. (c) Under intermediate selection both Type 1 and Type 2 events are possible. Shown is a Type 2 event. Although the deleterious intermediate aB does become fixed, the eventual path to fixation of the double mutant ab begins with AB once again fixed and does not have the deleterious intermediate being fixed.

The dynamics of the compensatory substitution process for both Type 1 and Type 2 events is the focus of this work. The rates of these alternate events are derived, representing the first derivation of the direct double-substitution rate under weak selection and reversible mutation. The interactions among selection strength, population size, mutation, recombination, and their effect on the probability of the alternative compensatory substitution pathways are explored using simulation. This work also shows that the expected time to compensatory substitution may depend upon the fixation pathway. These results yield insights into how fitness valleys can be crossed when selection is more moderate.

Methods

Theory

Substitution Model

Consider a population of N diploid (or 2N haploid) individuals and a two-locus, two-allele mutation model as shown in Figure 3A. The first locus consists of alleles A and a, the second of alleles B and b. Let µ be the bidirectional mutation rate between alleles A and a, and let this same rate correspond to mutations between B and b. Under a haploid (genic) selection mechanism, the fitnesses of the four possible haplotypes AB, aB, Ab, and ab are 1, 1 - s, 1 - s, and 1 respectively. This situation corresponds to the case of symmetric neutral compensatory evolution [13] but with reversible mutation [18]. For this derivation assume that the loci are tightly linked, such that recombination between loci does not occur.

Figure 3

Graphical model representations. (A) The mutation-selection model used. The four haplotypes and their respective fitnesses are depicted within circles, with the mutations rates labeling the arrows. (B) Continuous-time Markov model of the compensatory substitution process. Rates between states are shown on arrows. Each state represents the first fixation event for the given haplotype. (C) The discrete-time Markov chain, or jump chain, embedded within the continuous-time Markov chain. The relative probabilities of moving between fixed states are shown on arrows. Deleterious states are considered together for simplicity, and ab is considered an absorbing state to match the simulation procedure.

Define a four-state continuous-time Markov chain of the substitution process, as shown in Figure 3B, in which each state represents the novel fixation of the population for each of these four haplotypes. The notion of novel fixation events is an important one: a population fixed for AB, becoming polymorphic then losing the polymorphism to become again fixed for AB does not represent a novel fixation event and no change of state will have been observed. In this model the population may move between two fixed states if (1) the haplotypes differ at only one of the two loci, or (2) both haplotypes are of high fitness.

Direct movement between states Ab and aB is disallowed, as such a move is highly unlikely unless intermediates are only very slightly deleterious. This is because if both of these states are of low fitness, the process must actually avoid a fitness peak while making two substitutions. If only single mutation events are allowed, then one of the fitter haplotypes must be created in the process, and is likely to fix via selection. If double mutations are allowed, these will occur at rate µ² and must then fix by drift, and it is far more likely that a single mutant of higher fitness will arise before this happens.

Rates into deleterious states

Consider a population fixed for the AB haplotype. Following the formulation of [25], the rate at which new aB haplotypes arise that are destined to go to fixation in this population is the product of the rate at which new aB mutants arise in the population and the probability that such a new mutation becomes fixed in the population. The rate at which aB haplotypes arise in a population of 2N AB haplotypes, or equivalently the probability that a new aB mutant will arise in such a population in a single generation, is 2Nµ, where µ is the per-locus mutation rate. Recalling that the aB haplotype has fitness 1 - s, the probability of fixation of such a new mutation is given by $\frac{1-e^{2s}}{1-e^{4Ns}}$[26]. Consequently the rate from AB to aB is

$$r_1=2N\mu \frac{1-e^{2s}}{1-e^{4Ns}}.$$

(1)

If the assumption is made that mutation rates are uniform across loci and between all alleles, then Equation 1 will hold for leaving any state with fitness 1 {AB, ab} and entering any state with fitness 1 - s {aB, Ab}.

The probability of fixation of a deleterious haplotype does not depend upon whether or not another fit haplotype arises in the process. Consider a population of AB individuals in which an aB haplotype has arisen. Before it becomes fixed, suppose a mutation occurs creating a new ab haplotype from one of the copies of aB. Because the AB and ab haplotypes are of equal fitness, the emergence of the new haplotype affects the probability of fixation of the intermediate exactly as would a reverse mutation. Alternatively, the introduction of a new Ab haplotype via mutation from AB will serve to increase slightly the probability of fixation of the aB haplotype by lowering the overall population fitness, but this effect is expected to be slight unless the mutation rate is very high.

Rates out of deleterious states

Consider a population fixed for the aB haplotype. If the population fitness is renormalized to be 1, then a new ab haplotype arising in such a population will have fitness $1+t=\frac{1}{1-s}$, and will have selection coefficient $t=\frac{s}{1-s}$ where t > 0 implies selection in favor of this haplotype. In this case the probability of fixation of such a new haplotype given that it has arisen is $\frac{1-e^{-2t}}{1-e^{-4Nt}}$, and such haplotypes arise by mutation at rate 2Nµ per generation, yielding an overall rate of

$$r_2=2N\mu \frac{1-e^{-2t}}{1-e^{-4Nt}}.$$

(2)

As before, under symmetric mutation this rate r₂ will hold for leaving states with fitness 1 - s {aB, Ab} and entering any state with fitness 1 {AB, ab}.

Rates of double substitution

For Equations 1 and 2 the assumption was made that new haplotypes arise in populations that are fixed for a particular state. It is possible however that a new haplotype could arise while the population is still polymorphic, and the identities and frequencies of the existing haplotypes could vary greatly. This will be considered when deriving the rate of entering the ab state.

With four possible haplotypes, obtaining formulae for the change in frequency of these haplotypes becomes a three-dimensional diffusion problem and is quite challenging. However, in the case in which some haplotype fitnesses are equivalent the problem may be simplified. By grouping states with equal fitness the dimensionality of the problem can be reduced [19]. Let x be the combined frequency of the two deleterious haplotypes in the population (each with fitness 1 - s), and y = 1 - x be the combined frequency of the remaining haplotypes (each with fitness 1). In this case, the problem reduces to a one-dimensional diffusion process on x. Sewall Wright [1] derived the equilibrium distribution of a mutant allele class resulting from the actions of reversible mutation and selection in a finite population as

$$\varphi \left(x\right)=C{e}^{4Nsx}{x}^{2\theta -1}{\left(1-x\right)}^{2\theta -1}$$

(3)

where θ = 4Nµ and C is a normalizing constant such that ${\int }_0^1\varphi \left(x\right)dx=1$. Equation 3 differs slightly from Wright's original because there are two loci and mutations may occur at either locus. This provides the distribution of the frequency of deleterious intermediate haplotypes {Ab, aB}.

Consider a population lacking the ab haplotype. If the mutation rate is sufficiently low that the possibility of a double mutant ab arising from an AB haplotype in a single generation can be ignored, then ab haplotypes can only arise from Ab or aB haplotypes via single mutation. Therefore, conditional on the deleterious frequency x, the probability that a new ab haplotype enters the population this generation is 2Nµx.

Under the assumption that the mutation rate is very small relative to the population size, the fixation probability of the y group can be approximated by

$$p_y=\frac{1-e^{-4Ns{y}^{\prime }}}{1-e^{-4Ns}}$$

(4)

where $y^{\prime }=y+\frac{1}{2N}$ is the frequency of the y group after a new ab mutant has entered the population [26]. The probability of fixation of the new ab haplotype is

$$p_{ab}=\mathbb{P}\left(ab\mathsf{\text{ fixes}}|AB\mathsf{\text{ or }}ab\mathsf{\text{ fixes}}\right)\mathbb{P}\left(AB\mathsf{\text{ or }}ab\mathsf{\text{ fixes}}\right)=\frac{p_y}{2N{y}^{\prime }}$$

(5)

since AB and ab have equal fitness. Taken together, and integrating over all possible values of the deleterious intermediate frequency, the average probability α in each generation that a new ab haplotype will arise that is destined to go to fixation is

$$\begin{array}{l}\alpha ={\displaystyle {\int }_0^1\mathbb{P}\left(ab\text{ arises}|x\right)\mathbb{P}\left(ab\text{ fixes}|x,ab\text{ arises}\right)\mathbb{P}\left(x\right)dx}\\ =2N\mu {\displaystyle {\int }_0^1{p}_{ab}x\varphi \left(x\right)}dx.n\end{array}$$

(6)

Equations 4 - 6, from [19] but included here that this manuscript may be self-contained, give the rate α at which new successful ab haplotypes enter the population. By marginalizing over the stationary frequencies of deleterious intermediates at equilibrium, this key insight [19] implicitly considers both pathways by which the new ab haplotype can fix (with or without a deleterious intermediate ever fixing in the population). Previous work on rates of compensatory substitution only considered the case in which the deleterious intermediate could never fix because selection was too strong [13, 14, 18].

Considering the Markov substitution model of Figure 3B, another way to describe Equation 6 is as the total rate into state ab

$$\begin{array}{l}\alpha ={\displaystyle \sum _{k\ne ab}{\pi }_k{r}_k\to ab}\\ =({\pi }_{aB}+{\pi }_{Ab})r_2+{\pi }_{AB}{r}_3\end{array}$$

(7)

where π _k is the stationary probability of being in state k and r₂ is given by Equation 2. So the total rate into state ab is the sum of rates into ab from each other state, weighted by the stationary probability of being in that particular state. This assumes that the detailed-balance conditions are met, meaning that the process is at stationarity. This is not such a strong assumption, as the process reaches something very close to stationarity quickly [19]. Solving for r₃ yields

$$r_3=\frac{\alpha -\left({\pi }_{aB}+{\pi }_{Ab}\right)r_2}{{\pi }_{AB}}.$$

(8)

Because the model is symmetric, α also describes the total rate into state AB, and therefore

$$r_4=\frac{\alpha -\left({\pi }_{aB}+{\pi }_{Ab}\right)r_2}{{\pi }_{ab}}.$$

(9)

The expected stationary frequency of the fixed states can be derived from Equation 3, which provides the distribution of the frequency of the pooled deleterious intermediate haplotypes (and conversely the pooled high-fitness haplotypes as well). Alternatively these frequencies can be estimated from simulations. Because the model is symmetric, π _AB = π _ab and π _Ab = π _aB .

Equations 8 and 9 represent the rates of direct double-substitution, without first fixing the deleterious intermediate. This is the first derivation of a direct double-substitution rate when selection is weak.

A discrete-time Markov model

Embedded within the continuous-time Markov model of Figure 3B describing the rates of moving between states is a discrete-time Markov chain representing the transition probabilities between states (Figure 3C). Consider, as before, only novel fixation events such that each state is the population becoming newly fixed for that state. The amount of time required for each fixation event to occur is not considered here, merely the state changes themselves. This model will assist in understanding some of the properties of the compensatory substitution process.

For simplicity, the deleterious haplotypes aB and Ab, both of relative fitness 1-s, are considered together as a single state representing the population being fixed for one of the deleterious haplotypes. Because both the AB and ab haplotypes have fitness 1 and because mutation is symmetric, the probabilities of moving from state aB/Ab to either AB or ab are both $\frac{1}{2}$. Since each simulation replicate ends when the population becomes fixed for the ab haplotype, ab is considered an absorbing state; a model of the true process differs only in that ab is not absorbing and would be completely symmetric.

In this Markov model there is some probability β of moving from state AB directly to state ab and some probability 1 - β of moving to one of the deleterious intermediate states. When β = 0 the compensatory substitution can only occur via the fixation of the intermediate, and when β = 1 the intermediate can never become fixed.

Simulations

The forward evolution of a population of 2N haploid individuals is simulated using Wright-Fisher sampling. Each simulation replicate begins with the population fixed for the AB haplotype and ends when the ab haplotype becomes fixed in the population. The mutation-selection model used for simulation is the twolocus, two-allele model of [18], shown in Figure 3A. This is the simplest possible model of compensatory evolution that allows reversible mutation.

In each generation the number of mutations occurring in the population in a generation is drawn as a Poisson(4Nµ) random variable, where µ is the per generation per individual per locus mutation rate. If more than one mutation is to be introduced, the haplotypes to be mutated are selected uniformly without replacement, meaning that two mutations cannot occur to the same copy. Mutations to a haplotype occur at the first or second locus with equal probability. The number of recombination events in each generation is drawn as a Poisson(2N ρ) random variable, where ρ is the recombination rate per generation per individual. For each recombination event, two haplotypes are chosen without replacement with probability proportional to their frequency, the two selected haplotypes recombined, and the two recombinant haplotypes returned to the population. Following mutation and recombination, one generation begets the next by resampling with replacement with probabilities proportional to the haplotype fitnesses (see Supporting Information for full details of population resampling).

After each round of mutation-recombination-resampling the composition of the population is evaluated to determine the allele frequencies. Whenever a haplotype becomes fixed in the population that is different from the last fixed haplotype, the generation and haplotype are recorded. If the last fixed haplotype before the simulation is completed was AB, it is a Type 2 event (Figure 2B). If the last fixed haplotype before completion was a deleterious intermediate (aB or Ab), it is a Type 1 event (Figure 2A).

Type 2 events encompass two scenarios: either the deleterious intermediate haplotype is lost from the population before the AB haplotype, leaving only the two high-fitness haplotypes, or the deleterious intermediate persists until after the AB haplotype is lost. While the latter case may result in one of the two new alleles becoming fixed before the other, as long as the intermediate haplotype does not fix it remains a Type 2 event (Figure 2C).

Simulations were performed under combinations of a range of mutation rates, recombination rates, and selection coefficients (Table 1). The simulation values chosen are similar to those of [19], and when scaled by the population size reflect reasonable estimates for a variety of natural populations. One thousand simulation replicates were performed for each combination of parameters (see Supporting Information for exceptions).

Table 1 Simulation Parameters

Results and Discussion

Direct double substitution probability

The direct double substitution probability β (Eq. 10) can be easily estimated from simulations. When fixed for the AB haplotype, the population can subsequently become fixed for either the ab haplotype or for one of the deleterious intermediate haplotypes. β can be estimated as the proportion of times in which the population becomes fixed for the ab haplotype after having been fixed for the AB haplotype.

Estimates of β under different simulation conditions are plotted as the points in Figure 4. For small mutation rates and low selection coefficients there is very little chance of observing a direct double substitution, while for high mutation rates and large selection coefficients there is little chance of observing anything but direct double substitutions.

Figure 4

The direct double substitution probability β. The analytical predictions using Equation 10 (lines) are consistent with simulations (points) for most mutation rates (colors), except when the mutation rate is very large.

The expected value of β can be derived from the rates of the continuous-time Markov model of Figure 3B. Because the rates out of state AB are given by Equations 1 and 8,

$$\beta =\frac{r_3}{2r_1+r_3}.$$

(10)

Analytical expectations for β from Equation 10 are shown in Figure 4 for comparison to the simulation results. The analytical result is a very good approximation to the observed direct double substitution probability across the range of selection coefficients, especially when the mutation rate is relatively low. The poor approximation of the analytical prediction when the mutation rate is very high is to be expected since the derivation of α assumes relatively low mutation rates.

Alternate compensatory substitution pathway probabilities

The probability that compensatory substitutions occur without the deleterious intermediate first becoming fixed (Type 2) is a function of selection, mutation, and recombination, all scaled by population size. Importantly, it is not the same as the direct double-substitution probability β. This is because the process may visit other states before proceeding from AB to ab directly. Consequently, β serves as a lower bound on the probability that a compensatory substitution event is of Type 2.

A more precise estimate for the probability of a compensatory substitution event being of Type 2 can be derived from β as well. Given that the population is in state AB, the probability that it moves directly to ab without returning to AB is β. The probability that it moves from AB to ab via the deleterious intermediate, without returning to AB is (1 - β)/ 2. Thus the relative probability that it moves from AB to ab directly before moving to ab via the deleterious intermediate is

$$\mathbb{P}\left(\mathsf{\text{Type}}2\right)=\frac{\beta }{\beta +\left(1-\beta \right)2}=\frac{2\beta }{1+\beta }.$$

(11)

The population could of course have returned to the initial state AB from the deleterious intermediate as well. But because the process is Markovian, the probability of observing a Type 2 event given that the current state is AB is still given by Equation 11.

The analytical prediction of Equation 11 fits the results of simulations in the absence of recombination quite well across the range of selection coefficients explored, except (unsurprisingly) when the mutation rate is very high. Like the direct double substitution probability β, the proportion of Type 2 compensatory substitution events increases with the strength of selection (Figure 5). In addition, Type 2 compensatory substitution events are very rare when the mutation rate is relatively low, such that the vast majority of compensatory substitution events involve the fixing of the deleterious intermediate haplotype.

Figure 5

Proportion of Type 2 compensatory substitutions. The analytical prediction of the probability of Type 2 events from Equation 11 (black lines) compares favorably to simulations without recombination (blue circles). Simulations with recombination are depicted as well (red triangles). (A) θ = 1 (B) θ = 0.1 (C) θ = 0.01 (D) θ = 0.001

Recombination is often thought of as having the positive effect of putting together beneficial mutations [27, 28]. Even negatively-epistatic alleles that are beneficial when paired will be aided by recombination [15]. But when the two mutants are individually deleterious and jointly neutral, recombination increases the time required for compensatory evolution [13]. Empirically, fewer compensatory substitutions are observed as physical distance (and hence recombination) between interacting loci increases [29, 30].

In this context, recombination reduces the proportion of Type 2 substitutions observed (Figure 5). Since the a and b alleles spend most of their time at low frequency due to selection against them, recombination serves to break apart ab haplotypes much more than it creates them. Consider a population with a single copy of each of the a and b alleles. A single ab haplotype will be broken up by recombination if it is one of the two selected, but an ab haplotype will only be created by recombination if both of the haplotypes bearing the deleterious alleles are selected. This tendency to separate rare alleles increases the chance that at least one of the alleles will be lost, decreasing the probability of Type 2 events. This effect of recombination is obviated when the mutation rate is very high, where Type 2 events dominate over a wide range of selective strength.

Taken together these paint a picture of the dynamics of the compensatory neutral substitution process. If mutation rates are relatively low, then times to compensatory substitution will be very large, particularly when selection is strong [17]. In this case the process will likely proceed through the deleterious intermediate fixed state, unless selection is strong (Figure 5). Recombination, if strong, will increase both the time to fixation [19] and the probability of passing through the deleterious intermediate state (Figure 5). Conversely, high mutation and strong selection will lead to faster Type 2 substitution events [13], and high mutation rate alone provides immunity from the effect of recombination (Figure 5). This suggests potentially very different compensatory substitution dynamics in different species.

Reversions to the ancestral state

It was noted earlier that the population, after fixing a deleterious intermediate state, could revert to the original fixed state AB rather than proceeding to the alternative state ab. In fact, it can cycle back a number of times before eventually becoming fixed for ab, either via the intermediate or directly. Each time in state AB the probability of "failure" (leaving then returning to AB) is (1 - β)/ 2, and the probability of "success" (reaching ab by any means) is (1 + β)/ 2. The number of returns to AB before ultimately reaching ab is therefore a geometrically distributed random variable on the space {0,1,2,...} with mean (1 - β)/(1 + β). Because only the final path from AB to ab determines if an event is of Type 1 or Type 2, the expected number of reversions is the same for both pathways.

Figure 6 shows the average number of times in simulation that the population fixes a deleterious intermediate state and then returns to state AB before fixing in state ab via any pathway. The number of returns to the initial state is lower when the mutation is high. This is because high mutation increases the rate of appearance (and subsequent fixation) of ab haplotypes (Figure 5); if intermediates do not become fixed, there can be no return to the initial condition. Recombination had little effect on the number of reversions, slightly increasing the number of reversions when selection is moderate to large and mutation is low (Additional file 2). Increasing selection against intermediates reduces the number of reversions by removing deleterious alleles from the population; by preventing deleterious mutants from fixing in the population it also prevents reversions (as well as Type 1 fixation events).

Figure 6

Reversions to the ancestral state. The average number of times returning to state AB after having been fixed for a deleterious intermediate state are shown under four mutation rates (panels) across a variety of selection coefficients. The analytical predictions (lines) from β of the number of recursions match well with results from simulations without recombination (points), except when the mutation rate is very large.

Relative time of the two pathways

An important question is whether the expected amount of time to compensatory substitution differs for Type 1 and Type 2 events. It is not obvious that these two events should take the same amount of time, or if differences between them should remain constant across a wide range of selection, mutation, and recombination. Because the expected number of reversions to the ancestral state is the same regardless of the pathway by which fixation is ultimately achieved, the power to detect differences between Type 1 and Type 2 events is maximized by comparing the time required for the final path from AB to ab.

The number of generations from the last time the population enters state AB until the time it enters state ab was calculated for each simulation replicate. If there is genuinely no difference in the expected time for Type 1 and Type 2 fixation events, then the difference in the two sample means is expected to be zero. Under certain assumptions about the samples (normally distributed and equal variances), the difference in means, normalized by the pooled variance, will be a standard normal random variable. Large deviations in this test statistic imply significant differences in the expected times of the two pathways.

The important observation is that natural selection against deleterious intermediates eliminates the difference in time required between the two paths observed, and this appears true across all mutation rates examined (Figure 7). With increased selection the probability of Type 2 events increases (Figure 5), but although Type 1 events are less likely to occur under strong selection, when they do occur they take the same amount of time as Type 2 events. This implies that the time to fixation alone is insufficient for determining the pathway by which fixation occurred if selection has been moderate. However, if selection has been very weak and the mutation rate is very high, it might be possible to use time to fixation to determine the fixation pathway.

Figure 7

The normalized time difference between Type 1 and Type 2 substitution events. The test statistic is the difference in mean times for the two paths, normalized by the pooled standard deviation, and is shown under four mutation rates (panels) across a range of selection coefficients. Shown are simulations in the absence (blue) and presence (pink) of recombination.

Recombination had little effect on either the differences observed under weak selection or the lack of differences observed at stronger selection (Figure 7). Because recombination most strongly affects the probability of observing Type 2 events when the mutation rate is low and selection is strong (Figure 5), if recombination is to have an effect on the difference in path times it should be most evident under these conditions. However, it is under precisely these conditions that the difference between path times is minimized (Figure 7).

Extensions

Model asymmetry

The total rate of compensatory substitution α in the more general case in which mutation for the two loci and selection coefficients against the two deleterious haplotypes could differ was previously derived [19]. This is attractive because it would allow for more biologically realistic scenarios, such as the semi-stable U-G intermediates in RNA stem pairing [23]. However, in deriving the more general case they made the assumption that Ns was very large in order to rule out the possibility of the coexistence of both deleterious haplotypes. This assumption might be problematic for many systems in which either population sizes or selection coefficients are not expected to be very large. Nevertheless, it should be straightforward to incorporate into the work presented here.

The case in which the fitness of ab is far greater than that of any other state, such as when normal cells mutate to cancerous ones, was considered by [22]. Under these conditions they derived a rate of double substitution, giving it the colorful term "stochastic tunneling." However, unlike the approach shown here, these authors assumed that the probability of fixation of a new double mutant was independent of the population composition at the time it arises. This assumption simplifies matters greatly, but such an approximation may be limited to the case in which the ab state is of exceptionally high fitness.

Expanding the state space

The two-locus, two-allele model considered here will be insufficient for some important biological cases, such as the evolution of RNA stem positions. In this case, either the state space of the model can be expanded to a two-locus, four-allele model, or the state space of RNA can be compressed somehow to fit the current model. The latter approach is simpler, but risks losing valuable information about the evolutionary process by grouping states together. The former approach for the case of strong selection and large populations was explored by [18].

A Phylogenetic substitution model

Because the rates derived in this work are rates of the substitution process, a phylogenetic model may be defined using these rates. In this case, the elements of the instantaneous rate matrix Q = {q _ij } of the continuous-time Markov chain can be defined as

$$q_{ij}=\left\{\begin{array}{cc}u{\pi }_j{r}_1\hfill & \mathsf{\text{if}}i\in \mathbb{Y}\mathsf{\text{and}}j\in \mathbb{X}\hfill \\ u{\pi }_j{r}_2\hfill & \mathsf{\text{if}}i\in \mathbb{X}\mathsf{\text{and}}j\in \mathbb{Y}\hfill \\ u{\pi }_j{r}_3\hfill & \mathsf{\text{if}}i=AB\mathsf{\text{and}}j=ab\hfill \\ u{\pi }_j{r}_4\hfill & \mathsf{\text{if}}i=ab\mathsf{\text{and}}j=AB\hfill \\ 0\hfill & \mathsf{\text{if}}i,j\in \mathbb{X}\mathsf{\text{and}}i\ne j\hfill \\ -\sum _{i\ne j}{q}_{ij}\hfill & \mathsf{\text{if}}i=j\hfill \end{array}\right.$$

(12)

or alternatively, in terms of the population genetic parameters as

$$q_{ij}=\{\begin{array}{ll}u{\pi }_{j}2N\mu \frac{1-e^{2\text{s}}}{1-e^{4N\text{s}}}\hfill & \text{if}i\in \mathbb{Y}\text{and}j\in \mathbb{X}\hfill \\ u{\pi }_{j}2N\mu \frac{1-e^{-2\text{s}}}{1-e^{-4N\text{s}}}\hfill & i\text{f}i\in \mathbb{X}\text{and}j\in \mathbb{Y}\hfill \\ u\frac{{\pi }_j}{{\pi }_i}2N\mu \left({\displaystyle {\int }_0^1\frac{x(1-e^{-2\text{s}(2N(1-x)+1))}}{(1-e^{-4N\text{s}})(2N(1-x)+1)}}\frac{e^{4N\text{s}x}{(x(1-x))}^{8N\mu -1}}{{\displaystyle {\int }_0^1{e}^{4N\text{s}z}}{(z(1-z))}^{8N\mu -1}dz}dx-\frac{(1-e^{-2\text{s}}){\pi }_{\mathbb{X}}}{1-e^{-4N\text{s}}}\right)\hfill & i\text{f}i,j\in \mathbb{Y}\text{and}i\ne j\hfill \\ 0\hfill & i\text{f}i,j\in \mathbb{X}\text{and}i\ne j\hfill \\ -{\displaystyle \sum _{i\ne j}{q}_{ij}}\hfill & i\text{f}i=j\hfill \end{array}$$

(13)

where $\mathbb{X}=\left\{aB,Ab\right\},\mathbb{Y}=\left\{AB,ab\right\}$, and u is a rate scaling factor to ensure that branch lengths are in terms of expected number of substitutions per site.

Phylogenetic models of the compensatory substitution process have been developed previously [31, 32], but these models have double and single substitution rates that are uncoupled from each other and that are purely descriptive in nature. Equation 13 on the other hand represents a fully-defined phylogenetic substitution model for the analysis of compensatory substitution specified entirely in terms of the underlying population genetic parameters. This is appealing because it moves towards models that are interpretive, in which parameters have meaning with respect to the process of evolution, rather than models that are simply descriptive of patterns observed.

The model described by Equation 13 presents a potentially powerful way to investigate the evolution of compensatory evolution using interspecific data, allowing the estimation of population parameters of interest. Application of the basic model described here to interspecific data would be making key assumptions about the evolutionary process, namely that each branch of the phylogenetic tree would share a common population size, selection coefficient, and mutation rate. Relaxing this constraint to allow different branches to have their own parameters has been done in different contexts using the Dirichlet process [33, 34] and could readily be applied to this model as well.

Origination vs. fixation process

The total compensatory substitution rate α is actually the rate at which successful compensatory haplotypes arise in the population, the inverse of which being the time until such haplotypes arise [19]. This is slightly different from what was calculated in simulations, which was the time until the actual fixation event itself. The difference in these times is very small: A compensatory (second) mutation that is destined to go to fixation, from the time it arises in the population, takes on the order of 4N generations until it becomes fixed [19, 35], which is trivial compared to the time typically spent waiting for the eventually-successful mutation to arise. Consequently the simulations should be expected to yield ever-so-slightly longer fixation times than predicted by α.

These two processes are distinguished as origination and fixation processes, respectively [36]. Origination processes can often have more desirable statistical properties than fixation processes. For example, in the case of neutral linkage the origination process remains Poisson while the fixation process is something more complex and overdispersed [37]. While substitution models are rarely explicit about the process being modeled, the theoretical results and phylogenetic substitution model presented here pertain to the origination process rather than the fixation process.

Conclusions

There are many cases in which loci interact epistatically and natural selection is weak to moderate, and this work derived rates for the compensatory neutral substitution process along the different possible paths to fixation under these conditions. These were shown to have good predictive power regarding the probability of observing compensatory substitution via alternate paths under a range of conditions. Because there is little difference in the time to fixation between the paths when selection is not negligible, the relative probability of these paths may be important for distinguishing the evolutionary history of compensatory substitutions. To this end, a phylogenetic substitution model parameterized in terms of these population genetic parameters was presented that could be used for inference of population histories using interspecific data.