Copyright Ó 2011 by the Genetics Society of America
DOI: 10.1534/genetics.110.125997

Epistasis Increases the Rate of Conditionally Neutral Substitution
in an Adapting Population
Jeremy A. Draghi,*,1 Todd L. Parsons* and Joshua B. Plotkin**Department of Biology, University of Pennsylvania, Philadelphia, Pennsylvania 19104
Manuscript received December 16, 2010
Accepted for publication January 24, 2011
ABSTRACT
Kimura observed that the rate of neutral substitution should equal the neutral mutation rate. This classic
result is central to our understanding of molecular evolution, and it continues to influence phylogenetics,
genomics, and the interpretation of evolution experiments. By demonstrating that neutral mutations
substitute at a rate independent of population size and selection at linked sites, Kimura provided an
influential justification for the idea of a molecular clock and emphasized the importance of genetic drift inshaping molecular evolution. But when epistasis among sites is common, as numerous empirical studies
suggest, do neutral mutations substitute according to Kimura’s expectation? Here we study simulated,asexual populations of RNA molecules, and we observe that conditionally neutral mutations—i.e.,mutations that do not alter the fitness of the individual in which they arise, but that may alter the fitness
effects of subsequent mutations—substitute much more often than expected while a population isadapting. We quantify these effects using a simple population-genetic model that elucidates how thesubstitution rate at conditionally neutral sites depends on the population size, mutation rate, strength of
selection, and prevalence of epistasis. We discuss the implications of these results for our understanding of
the molecular clock, and for the interpretation of molecular variation in laboratory and natural
populations.

K

IMURA’S observation that the rate of substitution at
a neutral site should equal the neutral mutation
rate is one of the most elegant and widely applied resultsin population genetics (Kimura 1968; Kimura and Ota
1971; Bromham and Penny 2003; Hughes 2008; Neiet al. 2010). This theory performs well for sites in agenome that can be classified as unconditionally neutral:that is, sites at which the fitness effects of mutations are
negligible in any environment, and in combination withany genetic background. But what does neutral theory
predict about the fate of a mutation that is known to be
neutral only in the genetic background in which it arose?
Such mutations may interact epistatically with subsequent mutations at other loci and are thus called
conditionally neutral. In light of recent studies supporting a constructive role for such epistatic neutral variationin adaptive evolution (Schuster and Fontana 1999;
Depristo et al. 2005; Koelle et al. 2006; Amitai et al.
2007; Cowperthwaite and Meyers 2007; Wagner
2008a; Bloom and Arnold 2009; Draghi et al. 2010), we
ask whether Kimura’s foundational result extends toconditionally neutral mutations.

Supporting information is available online at http://www.genetics.org/
cgi/content/full/genetics.110.125997/DC1.
1
Corresponding author: Department of Biology, University of Pennsylvania, 204K Lynch Labs, 433 S. University Ave., Philadelphia, PA 19104.
E-mail: jdraghi@gmail.com
Genetics 187: 1139–1152 (April 2011)

To understand the generality of Kimura’s result, it is
helpful to consider an informal derivation. Imagine an
idealized population of N haploid individuals, one ofwhich will eventually be the ancestor of the future
population. If unconditionally neutral mutations occur
at rate m per replication, then on average Nm mutations
will arise in the population each generation. Because
these mutations can never affect fitness, they cannot
affect the eventual fate of the lineages in which theyarise. Therefore, each unconditionally neutral mutation
will arise in the eventual common ancestor with probability 1/N; otherwise, it will be lost. The average rate ofneutral substitution, k, therefore, equals the rate of(unconditionally) neutral mutation times the fixation
probability of each mutant:1
k ¼ N m ¼ m:
N

ð1ÞThe reasoning behind Equation 1 is compelling, andmany studies have argued that this result holds for sexualand asexual species, for neutral mutations linked topositively or negatively selected sites, and for populations of varying sizes (Kimura and Ota 1971; Birky andWalsh 1988; Gillespie 2000; Bromham and Penny
2003). As a result, the rate of substitution at neutral sites
is now viewed as one of the most robust and wellunderstood features of molecular evolution. Extensions

1140

J. A. Draghi, T. L. Parsons and J. B. Plotkin

to the neutral theory have mainly focused on the apparent
overdispersion of neutral substitutions (Gillespie 1986,
1993; Takahata 1987; Bastolla et al. 1999, 2002,
2003; Cutler 2000; Wilke 2004; Bloom et al. 2007;
Raval 2007). With the exception of a few studies thatpredict small deviations in models with lethal mutations
and stabilizing selection (Bastolla et al. 1999; Bloom
et al. 2007), most work has confirmed or, more often,
tacitly assumed that Equation 1 accurately describes
the mean substitution rate. These studies have largely
ignored the impact of conditionally neutral mutations:mutations that are neutral on the genetic backgroundin which they arise, but that may alter the fitness effectsof subsequent mutations. If neutral mutations haveepistatic interactions of this sort, then it is unclear
whether Kimura’s equation describes their substitution
rate.
A diverse array of recent computational and empiricalstudies has demonstrated the importance of neutral
mutations with epistatic effects (reviewed in Wagner
2008a). Evolutionary simulations with RNA folding
algorithms (Huynen 1996; Huynen et al. 1996; Fontana
and Schuster 1998; Ancel and Fontana 2000; Wagner
2008b) and model gene networks (Bergman and Siegal
2003; Ciliberti et al. 2007) indicate that neutral changes
may often be prerequisites for adaptive substitutions
and that the interactions between neutral and adaptivechanges can lead to complex dynamics of phenotypicevolution; theoretical developments have generalized and expanded these results (van Nimwegen
and Crutchfield 2000; Lenski et al. 2006; Wagner
2008a,b; Weissman et al. 2009; Draghi et al. 2010).
Additional evidence comes from laboratory evolution
experiments with proteins, in which apparently neutral
mutations permit future adaptations by changing
thermodynamic stability, codon usage, or promiscuous
protein–ligand interactions (Depristo et al. 2005;
Bloom et al. 2006; Amitai et al. 2007; Cambray andMazel 2008; Bloom and Arnold 2009). The epistastic
effects of nearly neutral mutations can even explain theevolution of consequential innovations, such as adaptive
expansion into a new niche (Blount et al. 2008), thesudden escape of a pathogen from population immunity(Koelle et al. 2006; van Nimwegen 2006; Kryazhimskiy
et al. 2011) or susceptibility to a drug (Bloom et al. 2010;
Kryazhimskiy et al. 2011).
If some neutral mutations can facilitate future adaptation through epistatic interactions, selection might drivethese neutral mutations to fixation by hitchhiking—that is, by linkage to subsequent beneficial mutations.
However, other neutral mutations will impede future
adaptive changes, and fixation of these neutral mutations would be disfavored by selection. In each case, theeffects of a mutation on an individual’s evolvability—that is, its capacity for adaptation—causes its probabilityof fixation to be larger or smaller than that of an
unconditionally neutral mutation. Naively, one might

expect that conditionally neutral mutations would be no
more likely to enhance evolvability than to diminish it.
Consequently, the effects of evolvability on the fixationof these mutations might average out, and Equation 1
might accurately describe the substitution rate ofepistatic neutral mutations. Here we show that this naive
expectation is incorrect. Instead, ‘‘neutral epistasis’’ inan asexual, adapting population causes a significant
elevation of the substitution rate at conditionally neutral
sites, compared to Kimura’s classical expectation for
unconditionally neutral sites. We first demonstrate thesedepartures from the conventional substitution rate insimulated populations of replicating RNA molecules,
and we confirm that the substitution rate is caused by the
epistatic effects of neutral mutations. We then explore asimple population-genetic model that quantifies how
epistasis, population size and mutation rate, and selection coefficients jointly determine the substitution rate
at conditionally neutral sites in adapting populations.
Finally, we discuss the implications of these results forthe molecular clock and for the inference of evolutionary processes in natural and laboratory populations ofnonrecombining organisms and chromosomes.

METHODS

Measuring substitution rates: Following Gillespie andothers (Gillespie 1993), we distinguish between two
types of substitution events. An origination event is thefirst appearance of a genotype which will later beancestral to everyone in the population. Inspecting a
genealogy, we say that such mutant genotypes were‘‘destined’’ to fix; when they do fix, we mark a fixation
event. Each fixation event corresponds to an earlier
origination event. Because of this correspondence, themean rate of origination events and the mean rate of
fixation events will eventually converge, and such convergence defines a steady-state population with respect tosubstitution. We illustrate transient and steady-state substitution dynamics in supporting information, Figure S1.
Linkage and changes in mutation rate may cause the rateof one process to temporarily exceed the other, in whichcase we may choose to inspect whichever process is more
informative.
We quantified substitutions using origination events.
Specifically, we choose an individual at random from thefinal generation and trace its ancestry back to the initial
generation. This lineage eventually passes through themost recent common ancestor (MRCA) of the final
population. Each mutation encountered before the
MRCA represents an origination event. If neutral mutations were unconditionally neutral, we would expect
them to arise on this lineage at their mutation rate. We
therefore measured the neutral mutation rate of eachparent along the lineage as our null expectation of theneutral substitution rate. Using the origination process

Episistasis and Substitution Rates

instead of the fixation process allows us to calculate a null
expectation that is robust to variation in neutral mutation rates and to selective bottlenecks.
In our simulations we did not directly determine which
genotype is the MRCA of the final population. For the
RNA results, we tracked a lineage back from a single
individual: the population size, N, and the rate of beneficial substitution determine how quickly this lineage
converges to the MRCA. As a consequence, we may havemissed some potential substitutions that arose late ineach simulation and erroneously classified some polymorphic variants as neutral substitutions. However, such
events could only cause the measured neutral substitution rate to tend toward the null expectation (i.e., theneutral mutation rate). We therefore expect that our
RNA results on the neutral substitution rate may besomewhat conservative when t is close to the end of
the simulation (30,000 generations). In the simple
population-genetic model we tracked substitutions inreal time, and we ran each simulation until the firstorigination event after the nominal ending time, T. We
could then be assured that all mutations arising before T
had been fixed or lost and that the origination events
measured in these T generations were an unbiased
estimate of the substitution rate. In both cases, themeasured substitution rates correspond to those expected with infinite sites; substitutions that occur at thesame base in the RNA model would each be counted.
Simulation methods: All simulations used the
Wright–Fisher model to evolve populations of discrete,
asexual, and haploid individuals. RNA simulations used
the Vienna RNA folding package (version 1.6.1) withdefault folding parameters. These simulations considered RNA sequences of 72 bases in length and computed
the minimum-free-energy structure to determine thephenotype. Fitness was then calculated as a function of
the tree edit distance, d, between an organism’s phenotype and a defined optimal phenotype. The tree edit
distance algorithm, included in the Vienna package,determines the minimum number of steps from a groupof edit operations that are needed to transform one
structure into another. Fitness is related to structural
distance as (1 1 s)d. Here s quantifies the strength of
selection; it is equivalent to the multiplicative selective
coefficient associated with a mutation that changes d bya single unit. The initial genotype was drawn randomly,
and its phenotype was defined as the optimum for the
initial period of stabilizing selection. The second optimum phenotype, used to impose directional selection,
was also created by randomly drawing genotypes anddiscarding those whose minimum-free-energy structureis the trivial, unfolded state. The second optimum wasalso required to be 40 units from the first optimum, sothat the pressure to adapt in the new environment wasstrong and uniform across replicates. Populations were
initially clonal, and each replicate began from an
independently drawn genotype and later adapted to1141

an independently drawn optimum. Substitutions were
measured as described above.
We also performed simulations of the simple
population-genetic model. Selection and reproduction are modeled identically in the two kinds of simulations. Again, the parameter s quantifies the strength of
selection by specifying the multiplicative fitness advantage associated with each beneficial mutation.
In both the abstract and RNA models, mutations
occur as a Bernoulli process: either no mutations occur,or a single mutation arises with probability U. Although
this simplification is biologically unrealistic, it enables us
to unequivocally assign a fitness effect to each mutation.
Because a neutral mutation cannot arise at the sametime as a beneficial mutation, unconditionally neutral
mutations may actually substitute at a slightly lower rate
than their mutation rate. This downward deviation from
Kimura’s prediction was found to be negligible.
In a separate series of RNA simulations, we measured
the evolvability and robustness of the population as awhole, as well as the effects of neutral mutations onthese properties. The most straightforward method,which is to measure these properties in every individualand every neutral mutant, is seriously distorted by thepresence of deleterious variation in the evolving populations: individuals harboring deleterious mutations
may appear to be evolvable, but most often can produce
only mutants as fit as existing genotypes. We therefore
chose to measure robustness and evolvability only inthose individuals whose fitness exceeded the mean of
their population; we also assayed the effects of neutral
mutations only when those mutations arose from suchindividuals. In practice, due to the large disadvantage(at least 10%) of deleterious variants, the distribution of
fitnesses in a population is left skewed, and this heuristic
typically excludes much less than half of the individuals
in a typical population.

RESULTS

Neutral substitutions in simulated RNA populations:
Kimura’s classic result and its subsequent elaborations(e.g., Charlesworth et al. 1993; Gillespie 2000) apply
to unconditionally neutral sites. Although unconditionally neutral sites are a convenient theoretical concept,
empirical measurements of fitness effects could never
suffice to identify such sites in practice (Wagner 2005).
Even in the simplest genome, there are far too many
possible genetic backgrounds to feasibly test the fitness
consequences of a mutation at a given site on all backgrounds, as well as in a range of possible environments.
In practice, we must use an empirically tractabledefinition of neutrality. We classify a mutation as neutral
if it has no effect on fitness in the genetic backgroundand environment in which it arises. Two considerations
motivated us to study substitutions at sites satisfying this1142

J. A. Draghi, T. L. Parsons and J. B. Plotkinbroad definition of neutrality. First, the initial fitness
effect of a mutation strongly determines its eventual
fate, especially in large populations with little recombination. Second, this broad definition can be easily
applied in evolution experiments. The use of wholegenome sequencing in experimental evolution has
made neutral theory testable in new ways, and several
studies have attempted to combine experimental fitness
measurements with ideas from neutral theory (e.g.,
Wichman et al. 2005; Barrick et al. 2009). By considering a broad definition of neutrality, which includes
conditionally neutral sites, we examine the dynamics ofthe types of mutations that an experimentalist would, inpractice, typically classify as neutral.
To study substitutions of such (conditionally) neutral
mutations, we simulated an asexual population of 72-base
RNA molecules evolving according to the Wright–Fisher
model (see methods). We used the Vienna RNA software
package to assign a structure to each sequence and tomeasure differences between structures. We assigned
fitness to each sequence according to its corresponding
structure’s similarity to a chosen optimal structure (seemethods). Populations were subject to stabilizing selection for T1 generations; the optimal phenotype was thensubstantially altered, and the populations adapted to thisnovel environment for T2 generations. To precisely
compare our substitutions to Kimura’s expectation, we
measured substitutions by recording when they initiallyarise as mutations (see methods).
As Figure 1 shows, the rate of (conditionally) neutral
substitutions in an evolving population of RNA molecules departs significantly from Kimura’s expectation. In
particular, the substitution rate spikes sharply when theenvironment (i.e., the target phenotype) shifts, coincident with a period of rapid adaptation. For a populationof size N ¼ 5000, the neutral substitution rate is about
twice as large as the neutral mutation rate shortly after
the environmental shift, and it remains elevated forthousands of subsequent generations. Furthermore,
populations of different sizes exhibit distinct rates ofneutral substitution, again in contrast to Kimura’s expectation for nonepistatic neutral mutations.
The first clue to the mechanism causing this high rateof neutral substitution lies in the dynamics of adaptationin simulated RNA populations. While the ensemble mean
fitness graphed in Figure 1 increases steadily, adaptationin the individual populations is much less predictable. A
population of replicating RNA molecules in this regime
exhibits punctuated evolution: periods of constant fitness
are disrupted by the sudden appearance and fixation ofa beneficial substitution (Fontana and Schuster 1998;
Ancel and Fontana 2000). In our study, the ensemble
mean population fitness rapidly increased immediately
following the environmental shift (see Figure 1), because
many beneficial mutations were available at the time ofthe shift. Moreover, the ensemble mean population
fitness continued to grow, at a slower pace, throughout

Figure 1.—The observed rate of substitution at (conditionally) neutral sites depends on the population size and therate of adaptation in simulated RNA populations. Eachpoint reflects neutral substitutions originating within a 100generation bin, averaged over at least 20,000 replicates. Points
are calculated as the observed number of substitutions minus
the expected number, divided by the expected number; this isan estimate of origination rate relative to Kimura’s expectation for unconditionally neutral mutations. Solid circles correspond to N ¼ 5000, while open circles correspond to N ¼
800; for both, the genomic mutation rate, U, is equal to0.004 and s ¼ 0.1. The lines depict ensemble means of mean
fitnesses within each replicate; the solid line corresponds to
N ¼ 5000; and the dashed line corresponds to N ¼ 800.

the subsequent 15,000 generations—reflecting the rare
discoveries of beneficial mutations within the ensembleof replicates. Since beneficial mutations were strongly
selected and quite rare in this period, we expect thatthose neutral mutations with positive epistatic consequences (i.e., those neutral mutations that increased agenotype’s chance of receiving a subsequent beneficial
mutation) were effectively selected and fixed throughhitchhiking (i.e., by linkage to a beneficial mutation).
In light of the considerations above, it is important tonote that the epistatic consequences of neutral mutations
were, indeed, mostly positive in our simulations. In
particular, we observed that neutral mutations tended
to increase evolvability. We define the evolvability of agenotype as the fraction of its point mutations that
improve fitness. Although we observed many neutral
mutations that decreased evolvability, the average effect
was positive in all generations after the environmental
shift (Figure 2a). Analysis of the median effects confirmed that neutral mutations were more likely to increase evolvability than to decrease it. As many studies
have shown, mutations that increase evolvability arefavored by selection, especially in nonrecombining
populations adapting under strong selective pressures(Kashtan and Alon 2005; Meyers et al. 2005; Draghi
and Wagner 2008, 2009). These results suggest thatneutral mutations substitute more often than Kimura’s
expectation because of their biased effects on evolvability.

Episistasis and Substitution Rates

1143

Figure 2.—Neutral mutations have a mean
positive effect on evolvability in evolving RNA
populations. (A) The mean change in evolvability caused by a neutral mutation. Each point represents mutations observed in 3000 replicate
simulations, averaged into 500-generation bins.
(B) The fraction of replicate populations without
any available beneficial point mutations. Points
represent averages from observations at 500generation intervals, averaged across over 2000
replicates. N ¼ 5000, U ¼ 0.004, and s ¼ 0.1, asabove.

To understand why neutral mutations typically increased evolvability in our simulations, we quantified
how often the populations became ‘‘stuck’’ on fitness
plateaus, unable to further adapt. Figure 2b shows thefraction of replicate populations in which no individuals
can produce any beneficial point mutation. This fraction
is initially small at the time of the environmental shift but
grows quickly to over 80%. This behavior is not caused bya large fraction of the replicate populations reaching a‘‘local peak’’ on their fitness landscape and thereafter
remaining stuck at that fitness indefinitely. Instead, we
observed that most populations became temporarily
stuck on fitness plateaus and then subsequently evolve
the potential to find further beneficial mutations. In
particular, we measured evolvability along lineages thatsurvive until the end of the simulation. Such lineages
passed through a genotype with zero evolvability 4.11
times, on average, during a 15,000 generation period ofadaptation. Neutral mutations, as opposed to deleterious
ones, were primarily responsible for potentiating adaptation in populations that had become stuck (in 20
replicates of the above simulations, 55 of the 57 observed
substitutions that increased evolvability in unevolvable
populations were neutral, as opposed to deleterious).
This result supports the decision to measure evolvability
only in genotypes with fitness exceeding the population
mean (see methods). Thus, the most important consequence of neutral mutations was to facilitate adaptationin populations that had become temporarily stuck—aphenomenon that we explore in greater generality below.
A simple population-genetic model: In the RNA
simulations described above, populations often reachedunevolvable states, which were subsequently alleviatedby epistatic neutral mutations. To explore this process
more generally, we have devised a population-genetic
model that links neutral mutations, beneficial mutations, and evolvability. In addition to the RNA simulations, our modeling approach was influenced by ourprevious work on robustness and evolvability (Draghi
et al. 2010) and by Gillespie’s pseudohitchhiking model(Gillespie 2000). We model a population of N haploid
individuals evolving according to the Wright–Fisher
model. We consider all sites in which mutations could

be neutral, conditionally neutral, or beneficial, anddefine U as the total mutation rate of these sites. We
assume two categories of genotypes. Some genotypes are
evolvable: a proportion pB . 0 of their mutations are
beneficial, conferring a multiplicative fitness advantage
1 1 s. The other class of genotypes are unevolvable—i.e.,
no beneficial mutations are (directly) available to them.
When a beneficial mutation occurs in an evolvablegenotype, the mutant genotype becomes unevolvablewith probability a. Finally, a proportion pNE of mutations
in both categories of genotypes are neutrally epistatic,
and these mutations switch an individual from one
category to another (i.e., from evolvable to unevolvable,or conversely) without changing its fitness.
Under this simple model, the epistatic effects ofmutations are described by two parameters: the parameter a quantifies the chance that a beneficial mutation
arrives at a fitness plateau, with no immediate potential for further adaptation; whereas pNE quantifies the
chance that a neutral mutation is epistatic and therefore
alters an individual’s evolvability. Crucially, the epistatic
neutral mutations are not constrained to increase evolvability. In fact, if they arise in evolvable and unevolvable
backgrounds equally, neutral mutations will have nomean effect on evolvability. For this section we focus on
the substitution rate at the conditionally neutral sites—
i.e., the proportion pNE of mutations that are neutral butthat alter an individual’s evolvability. The remaining
fraction 1  pNE  pB of uncategorized mutations areunconditionally neutral and substitute at a rate equal
to their mutation rate.
Monte Carlo simulations of this simple model confirm that neutral epistasis is a sufficient mechanism withwhich to explain an elevated rate of neutral substitution.
Figure 3 shows the neutral substitution rate measured inour simplified model, in simulations using mutation
parameters similar to the RNA model. When neutral
mutations have no epistatic effects (i.e., when we set
pNE ¼ 0), the neutral substitution rate is indistinguishable from the neutral mutation rate, as expected (seemethods). However, when genotypes can be evolvable
or unevolvable, and some neutral mutations switch
between these categories (i.e., pNE . 0), then the neutral

1144

J. A. Draghi, T. L. Parsons and J. B. Plotkin

substitutions, denoted t. To do so, however, we maketwo simplifying approximations. We assume that the
fixation of beneficial mutations occurs immediately andthat all epistatic neutral substitutions occur by hitchhiking, and not solely by neutral drift. (We expect theseapproximations to be accurate provided s is large andclonal interference is negligible.) These assumptions
lead to a recurrence: E½t ¼ aE½T un  1 ð1  aÞðE½T ev  1
E½tÞ, and thus the expected substitution rate at conditionally neutral sites is
k cond 

Figure 3.—The substitution rate of epistatic neutral mutations, relative to their mutation rate, is elevated under directional selection in our simple population-genetic model. The
solid circles show the mean origination rate measured in simulations that allow neutral mutations to change an individual’s
evolvability (pNE ¼ 0.01). The open circles show control simulations in which neutral mutations have no effect on evolvability (i.e., pNE ¼ 0; all genotypes are evolvable). In analogy
to our RNA experiments, populations evolved under stabilizing selection for 15,000 generations (i.e., pB was set to zerofor this period), and then under directional selection foranother 15,000 generations. N ¼ 5000, U ¼ 0.004, s ¼ 0.1,
and pB ¼ 0.001.

substitution rate is significantly higher than Kimura’s
expectation.
We have derived an analytic expression that predicts
the substitution rate at conditionally neutral sites in thissimple population-genetic model (See appendix:
derivation of equation 2). We first derived an expression for the expected waiting time until the first
beneficial mutation that is destined to fix arises, startingfrom a population of N unevolvable individuals:
E½Tun  ¼

pﬃﬃﬃﬃﬃ 22uNE 1 GðuNE Þ2:
N pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2uB pðsÞ Gð2uNE Þð2ÞHere uNE ¼ NUpNE, uB ¼ NUpB, and pðs Þ ¼ ð1  e 2s Þ=
ð1  e 2Ns Þ. Because an unevolvable population requires neutral changes in evolvability before it canadapt, we expect that one conditionally neutral mutation will hitchhike with the first beneficial substitution
in an initially unevolvable population.
If, on the other hand, a population is composed ofentirely evolvable genotypes, then the waiting time
before the next beneficial substitution is geometric withmean: E½T ev  ¼ 1=uB pðs Þ. Because all individuals in sucha population are already adaptable, we expect that noepistatic neutral mutations will fix during this interval;
such mutations would make their genotypes unevolvable and therefore could not fix through hitchhiking.
We can use our expressions for Tun and Tev, above, tocalculate the expected time between epistatic neutral1a¼:
E½t aE½T un  1 ð1  aÞE½T ev 

ð3ÞTo quantify the magnitude of the changes in theneutral substitution rate caused by epistasis, we simulated this simple model for several combinations ofparameters. Figure 4 shows the results of these simulations along with our analytical predictions (Equation3). In the RNA simulations, we were limited by computation time to modeling strong selection (i.e., s ¼ 0.1)
and evolution toward a distant optimum. Here, using
our simple population-genetic model we generalize theRNA results to explore weaker selection and lower rates
of beneficial mutation—that is, parameters that correspond to less dramatic environmental changes. We findthat even modest selective coefficients can drive substantial increases in the neutral substitution rate: linkage
to selected mutations conferring a 1% benefit can
produce a twofold increase in the substitution rate ofepistatic neutral mutations. Larger selective coefficients
can increase the substitution rate of neutral mutations
by 10-fold or more. Moreover, increasing the beneficial
mutation rate of evolvable genotypes greatly enhances
the elevation of the neutral substitution rate. Thissuggests that the magnitude of the change in evolvability
caused by neutral mutations strongly affects their substitution rate. Finally, we find that increasing the rate ofepistatically neutral mutations diminishes their percapita fixation probability, although this change is small
when uNE > 1 (Figure 4).
Figure 4 also plots predictions derived from theanalytical approximation above. In general, the predictions closely match the simulations, despite theapproximations we have used. Note that our simple
population-genetic model allows for clonal interference
that is, competition between beneficial mutations indifferent lineages in an asexual population—while ouranalytical approximation assumes that only one beneficial mutation is segregating at any time. Clonal interference may explain the departures of the prediction
from the simulated data, particularly when s and uB arelarge. This close agreement suggests that differences in
evolvability, combined with hitchhiking on selected
mutations, largely determines the dynamics of epistatic
neutral mutations over a broad range of parameters.

Episistasis and Substitution Rates

1145

Figure 4.—Simulated and predicted rates of substitution of epistatic neutral mutations. Open circles show mean origination rates across 4000 replicate simulations of 100,000 generations each. The dashed lines and X’s show analytical predictions calculated from equation 3. N ¼ 1000 and U ¼ 0.01 for all simulations; s ¼ 0.01 is the smallest selection coefficientsimulated.

Although we have ignored deleterious mutations inthis simple model, they could interact with beneficial
and neutral changes to influence our results in several
ways. We explored three scenarios for the parameters
used in Figure 3. First, we added a large fraction ofdeleterious mutations (pD ¼ 0.5) to the mutation
probabilities for both evolvable and unevolvable backgrounds. Therefore, in an evolvable background, afraction pD mutations are neutral, pNE are conditionallyneutral, pB are beneficial, and (1  pD  pNE  pB) areunconditionally neutral; in an unevolvable genotype, pD
mutations are neutral, pNE are conditionally neutral,none are beneficial, and (1  pD  pNE) are unconditionally neutral.
We found that adding the possibility of deleterious
mutations had little effect, regardless of whether thesemutations have large (s ¼ 0.2) or small (s ¼ 0.02) fitness
effects (see Figure S2). We also explored the possibility
that evolvable genotypes might be less robust to mutation, and so have a higher probability of deleterious
mutation. In this scenario, the rate of conditionallyneutral substitution is lower than in the absence ofdeleterious mutations, but still much higher than Kimura’s expectation (see Figure S3).
Finally, we modeled epistatic deleterious mutations.
These mutations had a direct fitness cost, but likeconditionally neutral mutation, would change unevolvable backgrounds into evolvable genotypes and viceversa. When such mutations are at least modestly
deleterious, the rate of conditionally neutral substitution
is still substantially elevated compared to Kimura’s
expectation (see Figure S4).
Neutral epistasis and stabilizing selection: In thissection we discuss another distinct anomaly in neutral
substitutions observed in our RNA simulations: the increase in substitutions prior to the environmental shift(Figure 1). This effect is quite striking in the large
population we simulated: an increased origination rateof neutral mutations destined to fix is apparent for as
many as 1000 generations prior to the environmental
shift. This pattern may in part be caused by mutations
that segregate neutrally in the first environment, butthat become beneficial in the later selective environment and are then directly selected and brought to1146

J. A. Draghi, T. L. Parsons and J. B. Plotkinbetween robustness and evolvability in a large population under stabilizing selection (Figure 6).

DISCUSSION

Figure 5.—Neutral mutations that arise in a population under stabilizing selection have a mean positive effect on evolvability in some future environment. Evolvability, on the y-axis,
is measured relative to the target phenotype in the future environment. Mutant effects are averaged in 500-generationbins from the 15,000 generation preceding the environmental shift and averaged across 3000 replicate simulations. N ¼
5000, U ¼ 0.004, and s ¼ 0.1, as above.

fixation. Another possibility is that neutral mutations
segregating prior to the environmental shift may increase evolvability after the shift, and thus sweep tofixation by hitchhiking on subsequent beneficial mutations. Figure 5 shows that this latter phenomenon does
indeed occur: the neutral mutations that arise before
the environmental shift tend to increase an individual’s
evolvability after the shift. But why should neutral
mutations arising in population under stabilizing selection be biased toward increasing an individual’s evolvability in some future, novel environment?
The key to understanding this unusual phenomenon
is that populations under stabilizing selection evolve to
reduce the negative effects of deleterious mutations. As
predicted by theoretical models of neutral networks(van Nimwegen et al. 1999) and observed in otherstudies (e.g., Plotkin et al. 2006; Draghi and Wagner2009), our large RNA populations evolved high amountsof robustness while under stabilizing selection (Figure6). Theory also predicts that when every genotype canproduce beneficial mutations, robustness and evolvability will be negatively correlated (Draghi et al. 2010);
indeed we find in our simulations that, under stabilizing
selection, the effects of a neutral mutation on robustness
in the current environment are negatively correlatedwith evolvability in the future, novel environment (R 2 ¼
0.993, Spearman rank correlation). When robustness
in the population is high, the typical neutral mutation
will tend to reduce robustness (in the current environment), which has the the consequence of increasing
evolvability (in the future environment). Thus, theelevated neutral substitution rate observed before theenvironmental shift can be explained by the trade-off

We have demonstrated that Kimura’s classic expectation for the substitution rate of strictly neutral mutations
substantially underestimates the substitution rate ofconditionally neutral mutations that are linked to
selected sites, when populations are adapting in amutation-limited regime. Although we initially used a
computational model of RNA to illustrate this phenomenon, our results are not specific to that genotype–phenotype model, or to the relatively strong selectionused in those simulations. In general, we have shownthat whenever adaptation in a population is limited by
the supply of beneficial mutations, genotypes will
typically have lower evolvability than mutationally adjacent genotypes with the same fitness, and so conditionally neutral mutations will tend to increase evolvability
and fix more often than their mutation rate. Thisbehavior will occur in any genotype–phenotype-fitness
landscape in which selectively equivalent genotypes arelinked by mutation and differ in the phenotypes foundin their mutational neighborhoods. There is ample
evidence that such landscapes, in which neutral ornearly neutral mutations can potentiate adaptation,
are common at the level of proteins (Depristo et al.
2005; Wagner 2005, 2008a; Bloom et al. 2006, 2010;
Koelle et al. 2006; van Nimwegen 2006; Amitai et al.
2007; Cambray and Mazel 2008; Bloom and Arnold
2009; Kryazhimskiy et al. 2011) and at the level of gene
networks and metabolism (Bergman and Siegal 2003;
Ciliberti et al. 2007; Blount et al. 2008). The prevalence of epistasis linking neutral and adaptive changes,
and the consequences we have documented for the
resulting patterns of molecular evolution, should motivate further effort to measure the epistatic consequences of neutral mutations.
There is a long history of considering the waiting timeto the arrival of a double mutant. The first approach(Gillespie 1984; Kimura 1985; Weinreich and Chao2005) assumes that the first mutation is deleterious andthe second is compensatory, under the so-called strong
mutation–weak mutation regime: s ? 1/N ? m. By neglecting back mutation, and assuming that a sufficiently large
population size and a sufficiently long time has elapsed
since the first appearance of the mutation, they argue
that the frequency of the deleterious mutation is wellapproximated by the mutation–selection balance—its
frequency in the deterministic, infinite population
approximation. Under these assumptions, they showthat the waiting time for the arrival of the double mutant
is O(m2). More recently, the problem has been considered in the context of carcinogenesis: Iwasa et al.
(2004a,b) and Weissman et al. (2009) consider theEpisistasis and Substitution Rates

1147

Figure 6.—Robustness increases during stabilizing selection, and so neutral mutations tend toincrease an individual’s evolvability in a future environment. (A) Mean robustness, measured as thefraction of mutations that are neutral in the current environment, for large RNA populations under stabilizing selection. Each point (solid circle)
is a population mean, averaged across over 2000
replicate populations. (B) The y-axis measuresthe mean effect on evolvability, measured in thenovel environment, of neutral mutations, plottedagainst the mean robustness of populations. Mutant effects are averaged in 500-generation binsfrom the 15,000 generation preceding the environmental shift, and averaged from 3000 replicate simulations. N ¼ 5000, U ¼ 0.004, and s ¼ 0.1, as above.

waiting time to the arrival of the second mutant for amonomorphic, wild-type population. Those authors
obtained approximations to the rate of production of
double mutants, again neglecting back mutation andassuming small mutation rates:
1:
Up NE ; Up B >N
Under these assumptions they found that the rate ispﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2uB pðsÞuNE
pﬃﬃﬃﬃﬃ:
ON1
In our regime of interest, UpNE, UpB ¼ O(N),ﬃﬃﬃﬃthis
ﬃ pexpression gives the correct asymptotic order, O N ,but consistently underestimates the expected waitingtime given in Equation 2. Assuming that the first mutation is neutral, the analysis of Iwasa et al. (2004a,b) was
made rigorous and extended to our regime of interest in
Durrett et al. (2009), where a waiting time distribution
analogous to (4) is derived for a Moran model, again
neglecting back mutation. Our analysis additionallyjustifies the neglect of back mutations.
Our analysis shows that Kimura’s classic result does
not hold for neutral mutations that have epistaticconsequences. This result illustrates, and complicates,
a familiar quandary: How can we assign causal mechanisms to observed molecular changes? If a neutral
mutation allows a change at another site to be beneficial,then the substitution of that neutral mutation clearlycauses the later substitution. This type of epistasis
depends on the genotype–phenotype-fitness relationship, but not on the population size or mutation rate.
Whether the beneficial mutation reciprocally explains
the fixation of the prerequisite neutral change, however,depends critically on whether the adaptation arose when
the first mutation was still polymorphic, and so depends
on demography and mutation rates. The statement that‘‘neutral changes facilitate adaptation’’ does not suffice
to predict an elevated neutral substitution rate without
an analysis of the population dynamics. Together withrecent work on robustness and evolvability (Draghi et al.

2010), the analysis here demonstrates that standard
population-genetic models combined with contemporary observations about epistasis can lead to surprising
departures from traditional evolutionary theory.
These results seem to lead to a conflict in thedefinition of neutral mutations. Kimura’s elegant
formula for the substitution rate applies only to unconditionally neutral mutations. But the locations ofunconditionally neutral sites, or whether any such sites
exist at all, cannot be resolved by empirical measurements (Wagner 2005)—principally because it is impossible to test fitness effects in all genetic backgrounds. By
contrast, the common, operational definition we use
here is empirically convenient: a mutation is neutral if it
does not alter fitness in the genome and environment inwhich it arises. Even though substitutions at such sites do
not follow Kimura’s simple expectation, we can nevertheless analyze and understand their behavior (Equation 3). Moreover, other results in population genetics
suggest that it is important to study the behavior of suchsites, even though the neutrality of such sites depends onthe genetic background and may be transient. The factthat such mutations are neutral when they arise isextremely important—allowing them to reach highfrequency and subsequently alter, via epistasis, thelong-term evolutionary trajectory of a population (vanNimwegen and Crutchfield 2000; Weissman et al.
2009; Draghi et al. 2010). Unfortunately, Kimura’s
simple expression for the substitution rate may not be
useful under this revised, operational definition ofneutrality.
Although our analysis relies on complete linkage
between sites, neutral mutations may play a role in theadaptation of sexual populations as well. However, thebasis of our approach here is to classify mutations asconditionally neutral by assessing their effect with respect to a genotype’s otherwise identical, asexual parent.
This method allowed a precise comparison with Kimura’s expectation for strictly neutral sites, but it is inapplicable to sexual populations, where parents andtheir offspring may have different genotypes even in theabsence of recent mutations. Future research could

1148

J. A. Draghi, T. L. Parsons and J. B. Plotkin

develop rigorous methods to measure the contributionof neutral variation to the rates of adaption and patternsof substitution in sexual organisms.
Kimura’s result is often cited as the basis for themolecular clock, and so it is reasonable to ask if ourresults undermine that basis. On the one hand, we haveseen that adaptive evolution may fix ostensibly neutral
mutations through their epistatic interactions with
beneficial changes. However, the foundation for themolecular clock has been the belief that periods ofadaptive change of a gene are much more rare thanperiods of stabilizing selection. Whether or not this belief
has a broad empirical basis is still an active controversy(Takahata 2007; Hahn 2008), which certainly cannot
be resolved by theoretical studies such as this. However,
our results can inform the design of methods to inferselection from molecular data and contribute to resolving empirical questions surrounding the molecular
clock. Recently, theory and experiments have begun tochallenge the notion that a constant rate of molecular
evolution is indicative of neutral evolution (Wichman
et al. 2005; Barrick et al. 2009; Kryazhimskiy et al. 2009).
To understand how these studies, and experimental
evolution in general, reveal mechanisms of molecular
evolution acting in nature we must expand population
genetics to accommodate the effects of epistasis, particularly those interactions that link neutral and beneficial
changes, on the patterns of substitution.
We thank Sergey Kryazhimskiy and Michael Desai for many productive discussions about epistasis and two anonymous reviewers forhelpful critiques. J.B.P. acknowledges funding from the Burroughs
Wellcome Fund, the David and Lucile Packard Foundation, the James
S. McDonnell Foundation, the Alfred P. Sloan Foundation, theDefense Advanced Research Projects Agency (HR0011-09-1-0055),
and the U.S. National Institute of Allergy and Infectious Diseases(2U54AI057168).

LITERATURE CITED
Abramowitz, M., and I. A. Stegun, 1965 Handbook of Mathematical
Functions with Formulas, Graphs, and Mathematical Tables. Dover
Publications, New York.
Amitai, G., R. D. Gupta and D. S. Tawfik, 2007 Latent evolutionary
potentials under the neutral mutational drift of an enzyme. Hfsp
J. 1(1): 67–78.
Ancel, L. W., and W. Fontana, 2000 Plasticity, evolvability, andmodularity in RNA. J. Exp. Zool. 288(3): 242–283.
Athreya, K. B., and P. Ney, 1972 Branching processes. SpringerVerlag, Berlin.
Barrick, J. E., D. S. Yu, S. H. Yoon, H. Jeong, T. K. Oh et al.,2009 Genome evolution and adaptation in a long-term experiment with Escherichia coli. Nature 461(7268): 1243–1247.
Bastolla, U., H. Roman and M. Vendruscolo, 1999 Neutral evolution of model proteins: diffusion in sequence space and overdispersion. J. Theor. Biol. 200(1): 49–64.
Bastolla, U., M. Porto, H. Roman and M. Vendruscolo,
2002 Lack of self-averaging in neutral evolution of proteins.
Phys. Rev. Lett. 89(20): 208101.
Bastolla, U., M. Porto, H. Roman and M. Vendruscolo,
2003 Statistical properties of neutral evolution. J. Mol. Evol.
57: S103–S119.
Bergman, A., and M. Siegal, 2003 Evolutionary capacitance as a
general feature of complex gene networks. Nature 424(6948):
549–552.

Birky, C., and J. Walsh, 1988 Effects of linkage on rates of molecular evolution. Proc. Natl. Acad. Sci. USA 85(17): 6414–6418.
Bloom, J. D., and F. H. Arnold, 2009 In the light of directed evolution: pathways of adaptive protein evolution. Proc. Natl. Acad.
Sci. USA 106: 9995–10000.
Bloom, J. D., S. T. Labthavikul, C. R. Otey and F. H. Arnold,
2006 Protein stability promotes evolvability. Proc. Natl. Acad.
Sci. USA 103(15): 5869–5874.
Bloom, J. D., A. Raval and C. O. Wilke, 2007 Thermodynamics ofneutral protein evolution. Genetics 175: 255–266.
Bloom, J. D., L. I. Gong and D. Baltimore, 2010 Permissive secondary mutations enable the evolution of influenza oseltamivir
resistance. Science 328(5983): 1272–1275.
Blount, Z. D., C. Z. Borland and R. E. Lenski, 2008 Historical contingency and the evolution of a key innovation in an experimental population of Escherichia coli. Proc. Natl. Acad. Sci. USA
105(23): 7899–7906.
Bromham, L., and D. Penny, 2003 The modern molecular clock.
Nat. Rev. Genet. 4(3): 216–224.
Cambray, G., and D. Mazel, 2008 Synonymous genes explore different evolutionary landscapes. PLoS Genet. 4(11): e1000256.
Charlesworth, B., M. Morgan and D. Charlesworth, 1993 The
effect of deleterious mutations on neutral molecular variation.
Genetics 134: 1289–1303.
Ciliberti, S., O. C. Martin and A. Wagner, 2007 Innovation androbustness in complex regulatory gene networks. Proc. Natl.
Acad. Sci. USA 104(34): 13591–13596.
Cowperthwaite, M. C., and L. A. Meyers, 2007 How mutational
networks shape evolution: lessons from RNA models. Annu.
Rev. Ecol. Evol. Syst. 38: 203–230.
Cutler, D., 2000 Understanding the overdispersed molecular clock.
Genetics 154: 1403–1417.
DePristo, M., D. Weinreich and D. Hartl, 2005 Missense meanderings in sequence space: a biophysical view of protein evolution. Nat. Rev. Genet. 6(9): 678–687.
Draghi, J., and G. R. Wagner, 2008 Evolution of evolvability in adevelopmental model. Evolution 62(2): 301–315.
Draghi, J., and G. P. Wagner, 2009 The evolutionary dynamics of
evolvability in a gene network model. J. Evol. Biol. 22(3): 599–611.
Draghi, J. A., T. L. Parsons, G. P. Wagner and J. B. Plotkin,2010 Mutational robustness can facilitate adaptation. Nature
463(7279): 353–355.
Durrett, R., 2005 Probability: Theory and Examples, Ed. 3. Thomson
Brooks/Cole, Belmont, CA.
Durrett, R., D. Schmidt and J. Schweinsburg, 2009 A waitingtime problem arising from the study of multi-stage carcinogenesis. Ann. Appl. Probab. 19(2): 676–718.
Feller, W., 1957 An Introduction to Probability Theory and Its Applications, Ed. 2. Wiley, New York.
Fontana, W., and P. Schuster, 1998 Continuity in evolution: on
the nature of transitions. Science 280(5368): 1451–1455.
Gillespie, J., 1984 Molecular evolution of the mutational landscape. Evolution 38(5): 1116–1129.
Gillespie, J., 1986 Natural selection and the molecular clock. Mol.
Biol. Evol. 3(2): 138–155.
Gillespie, J., 1993 Substitution processes in molecular evolution. 1.
Uniform and clustered substitutions in a haploid model. Genetics 134: 971–981.
Gillespie, J., 2000 Genetic drift in an infinite population: the pseudohitchhiking model. Genetics 155: 909–919.
Gradshteı̆n, I. S., and I. M. Ryzhik, 2000 Table of Integrals, Series,
and Products, Ed. 6. Academic Press, San Diego.
Hahn, M. W., 2008 Toward a selection theory of molecular evolution. Evolution 62(2): 255–265.
Hughes, A. L., 2008 Near neutrality leading edge of the neutral theory of molecular evolution. Year Evol. Biol. 1133: 162–179.
Huynen, M., 1996 Exploring phenotype space through neutral evolution. J. Mol. Evol. 43(3): 165–169.
Huynen, M., P. Stadler and W. Fontana, 1996 Smoothness withinruggedness: the role of neutrality in adaptation. Proc. Natl. Acad.
Sci. USA 93(1): 397–401.
Iwasa, Y., F. Michor and M. A. Nowak, 2004a Evolutionary dynamics of invasion and escape. J. Theor. Biol. 226(2): 205–214.
Iwasa, Y., F. Michor and M. A. Nowak, 2004b Stochastic tunnels in
evolutionary dynamics. Genetics 166: 1571–1579.

Episistasis and Substitution Rates
Kashtan, N., and U. Alon, 2005 Spontaneous evolution of modularity and network motifs. Proc. Natl. Acad. Sci. USA 102(39):13773–13778.
Kimura, M., 1968 Evolutionary rate at molecular level. Nature
217(5129): 624–626.
Kimura, M., 1985 The role of compensatory neutral mutations inmolecular evolution. J. Genet. 64(1): 7–19.
Kimura, M., and T. Ota, 1971 On the rate of molecular evolution.
J. Mol. Evol. 1(1): 1–17.
Koelle, K., S. Cobey, B. Grenfell and M. Pascual, 2006 Epochal
evolution shapes the phylodynamics of interpandemic influenza
A (H3N2) in humans. Science 314(5807): 1898–1903.
Kryazhimskiy, S., G. Tkacik and J. B. Plotkin, 2009 The dynamics
of adaptation on correlated fitness landscapes. Proc. Natl. Acad.
Sci. USA 106(44): 18638–18643.
Kryazhimskiy, S., G. J. Dushoff, G. A. Bazykin, and J. B. Plotkin,2011 Prevalence of epistasis in the evolution of influenza A surface proteins. PLoS Genet. 7(2): e1001301.
Lenski, R. E., J. E. Barrick and C. Ofria, 2006 Balancing robustness and evolvability. PLoS Biol. 4(12): e428.
Meyers, L., F. Ancel and M. Lachmann, 2005 Evolution of genetic
potential. Plos Comput. Biol. 1(3): 236–243.
Nei, M., Y. Suzuki and M. Nozawa, 2010 The neutral theory of molecular evolution in the genomic era. Annu. Rev. Genomics Hum.
Genet. 11: 265–289.
Pakes, A. G., 1972 Further results on the critical Galton–Watsonprocess with immigration. J. Aust. Math. Soc. 13: 277–290.
Plotkin, J. B., J. Dushoff, M. M. Desai and H. B. Fraser, 2006 Codon
usage and selection on proteins. J. Mol. Evol. 63(5): 635–653.
Raval, A., 2007 Molecular clock on a neutral network. Phys. Rev.
Lett. 99(13): 138104.
Schuster, P., and W. Fontana, 1999 Chance and necessity in evolution: lessons from RNA. Phys. D Nonlinear Phenomena 133(1-4): 427–452.

1149

Takahata, N., 1987 On the overdispersed molecular clock. Genetics 116: 169–179.
Takahata, N., 2007 Molecular clock: an anti-neo-Darwinian legacy.
Genetics 176: 1–6.
van Nimwegen, E., 2006 Influenza escapes immunity along neutral
networks. Science 314(5807): 1884–1886.
van Nimwegen, E., J. Crutchfield and M. Huynen, 1999 Neutral
evolution of mutational robustness. Proc. Natl. Acad. Sci. USA
96(17): 9716–9720.
van Nimwegen, E., and J. P. Crutchfield, 2000 Metastable evolutionary dynamics: Crossing fitness barriers or escaping via neutral
paths? Bull. Math. Biol. 62(5): 799–848.
Wagner, A., 2005 Robustness and Evolvability in Living Systems. Princeton University Press Princeton, Princeton, NJ.
Wagner, A., 2008a Opinion: Neutralism and selectionism: a networkbased reconciliation. Nat. Rev. Genet. 9(12): 965–974.
Wagner, A., 2008b Robustness and evolvability: a paradox resolved.
Proc. Biol. Sci. 275(1630): 91–100.
Weinreich, D. M., and L. Chao, 2005 Rapid evolutionary escape bylarge populations from local fitness peaks is likely in nature. Evolution 59(6): 1175–1182.
Weissman, D. B., M. M. Desai, D. S. Fisher and M. W. Feldman,2009 The rate at which asexual populations cross fitness valleys.
Theor. Popul. Biol. 75(4): 286–300.
Wichman, H., J. Wichman and J. Bull, 2005 Adaptive molecular
evolution for 13,000 phage generations: a possible arms race.
Genetics 170: 19–31.
Wilke, C., 2004 Molecular clock in neutral protein evolution. BMC
Genet. 5: 25.

Communicating editor: H. G. Spencer

APPENDIX: DERIVATION OF EQUATION 2

In what follows, we derive the distribution for Tun, the waiting time to the arrival of the first beneficial mutation,
assuming that the initial number of adaptable individuals is zero. In previous work (Draghi et al. 2010), we obtained
the expected waiting time for arbitrary initial numbers. Here we present an alternate approach that both gives thecomplete distribution, and is more conceptually clear, though is only applicable for more limited initial conditions. In
the interest of readability, we omit the details of convergence, which may be readily supplied, and adopt an informal
tone. We remark that Durrett et al. (2009) obtained an equivalent expression for the waiting-time distribution
starting from the Moran model with one-way mutation. Their arguments are similar, although ours are simplified byour reference to known results regarding the Galton–Watson process with immigration.
We will assume a Wright–Fisher model in a population of a fixed number of individuals, N, obtaining an asymptotic
expression in the limit N / ‘. Let X Nn be the number of adaptable individuals in the nth generation, and let T Nun be thetime of arrival of the first beneficial mutation.
The transition probabilities for the number of adaptable individuals is binomially distributed,

  
N
j
N
N
P X n ¼ jjX n1 ¼ i ¼
p i ð1  p i ÞN j ;
j
where

p i ¼ Up NE

i
1
N





i
1 ð1  Up NE ÞNis the probability that a random individual chosen from the previous generation produces an adaptable offspring. We
assume weak mutation—that is, we assume uNE ¼ NUpNE is constant.
Thus, if i > N, then




i
uNE
1 1
Np i ¼ uNE 1 
i/uNE 1 i
NN
as N / ‘, and using the Poisson approximation to the binomial distribution (see e.g., Feller 1957 or Durret 2005):

1150

J. A. Draghi, T. L. Parsons and J. B. Plotkin

lim X Nn ¼ X n ;N/‘

where X n  PoissonðX n1 1 uNE Þ and convergence is in distribution (technically, in the limit we take a stopped
process, to ensure X Nn >N .) Equivalently, each adaptable individual in generation n  1 produces a Poisson(1)distributed number of offspring, while a Poisson(uNE)-distributed number of adaptable offspring are born tononadaptable parents. Thus, for small i, we may approximate X Nn by a critical Galton–Watson process withimmigration (Athreya and Ney 1972), about which much is known.
More generally, we may consider the case when the adaptable types are weakly beneficial or weakly deleterious. This
corresponds to replacing pi with
ð1 1 sÞi
ð1 1 sÞi
1 ð1  Up NE Þ;p̃i ¼ Up NE 1 N  i 1 ð1 1 sÞiN  i 1 ð1 1 sÞi
where Ns ¼ s, for s constant. Proceeding as before, we find thatN p̃i /uNE 1 i
as N / ‘; i.e., incorporating weak selection leaves our Galton–Watson approximation unchanged.
We begin with a heuristic argument that we hope provides some intuition for our exact result; intuitively, since eachindividual has probability Up B ¼ uB =N of acquiring a beneficial mutation, and such mutation has probability p(s) offixing, we would expect that the first beneficial mutation that will eventually fix arises when the total number ofadaptable individuals who have ever lived is


1
:
O
Up B pðsÞ
Since each individual in the nth generation has on average one offspring, while migration is a Poisson process withmean uNE, we have thatE½X n11 jX n  ¼ X n 1 uNE ;so
E½X n11  ¼ E½X n  1 uNE :Iterating this relation, we get E½X n  ¼ E½X 0  1 nuNE . Hence,
"
E

n
X

#
Xk ¼

k¼0

¼

n
X
k¼0
n
X

E½X k kuNE 1 ðn 1 1ÞE½X 0 

k¼1¼

nðn 1 1ÞuNE 1 ðn 1 1ÞE½X 0 ;
2which, for large values of n, is asymptotically equivalent to n2uNE. Settingn 2 uNE 

1;
Up B pðsÞ

we getpﬃﬃﬃﬃﬃN1n  pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ :uNE Up B pðsÞuNE uB pðsÞAs we show below, this underestimates the true waiting time, but is asymptotically of the correct order.
We now proceed with a derivation of the waiting-time distribution; to simplify the exposition, we omit some of thedetails of convergence. In each generation, each adaptable individual may independently mutate to a beneficial typewith probability Up B ¼ uB =N, which is destined to fix with probability p(s). Thus, the waiting time to the first beneficial
mutation that will fix isEpisistasis and Substitution Rates

"

n 
Y

1151

 N#

Xk
uB
pðsÞ
N
k¼0
h uB
Pn N iN¼ E e ð N pðsÞ1e Þ k¼0 Xk ;PfT Nun . ng ¼ E

1

where


uB
uB
e ¼ ln 1 pðsÞ  pðsÞ ¼ OðN 2 Þ;NNN

and the expectation isP
taken over all possible sample paths X Nk nk¼0 .
As we have seen, E nk¼0 X k ¼ O ðn 2 Þ; the full distribution is asymptotically characterized via Theorem 5 in Pakes(1972)Pn
k¼0 X k/Y ;n2where Y has Laplace–Stieltjes transform
rﬃﬃﬃﬃ!
llY
2uNE;¼ sechE e2pﬃﬃﬃﬃﬃ
and convergence is in distribution. Motivated by our previous heuristic, we take n ¼ N t, so that
PfT Nun .

Ppﬃﬃﬃ
pﬃﬃﬃﬃﬃ
Nt NðuB pðsÞ1N eN Þt 2
Xk =Nt 2
k¼0:
N tg ¼ E eThus, taking
TN
ﬃﬃﬃﬃﬃ ;T un ¼ lim pun
N/‘ Ni.e., rescaling time bypﬃﬃﬃﬃﬃ
N and passing to the large population limit, we haverﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ !
hiuB pðsÞ
uB pðsÞt 2 Y
2u NE¼ sech
PfT un . tg ¼ E et :2Finally, using Equation 2 in Gradshteı̆n and Ryzhik (2000, Sect. 3.512), we haveð‘

PfT un . tgdt
rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ!
ð‘uB pðsÞ
2uNE
t dt¼ sech
2
0


1
1
p
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ¼
B ; uNE ;22uB pðsÞ

E½T un  ¼

0which, using Gauss’s duplication formula for the Gamma function, G (Abramowitz and Stegun 1965),


22z1
1
Gð2zÞ ¼ pﬃﬃﬃﬃ GðzÞG z 1;2pyields
22uNE 1 GðuNE Þ2
E½T un  ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ:
2uB pðsÞ Gð2uNE Þ
Using Stirling’s approximation to estimate the Beta distribution, B, we see that

ð4Þ

1152

J. A. Draghi, T. L. Parsons and J. B. Plotkin
rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
p
1
;
11O
E½T un  ¼
uNE uB pðsÞuNE
which is, up to a constant, the result of our heuristic argument. We also note that,1
1
B ; uNE ¼
1 2ln2 1 OðuNE Þ;
2uNE
so that12ln2
1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 OðuNE Þ:
E½T un  ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2uB pðsÞuNE
2uB pðsÞ
Thus, inpthe
limit
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ as uNE / 0 we approximately recover the result in Iwasa et al. (2004a,b), but
E½T un  . 1= 2uB pðs ÞuNE for all values of uNE.

GENETICS
Supporting Information
http://www.genetics.org/cgi/content/full/genetics.110.125997/DC1

Epistasis Increases the Rate of Conditionally Neutral Substitution
in an Adapting Population

Jeremy A. Draghi, Todd L. Parsons and Joshua B. Plotkin

Copyright Ó 2011 by the Genetics Society of America
DOI: 10.1534/genetics.110.125997

2 SI

J. Draghi et al.

FIGURE S1.—Comparison of three measures of substitution rates in populations with strictly neutral mutations. The
dashed, red line tracks the mean origination rate, or the rate at which mutations arise that will fix, in 10,000-generation
simulations; the dotted line tracks the rate at which mutations fix in the population. Both rates converge at steady-state,
but the origination rate is underestimated near the end of each simulation, and the substitution rate similarly
underestimated near the beginning. Note that these transitory periods are proportional to N, as expected fromcoalescent theory. The error in the origination rate measurement can be fixed by simulating populations for longer thanthe measurement period; the solid line shows the results of measuring origination rates over 10,000 generations, but
running the simulations until a mutation fixes that originated after 10,000 generations. This last method produces anaccurate estimate of the substitution rate for the entirety of the observation period.

J. Draghi et al.

3 SI

FIGURE S2.—Deleterious mutations have a negligible effect on the substitution rate of conditionally neutral
mutations. For all genotypes, a fraction pD = 0.5 of mutations are unconditionally deleterious with selective coefficientsD. The dashed line shows the mean results from Fig. 3 for comparison. N = 5000, U = 0.004, s = 0.1 for beneficial
mutations, and pB = 0.001.

4 SI

J. Draghi et al.

FIGURE S3.—The substitution rate of conditionally neutral mutations is substantially elevated, even when more
evolvable genotypes experience a greater deleterious mutation rate. For all genotypes, a fraction pD = 0.5 of mutations
are unconditionally deleterious with selective coefficient sD = 0.2. A separate fraction pD of mutations are deleterious
in evolvable backgrounds (genotypes for which pB = 0.001) and neutral in unevolvable backgrounds (pB = 0). WhenpD > 0 conditionally neutral mutations are shaped by selection against deleterious mutations and do not behave
neutrally even under stabilizing selection. However, conditionally neutral sites still substitute much more often than
Kimura's expectation when the population is adapting. The dashed line shows the mean results from Fig. 3 forcomparison. N = 5000, U = 0.004, s = 0.1 for beneficial mutations, and pB = 0.001.

J. Draghi et al.

5 SI

FIGURE S4.—Epistatic deleterious mutations reduce but do not eliminate the increase in the substitution rate ofconditionally deleterious mutations in adapting populations. For all genotypes, a fraction pD = 0.1 of mutations areunconditionally deleterious with the selective coefficients given in the figure. Conditionally neutral mutations occur at
one-tenth this rate; pNE = 0.01. Both kinds of epistatic mutations change evolvable genotypes to unevolvable ones andvice versa. The dashed line shows the mean results from Fig. 3 for comparison. N = 5000, U = 0.004, s = 0.1 forbeneficial mutations, and pB = 0.001.