INVESTIGATIONThe Structure of Genealogies in the Presenceof Purifying Selection: A Fitness-Class Coalescent
Aleksandra M. Walczak,*,1 Lauren E. Nicolaisen,†,1 Joshua B. Plotkin,‡ and Michael M. Desai†,2*CNRS-Laboratoire de Physique Théorique de l’École Normale Supérieure, 75231 Paris Cedex 05, France, †Department ofOrganismic and Evolutionary Biology, Department of Physics, and Faculty of Arts and Sciences Center for Systems Biology, Harvard
University, Cambridge, Massachusetts 02138, and ‡Department of Biology, University of Pennsylvania, Philadelphia, Pennsylvania
19104

ABSTRACT Compared to a neutral model, purifying selection distorts the structure of genealogies and hence alters the patterns ofsampled genetic variation. Although these distortions may be common in nature, our understanding of how we expect purifying
selection to affect patterns of molecular variation remains incomplete. Genealogical approaches such as coalescent theory have provendifﬁcult to generalize to situations involving selection at many linked sites, unless selection pressures are extremely strong. Here, we
introduce an effective coalescent theory (a “ﬁtness-class coalescent”) to describe the structure of genealogies in the presence of
purifying selection at many linked sites. We use this effective theory to calculate several simple statistics describing the expected
patterns of variation in sequence data, both at the sites under selection and at linked neutral sites. Our analysis combines a descriptionof the allele frequency spectrum in the presence of purifying selection with the structured coalescent approach of Kaplan et al. (1988),
to trace the ancestry of individuals through the distribution of ﬁtnesses within the population. We also derive our results using a more
direct extension of the structured coalescent approach of Hudson and Kaplan (1994). We ﬁnd that purifying selection leads to patternsof genetic variation that are related but not identical to a neutrally evolving population in which population size has varied in a speciﬁc
way in the past.

P

URIFYING selection acting simultaneously at many
linked sites (“background selection”) can substantiallyalter the patterns of molecular variation at these sites andat linked neutral sites (Hill and Robertson 1966; Kaplan
et al. 1988; Hudson and Kaplan 1994, 1995; McVean andCharlesworth 2000; Gordo et al. 2002; O’Fallon et al. 2010;
Seger et al. 2010). In recent years, evidence from sequence
data points to the general importance of these selective
forces among many linked variants in microbial and viral
populations and on short distance scales in the genomes ofsexual organisms (Comeron et al. 2008; Hahn 2008; Seger
et al. 2010). In these situations, existing theory does not
fully explain patterns of molecular evolution (Hahn 2008).

Copyright © 2012 by the Genetics Society of America
doi: 10.1534/genetics.111.134544
Manuscript received September 9, 2011; accepted for publication November 29, 2011
Supporting information is available online at http://www.genetics.org/content/
suppl/2011/12/01/genetics.111.134544.DC1.
1
These authors contributed equally to this work.
2
Corresponding author: Departments of Organismic and Evolutionary Biology and of
Physics, and FAS Center for Systems Biology, Harvard University, 435.20 Northwest
Lab, 52 Oxford St., Cambridge, MA 02138. E-mail: mdesai@oeb.harvard.edu

It is difﬁcult to incorporate negative selection at many
linked sites into genealogical frameworks such as coalescent
theory, because these frameworks typically rely on characterizing the space of possible genealogical trees before considering the possibility of mutations at various locations onthese trees. When selection operates, the probabilities ofparticular trees cannot be deﬁned independently of themutations, and the approach breaks down (Tavare 2004;
Wakeley 2009).
Despite this difﬁculty, a number of productive approaches
have been developed to predict how negative selectioninﬂuences patterns of molecular variation and to inferselection pressures from data. Charlesworth et al. (1993)
introduced the background selection model and showedthat strong purifying selection reduces the effective population size relevant for linked neutral sites (Charlesworth
1994; Charlesworth et al. 1995). However, weaker selection
also distorts patterns of variation, in a way that cannot becompletely described by a neutral model with any effective
population size (McVean and Charlesworth 2000; Comeron
and Kreitman 2002), a phenomenon often referred to as Hill–
Robertson interference (Hill and Robertson 1966). Several

Genetics, Vol. 190, 753–779 February 2012753

theoretical frameworks have been developed to analyze thissituation. The ancestral selection graph of Neuhauser andKrone (1997) and Krone and Neuhauser (1997) provides
an elegant formal solution to the problem, but unfortunately
it requires extensive numerical calculations (Przeworski et al.
1999). These limit the intuition we can draw from thismethod and make it impractical as the basis for inference
from most modern sequence data. An alternative approachis based on the structured coalescent and views the population as subdivided into different ﬁtness classes, tracing thegenealogies of individuals as they move between classes. Thisapproach was ﬁrst introduced by Kaplan et al. (1988) andfurther developed by Hudson and Kaplan (1994, 1995) inthe case in which ﬂuctuations in the size of each ﬁtness classcan be neglected. This structured coalescent approach has
been been the basis for computational methods developedby Gordo et al. (2002), Seger et al. (2010), and Zeng andCharlesworth (2011), and analytical approaches such as
those of Barton and Etheridge (2004), Hermisson et al.
(2002), and O’Fallon et al. (2010).
In this article, we build on the structured coalescentframework by introducing the idea of a “ﬁtness-class coalescent.” Rather than considering the coalescence process inreal time, we treat each ﬁtness class as a “generation” andtrace how individuals have descended by mutations throughﬁtness classes, moving from one generation to the next by
subsequent mutations. We show that the coalescent probabilities in this ﬁtness-class coalescent can be computed using
an approach based on the Poisson random ﬁeld (PRF)
method of Sawyer and Hartl (1992) or, equivalently, canbe derived as an extension of the structured coalescent approach of Hudson and Kaplan (1994).
Our ﬁtness-class coalescent theory can be precisely
mapped to a coalescence theory in which certain quantities(e.g., coalescence times) have different meanings than in the
traditional theory. We can then invert this mapping to determine the structure of genealogies and calculate statistics
describing expected patterns of genetic variation. This approach requires certain approximations, but it also has several advantages. Most importantly, we are able to derive
relatively simple analytic expressions for coalescent probabilities and distributions of simple statistics such as heterozygosity. Consistent with earlier work, we ﬁnd that theeffects of purifying selection are broadly similar to an effective population size that changes as time recedes into thepast. Our analysis makes this intuition precise and quantitative: we can compute the exact form of this time-varying
effective population size, as deﬁned by the rate of pairwise
coalescence. We also show that this intuition has important
limitations: for example, different pairs of individuals havedifferent time-varying effective population size histories,
meaning that in principle it is possible to distinguish selection from changing population size. Our approach also
makes it possible to calculate the diversity of selected alleles
themselves, which may be important when selection is common (Williamson and Orive 2002).

754

A. M. Walczak et al.

We begin in the next section by describing the ﬁtnessclass coalescent idea that underlies our approach. We thendescribe the details of our model and analyze two ways toimplement the ﬁtness-class coalescent. The ﬁrst relies on the
PRF method of Sawyer and Hartl (1992) to describe thefrequency distribution of distinct lineages within each ﬁtnessclass. We show how this lineage structure can be used tocompute coalescence probabilities in each ﬁtness class. The
second approach is based on tracing the ancestry of individuals in the order that events occur as described by Hudson
and Kaplan (1994) and implemented numerically by Gordo
et al. (2002). We show how we can sum over all possible
ancestral paths to compute equivalent coalescence probabilities in each ﬁtness class. The two approaches provide different and complementary intuitive pictures of the process
and depend on various approximations in somewhat different ways.
After computing coalescence probabilities with both approaches, we show how these probabilities can be used toanalyze the structures of genealogies, and we calculate various statistics describing genetic variation in these populations, which we compare to numerical simulations. We thendiscuss the relationship between our results, neutral theory,and earlier work on selection, and we explore how various
approximations limit our approach. The most important ofthese approximations is that we neglect ﬂuctuations in thesize of each ﬁtness class, analogous to earlier work (Hudson
and Kaplan 1994), which restricts our analysis to the case ofstrong selection (relative to inverse population size). Thisapproximation also means that we neglect Muller’s ratchet.
We describe this and related approximations and describe
their regime of validity in the Discussion. Finally, in the appendices we explore these approximations in more detail anddescribe how they inform the relationship between our workand earlier approaches.

The Fitness-Class CoalescentIn this section, we outline the main ideas underlying ourﬁtness-class coalescent approach. We begin our analysis byconsidering the balance between mutations at many linked
sites and negative selection against the mutants, which leads
to an equilibrium distribution of ﬁtnesses within a population (Haigh 1978). We illustrate this in Figure 1, for the casein which all deleterious mutations have the same ﬁtness
cost. Each individual is characterized by the number k ofdeleterious mutations it contains. Each ﬁtness class k contains many distinct lineages, each of which arose from deleterious mutations in more-ﬁt individuals, as illustrated in
Figure 2. Neutral mutations also occur, but we consider
these later.
Hudson and Kaplan (1994) observed that individuals
move between ﬁtnesses by deleterious mutations and thatwhen two individuals are in the same ﬁtness class they could
be from the same lineage and hence coalesce. Our ﬁtnessclass coalescent exploits this observation to deﬁne anFigure 1 The distribution of the fraction of the populationin each ﬁtness class. (A) The distribution of the number ofindividuals as a function of ﬁtness, where the most beneﬁcial class is arbitrarily deﬁned to have ﬁtness 1, and eachdeleterious mutation introduces a ﬁtness disadvantage ofs. Mutations move individuals to less-ﬁt classes, and selection balances this by favoring the classes more ﬁt than
average. The shape of the depicted steady-state distribution is a result of this mutation–selection balance. The inset (B) shows the processes that lead to this balance withina given ﬁtness class.

effective genealogical process that completely bypasses theancestral process in real time. Instead, we treat each ﬁtnessclass as a generation, and we count time in deleterious
mutations: each deleterious mutation moves us from one
generation to the next. In this way, we can trace the ancestryof individuals through the ﬁtness distribution. For example,
there is some probability that two individuals chosen fromﬁtness class k are genetically identical (i.e., come from thesame lineage). If not, they each arose from mutations withinﬁtness class k 2 1. If both those mutations occurred in individuals in the same lineage in ﬁtness class k 2 1, we say thetwo individuals “coalesced” in class k 2 1. If not, they camefrom different mutations from class k 2 2, and could havecoalesced there, and so on. In this way, we can constructa ﬁtness-class coalescent tree describing the relatedness of
two individuals, as illustrated in Figure 2.
In this article we show that the probability that two
randomly chosen individuals who are currently in ﬁtness classes k and k9 coalesce in class k 2 ℓ, Pck;k9/k2ℓ ; isapproximately
Pck;k9/k2ℓ ¼

1
Ak;k9 ;
2nk2ℓ sk2ℓ ℓ

individuals are in class k 2 ℓ at the same time. In otherwords, the probability that coalescence occurs in a classequals the inverse population size of the class times thenumber of generations lineages spend together in that class.
In the following sections of this article we derive Equation 1
in the two alternative ways mentioned in the Introduction:
by explicitly considering the lineage frequency distribution
and by following the path summation method of Hudson
and Kaplan (1994) and Gordo et al. (2002).
Our approach of treating mutation events as timesteps,
and computing coalescence probabilities at each timestep,
allows us to make a precise mapping to coalescence theory
in which certain quantities have a different meaning than in
the traditional theory. In this framework, we can calculate(1)

where nk is the population size of ﬁtness class k, sk is aneffective selection pressure against these individuals, and



k9
k
k2ℓ
k2ℓ
 :
¼ 
Ak;k9
(2)
ℓ
k þ k92ℓ þ k92 kThis coalescent probability is inversely proportional to the
population size of the ﬁtness class, nk2ℓ, and the effective
selection coefﬁcient within that class, sk2ℓ, modiﬁed by the
k9
combinatoric coefﬁcient Ak;ℓ : As we will see, this has a clear
intuitive interpretation. Fitness class k 2 ℓ has size nk2ℓ, so
the coalescence probability per real generation is 1=nk2ℓ . We
will see that each lineage spends of order 1/sk2ℓ generationsin that class, so the total coalescence probability in this class
has the form ð1=nk2ℓ Þð1=sk2ℓ Þ. This is multiplied by Aℓk;k9=2;
which we will show describes the probability that the two

Figure 2 Each ﬁtness class in the population is composed of many lineages, each of which was created by a single mutation and is (in ourinﬁnite-sites model) genetically unique. In the scheme each lineage is
depicted in a different color. The arrows denote an example of the ﬁtness-class coalescence process for two individuals sampled from classes 8
and 9. These individuals came from different lineages, and these lineages
were created by mutations from different lineages within the next mostﬁt class (as shown by the arrows). The arrows trace the ancestry of thetwo individuals back through the different lineages that successively
founded each other, until they ﬁnally coalesce in the class third fromright.

Coalescent Theory with Purifying Selection755

a simple analytic expression for the probability that two
lineages sampled from particular ﬁtness classes will coalescein any other ﬁtness class. These ﬁtness-class coalescenceprobabilities allow us to explicitly calculate the structure ofgenealogies in this “mutation time.” We can then computethe distribution of any statistic describing expected sequencevariation by averaging over the ﬁtness classes our originalindividuals come from. For a statistic x that depends on
genealogies between two individuals, for example, we write
expressions of the formPðxÞ ¼XHðk; k9Þ Prob k; k9 coalesce in k 2 ℓ Pðx j k; k9; ℓÞ;(3)where H(k, k9) describes the probability that two individuals
sampled at random from the population come from classes
k and k9, respectively.
From the form of these expressions and our simple result
for the coalescence probabilities, we can immediately see
the main effect of selection on the structure of genealogies.
The discussion following Equation 1 implies that the effectof negative selection is similar to that of an effective
population size that changes as time recedes into the distantpast—i.e., some Ne(t). This intuition has been suggested byearlier work (see, e.g., Seger et al. 2010). As we will see, ouranalysis describes the precise form of Ne(t): it follows the
distribution nk2ℓ as ℓ increases further to the past, modiﬁedby the coefﬁcient Ak;k9ℓ . We will also see that this picture oftime-varying population size has limits: different pairs ofindividuals have a different Ne(t). As is clear from Equation
3, these different histories are averaged according to the
distribution H(k, k9). While it is the average Ne(t) between
pairs that determines the distribution of pairwise statistics,
this suggests that statistical power may exist in larger samples to distinguish negative selection from neutral population expansion. We explore these general conclusions of ouranalysis in detail in the Discussion.
Note that in the standard neutral coalescent, one ﬁrst
calculates the distribution of coalescence times and then imagines mutations occurring as a Poisson process throughout the
coalescent tree, with rates proportional to branch lengths. In
our ﬁtness-class coalescent, by contrast, the coalescence times
are the mutations. To avoid confusion, from here on we will tothe effective generations in our model as “steps,” and refer tothe ﬁtness-class coalescent “times” as the “steptimes.” We reserve the word time to refer to the actual coalescent time,measured in actual generations.
After determining a ﬁtness-class coalescent tree, we caninvert our mapping to determine the structure of genealogies
in real time. We will do this by calculating how the steptime inour ﬁtness-class coalescent model translates into an actual
time in generations. This will allow us to relate the distribution of branch lengths in steptimes to an actual coalescent tree
in generations. We can then treat neutral mutations as is
usually done in the standard coalescent: as a Poisson processwith probabilities proportional to branch lengths.

756

A. M. Walczak et al.

Our ﬁtness-time coalescent requires a number of approximations that limit its applicability. Most importantly, we neglect
Muller’s ratchet, and more generally we ignore the effects of
ﬂuctuations in the size of each ﬁtness class. We discuss theseapproximations in more detail below. We ﬁnd that within
a broad and biologically relevant parameter regime they leadto systematic but small corrections to our results. Despite theselimitations, our approach also has several advantages relative
to previous work. The ﬁtness-time coalescent approach makes
many otherwise difﬁcult analytic calculations tractable, allows
us to compute the diversity at the selected sites in addition tolinked neutral sites, and may offer a useful basis for practicalmethods of coalescent simulation and inference.

Model
We imagine a ﬁnite haploid population of constant size N. Eachhaploid genome has a large number of sites, which begin insome ancestral state and mutate at a constant rate. Each mutation either is neutral or confers some ﬁtness disadvantage s(where, by convention, s . 0). We assume an inﬁnite-sites
framework, so there is negligible probability that two mutations segregate simultaneously at the same site. We assume
that there is no epistasis for ﬁtness and that each deleterious
mutation carries ﬁtness cost s, so that the ﬁtness of an individual with k deleterious mutations is wk = (1 2 s)k. Since weassume that s ,, 1, we often approximate wk by 1 2 sk.
The population dynamics are assumed to follow thediffusion limit of the standard Wright–Fisher model. Thatis, we assume that deleterious mutations occur at a genome-wide rate Ud per individual per generation (with deleterious mutations assumed to be decoupled from selection).
We deﬁne ud/2 [ NUd as the per-genome scaled deleterious
mutation rate. Similarly, neutral mutations occur at a rate Un
per individual per generation, and we analogously deﬁne
un/2 [ NUn. We assume that each newly arising mutation
occurs at a site at which there are no other segregating polymorphisms in the population (the inﬁnite-sites assumption).
We focus exclusively on the case of perfect linkage, where
we imagine that all the sites we are considering are in an
asexual genome or within a short enough distance in a sexual
genome that recombination can be entirely neglected. Although our model is deﬁned for haploids, this assumptionmeans that our analysis also applies to diploid populations
provided that there is no dominance (i.e., being homozygous
for the deleterious mutation carries twice the ﬁtness cost as
being heterozygous). In this case, our model is equivalent tothat considered by Hudson and Kaplan (1994).
We believe that this is the simplest possible model that isbased on a concrete picture of mutations at individual sites thatcan describe the effects of a large number of linked negatively
selected sites on patterns of genetic variation. It is essentiallyequivalent to the model described by Charlesworth et al.
(1993) and Hudson and Kaplan (1994), which has formed
the basis for much of the analysis of background selection(Charlesworth et al. 1993; Gordo et al. 2002; Seger et al. 2010).

Our analysis will develop a ﬁtness-class coalescent theory
that involves tracing the ancestry of individuals as theychange in ﬁtness by acquiring deleterious mutations. To dothis, we need to ﬁrst understand the distribution of ﬁtnesses
within the population. Since in our model all deleterious
mutations have the same ﬁtness cost s, we can classify individuals on the basis of their Hamming class, k, relative to thewild type (which by deﬁnition has k = 0). That is, individuals in class k have k deleterious mutations more than the
most-ﬁt individuals in the population. Note that not all individuals in class k have the same set of k deleterious mutations. Furthermore, k refers only to the number of deleterious
mutations an individual has; individuals with the same k canhave different numbers of neutral mutations. We normalize
ﬁtness such that by deﬁnition all individuals in class k =0 have ﬁtness 1. Individuals in class k then have ﬁtness
1 2 ks (Figure 1).
Haigh (1978) showed that the balance between mutationand selection leads to a steady state in which the fraction of
the population in ﬁtness class k, which we call hk, is given bya Poisson distribution with mean Ud/s,
hk ¼

 
e2Ud =s Ud k
:
k!
s

(4)

This means that the average ﬁtness in the population is 1 2
Ud and that k ¼ Ud =s.
Throughout our analysis, we assume that the population
exists in this steady-state mutation–selection balance. In
particular, we neglect the fact that in a ﬁnite population
there will be ﬂuctuations around this hk. This approximationis central to our approach, and we make it in subtly different
ways in both our lineage-structure and our sum of ancestralpaths calculations of the ﬁtness-class coalescence probabilities. It will typically be valid in the bulk of the ﬁtness distribution when selection is strong (Ns .. 1); our analysis islimited to this strong selection case and breaks down whenNs ≲ 1. We discuss this approximation in more detail in the
Discussion and in Appendix B. We note that this approximation also implies that we assume that Muller’s ratchet can beneglected. We will return to the question of the importanceof Muller’s ratchet in more detail in the Discussion.
We will later need to understand the distributions oftimings, Qkk21 ðtÞ; at which an individual mutates from classk 2 1 to class k. We can calculate this by noting that theprobability that an individual in class k arose from a mutation in an individual in class k 2 1 rather than a reproduction
event from an individual in class k is
NUd hk21 þ NUd hk21 :Nhk ½1 2 Ud 2 sðk 2 kÞ(5)Substituting in the steady state values for the hk, and notingthat these mutation events are a Poisson process, we ﬁnd
2skt
:
Qk21
k ðtÞ ¼ ske(6)Note that this calculation is identical to the equivalent distribution of mutation timings computed by Gordo et al. (2002)
following the approach of Hudson and Kaplan (1994).

Lineage Structure and the Fitness-Class
Coalescence Probabilities
In general, the individuals in a particular ﬁtness class k will
not be genetically identical. Rather, there will be a numberof different lineages within this class, each lineage createdby a deleterious mutation from class k 2 1. We now consider
the structure of lineage diversity among individuals withina given ﬁtness class in the mutation–selection balance. Notethat for our purposes here, we only consider deleterious
mutations in deﬁning lineages; we consider the diversity
at neutral sites separately below.
Consider a ﬁtness class k, which has an overall frequency
hk (Figure 1B). The frequency hk is maintained by a stochastic process in which the class is constantly receiving new
individuals from class k 2 1 due to deleterious mutations.
In our inﬁnite-alleles model, each such mutation createsa lineage which is an allele that is unique within the population. Each lineage ﬂuctuates in frequency for a while before eventually dying out, perhaps after acquiring additional
mutations that found new lineages in ﬁtness class k + 1. At
any given moment, there is some frequency distribution of
lineages in each class k (see Figure 2). While the identity ofthese lineages changes over time, there is a probability distribution that at any moment there is a given frequency
distribution of lineages. In steady state, this probability distribution does not change with time.
New lineages are founded in class k at a rate uk/2, where
uk ¼ 2Nhk21 Ud :(7)

These individuals are then removed from class k at a percapita rate :sk [ 2Ud 2 s k 2 k(8)
We refer to sk as the effective selection coefﬁcient against an
allele in class k, because it is the rate at which any particularlineage in class k loses individuals, and we deﬁne
gk ¼ Nsk :(9)

Using these deﬁnitions, we can compute the steady state
probability distribution of lineages using the PRF model of
Sawyer and Hartl (1992). The essential result is that the
number of distinct lineages in class k with a frequency between a and b (in the total population) is Poisson distributedÐwith mean ab fk ðxÞdx, wherefk ðxÞ ¼uk
1 2 e22gk ð12xÞ
:
xð1 2 xÞ 1 2 e22gk(10)Note that our Poisson Random Field result implies that on
average the sum of all the frequencies of all the alleles in

Coalescent Theory with Purifying Selection

757

Ð

ﬁtness class k is simply hk ¼ 01 x fk ðxÞdx, and that theprobability that two individuals chosen at the same time at
random from ﬁtness class k both come from the same lineÐ
age is 01 dxx2 fk ðxÞ=h2k .
We note that the PRF result involves various implicit
approximations, and is valid within a speciﬁc parameter
regime. Most importantly, we neglect ﬂuctuations in thesizes of each ﬁtness class. This has two main effects. First, itmeans that we neglect the corresponding ﬂuctuations in the
distribution of lineage frequencies fk(x). Second, it means
we are implicitly neglecting the fact that, given a lineage ofsize x exists in class k, the actual hk is on average not at its
steady state value (e.g., if a high-frequency lineage exists, hk
will tend to be larger). We explain these approximations indetail in Appendix B, and describe an alternative branchingprocess formulation for the lineage structure that corrects
for the second effect described above.
The ﬁtness-class coalescent probabilities

We can now calculate the degree of relatedness between
two individuals sampled from the population. Our goal is tounderstand the probability distribution of the ﬁtness-classcoalescence steptimes for two individuals chosen at random
from the population. We begin by calculating the coalescence probability in each step.
First, imagine that by chance we pick two individuals
from the same ﬁtness class k. If the two individuals are fromthe same lineage, they coalesce within this class. In this case,they are genetically identical and the coalescence steptime is
0. If not, we want to calculate the probability they coalesce: If the lineage of individual A in classin class k 2 1, Pk;k/k21c
k was founded by a mutation from class k 2 1 a time t1 ago,
and the lineage of individual B in class k was founded by
a mutation a time t2 ago, the probability the two individuals
came from a common lineage in class k 2 1 isð
xfk21 ðxÞ y
Pck;k/k21 ¼ dt1 dt2 Qk21
G ðy/x; jt2 2 t1 jÞ:k;k ðt1 ; t2 Þ
hk21 hk21 k21(11)Here Qk21
k;k ðt1 ; t2 Þ is the joint distribution of t1 and t2, x/hk is
the probability one of the individuals came from a lineage ofsize x given that the lineage exists, fk(x) is the probability
that the lineage exists, and Gk21 ðy/x; jt2 2t1 jÞ is theprobability a lineage in class k 2 1 changes in frequency
from x to y in time |t2 2 t1| (where y could be 0, corresponding to a lineage that has already mutated back to classk 2 2 by the time the second individual mutates to classk 2 1). The forms of Q and G are described in Appendix A.
If the two individuals coalesced in this ﬁrst step, thecoalescent steptime is 1. If not (which occurs with proba), we have to consider the probability theybility 12Pk;k/k21
ccoalesce at the next step (i.e., in the mutations that tookthem from class k 2 2 to k 2 1), Pck;k/k22 , and so on.
So far we have imagined that both individuals that weoriginally selected from the population came from the sameclass k. This will not generally be true. Rather, when we pick758

A. M. Walczak et al.

two individuals at random, they will come from classes k
and k9 with probability
2hk hk9 if k 6¼ k9(12)
Hðk; k9Þ ¼
h2k
if k ¼ k9
to be
For convenience we choose k # k9. We deﬁne Pk;k9/k2ℓc
the probability that two individuals from classes k and k9
coalesce in class k 2 ℓ. Note that Pck;k9/k2ℓ ¼ 0 for ℓ , 0. Forℓ $ 0 we haveð
xfk2ℓ ðxÞ yGk2ℓ ðy/x; jt2 2 t1 jÞ
Pk;k9/k2ℓ¼ dxdydt1 dt2 Qk2ℓ:ck;k9 ðt1 ; t2 Þ
hk2ℓ
hk2ℓ
(13)

From the set of coalescence probabilities Equation 13, we
can calculate the probability distribution of coalescencesteptimes between two individuals. We describe thesesteptimes by the distribution of classes in which coalescenceoccurs; given that we pick two individuals from classes k andk9 (with k # k9 by convention) the probability that they
coalesce in class k 2 ℓ is simply
k; k9/k2ℓ
fk9
k ðℓÞ ¼ Pc

ℓ21 h
i
Y
1 2 Pck; k9/k2j :

(14)

j¼0

We note that this expression contains an implicit approximation, as described in Appendix A.
Computing the coalescence probabilities

We now have a formal structure describing the structure ofcoalescent genealogies in the presence of negative selection.
It remains, however, to evaluate the coalescent probabilities
in each step by evaluating the integrals in Equation 13. We
explain the details of this calculation in Appendix A. We ﬁnd
Pck; k9/k2ℓ ¼

1
Ak; k9 ;1 þ 2Nhk2ℓ sðk 2 ℓÞ ℓ

(15)

k9
is a numerical coefﬁcient which depends on k,where Ak;ℓk9, and ℓ but not on the population parameters,





k9
k
k2ℓ
k2ℓ
 :
Ak;k9
¼ 
ℓ
k þ k92ℓ þ k9 2 k(16)

In Figure 3 we show examples of these coalescence probabilities for different population parameters. We see that theprobability of coalescence decreases with increasing selection coefﬁcients and population size.
Equation 15 is the complete solution for coalescentprobabilities in the nonconditional approximation. This
general form for the coalescence probabilities makes intuitive sense. Nhk2ℓ is the population size of class k 2 ℓ, and
1=sðk2ℓÞ is the average number of generations that anindividual spends in class k 2 ℓ before mutating away. SinceFigure 3 Examples of the coalescenceprobabilities Pk;k9/k2ℓ
for two individualscsampled from ﬁtness classes k and k9 tocoalesce in class k 2 ℓ, shown as a function of ℓ. Here Ud/s ¼ 8, s ¼ 1023, andresults are shown for Ns ¼ 10 (dottedlines), Ns ¼ 50 (dashed lines), and Ns ¼
100 (solid lines).

the per-generation coalescent probability in a population ofsize n is proportional to 1=n; it makes sense that the
coalescent probability in class k 2 ℓ is approximately proportional to one over the population size of this class times thenumber of generations individuals spend in this class. The
additional 1 in the denominator captures the fact that theindividuals might mutate away from the class before coalescing there (which reduces the average time they spend in theclass together). The numerical factor multiplying this basic
scaling, Ak;k9ℓ , comes from the integrals over the probability
distribution of mutant timings (i.e., the dt1 and dt2 integrals).
It reﬂects the probability that the ancestors of the two individuals we are considering were both in class k 2 ℓ at thesame time, since they could not otherwise coalesce there.
From this result, we can also form an intuitive picture of
the shape of genealogies in the presence of negative
selection. We have just seen that the coalescence probability
per actual generation depends on the parameters as
1=Nhk2ℓ , where the relevant value of ℓ increases as we goback in time. Thus the structure of genealogies in the presence of negative selection is similar to having a variable
population size as we go back in time. The precise natureof this variable population size is encoded in the ﬁtness
distribution hk2ℓ. For example, if we imagine sampling two
individuals from the same below-average ﬁtness class, theprobability distribution of their genealogies is like having
a population size that initially increases and then decreases
as we look backward in time. Of course, this analogy only
goes so far. Most importantly, the coalescent steptimes are
related to the statistics describing genetic diversity in a different way from how normal coalescent times are usuallyrelated to these statistics. We return to this point in thesection on the structure of genealogies below.

A Sum of Ancestral Paths Approach
We have just computed the ﬁtness-class coalescence probabilities by considering the lineage structure within each ﬁtnessclass. Kaplan et al. (1988) proposed a somewhat different
way to look at the same problem: they considered a sample
of individuals and, without explicitly describing lineage structure, computed the relative probabilities that the next event tooccur backward in time would involve a mutation or coalescent event. For example, if two individuals are in the sameﬁtness class, the next event could be either coalescence withinthat class or a mutation event. The rates at which these eventsoccur determines their relative probabilities.
In its original form, this approach used diffusion equations
to account for ﬂuctuations in the frequencies of each ﬁtnessclass hk. Barton and Etheridge (2004) used this framework toprovide a complete solution for the effect of selection at a single site on the structure of genealogies. However, it has not
yet proven possible to solve these equations in the more general case of selection at many linked sites. Instead, Hudson
and Kaplan (1994) made progress by neglecting ﬂuctuationsin the frequencies hk, the same approximation that is central
to our approach. Using this approximation, they derived a recursion relation for the mean time to a common ancestor,
their Equation 12. Gordo et al. (2002) used this equation as
the basis for a coalescent simulation.
Recursion relations of the Hudson and Kaplan (1994)
form can be solved numerically, and have been used to generate data describing coalescent statistics, but have not yet
led to an analytic description of the structure of genealogies.
We now demonstrate that these numerical methods areequivalent to our lineage-based formalism above, by showing that the Hudson and Kaplan (1994) approach can be

Coalescent Theory with Purifying Selection759

used to derive identical analytical formulas for the ﬁtnessclass coalescent probabilities. We refer to this as a “sum ofancestral paths” approach, because it relies on summing overall possible paths of individual ancestry through the ﬁtness
distribution. The equivalence of this approach to our lineagestructure calculations means that our analytical results in thisarticle match earlier numerical and simulation results basedon the Hudson and Kaplan (1994) formulation.
In order to calculate the coalescence probabilities fora sample of two individuals, we consider the set of all possible
ancestral paths these individuals may have followed. Eachpath is represented by an ordered set of events, backward intime. These events may either be deleterious mutation events,
which move one of the ancestral lineages to the previous
ﬁtness class, or coalescence events, which merge the twoancestral lineages. In order for two individuals to coalesce inclass k 2 ℓ, each ancestral lineage must undergo a series ofdeleterious mutation events, bringing them from their initialclasses to class k 2 ℓ. The lineages must then coalesce before
any additional deleterious mutations occur. For example, inorder for two individuals sampled from class k to coalesce inclass k 2 1, the ﬁrst event, backward in time, must be a deleterious mutation. This mutation can occur in either individual. After this event, one of the ancestral lineages is still inclass k, while the other is in class k 2 1. The second event,
backward in time, must be a deleterious mutation event in theancestral lineage that remains in class k. Both ancestral lineages are now in class k 2 1. Finally, the third event must bea coalescent event. Note that there are a total of two paths,
since either individual may have been the ﬁrst to mutate.
The probability of any particular ancestral path is the
product of the probability of each event in the path. We saw
above that deleterious mutations occur in an individual inclass k at rate sk. If the two individuals are in different
classes, they are not able to coalesce as the next event. Thus
the probability of each possible event is simply:
Pð1st event is del: mut: in kjk; k9Þ ¼

sk
sk þ sk9(17)

Pð1st event is del: mut: in k9jk; k9Þ ¼

sk9:
sk þ sk9(18)

If the two individuals are in the same class, the next event
may either be a coalescent event or a deleterious mutation.
Within each class, coalescence is a neutral process that
occurs with rate 1/Nhk. Therefore, we havePð1st event is coal:jk; kÞ ¼

Pð1st event is del: mut:jk; kÞ ¼

1=ðNhk Þ
1¼
sk þ sk þ 1=ðNhk Þ1 þ 2Nhk sk(19)2sk
2Nhk sk¼:
sk þ sk þ 1=ðNhk Þ 1 þ 2Nhk sk(20)

These probabilities are analogous to those used by Gordo
et al. (2002), derived from the framework of Hudson and
Kaplan (1994).

760

A. M. Walczak et al.

Using these probabilities, we can easily calculate theprobability of any particular path. In general, in order for
two individuals sampled from classes k9 and k to coalesce inclass k 2 ℓ, the ancestral paths must consist of some order of
k9 2 k + 2ℓ events which include k9 2 k + ℓ deleterious
mutation events in the ancestral lineage that began in k9,
and ℓ deleterious mutation events in the ancestral lineage
that began in k. The path must then conclude with a ﬁnal
coalescent event. Note that there are a total of
k9 2 k þ 2ℓlpossible paths, reﬂecting the number of ways to order the
mutation events in one lineage with those in the other. To
calculate the coalescence probability, we sum the probabilities of each path that results in this particular coalescence
event.
We can carry out this sum in the general case by dividing
up the ð k92kl þ 2ℓ Þ possible paths according to whether or not
the ancestral lineages ever coexisted in each class before
class k 2 ℓ. Each case leads to a different path probability,
and these probabilities can be exactly summed. We carry out
this calculation in detail in Appendix A. We ﬁnd that to1, we haveleading order in1 þ 2 Nhk2ℓ sðk2ℓÞ

Pk;k9/k2ℓ¼
c

1
Ak;k9 ;
1 þ 2 Nhk2ℓ sðk 2 ℓÞ ℓ(21)

which exactly matches our expression for the coalescenceprobabilities in our PRF approach, Equation 15.
We note that in deriving this result, we have made thesame approximations we used in our lineage structure based
approach. Thus the results from the PRF method and thesum of ancestral paths are exactly equivalent in the regime
where they are valid. However, there are subtle differencesin the results to higher orders of the approximations, which
provide useful intuition about the process. For example, inthe sum of ancestral paths approach it is more natural tok;k9/k2ℓ;
calculate fk9
k ðℓÞ directly, without ﬁrst calculating Pc
and doing so allows us to compute certain higher-ordercorrections to the coalescence probabilities. We discussthese details of the correspondence between the approximations used in the two methods in Supporting Information, File S4.

The Structure of Genealogies and Statistics of Genetic
Diversity
We can now use the coalescence probabilities describedabove to calculate the structure of genealogies in the presenceof negative selection. We can then use these genealogies tocalculate various statistics describing the genetic diversity
within the population. We know the coalescent probabilities
in each step of our ﬁtness-class coalescent process, so in
principle we can calculate the probability of any genealogy
relating an arbitrary number of individuals using methods

analogous to those used in standard neutral coalescent theory. This would then allow us to calculate the distribution of
any statistic describing the genetic diversity among these
individuals, again using methods analogous to neutral coalescent theory.
Here we will focus on the simplest genealogical relationship: the distribution of the time to the most recent common
ancestor of two individuals, which demonstrates the mainideas in the simplest context. This allows us to calculate the
distribution of the per-site heterozygosity p. This is the only
statistic relevant to a sample of two individuals. In larger
samples, the coalescent probabilities between any pair ofsampled individuals are independent of those between any
other pair that does not share the same most recent common
ancestor, so the distribution of per-site heterozygosity we
expect within such a sample is closely related to the ensemble distribution of p we calculate here.
In our ﬁtness-class coalescent framework, it is natural toconsider diversity at the negatively selected sites separatelyfrom diversity at linked neutral sites. We focus ﬁrst on the
distribution of coalescent steptimes and pd, the per-site heterozygosity at negatively selected sites alone, ignoring neutral mutations. We will then turn to the connection between
steptimes and actual times in generations, which will enable
us to calculate the distribution of neutral diversity, including
the per-site heterozygosity at neutral sites pn. In analyzing
data, we will of course typically not know a priori whichsites are neutral and which are negatively selected. In sucha situation, we merely add up the expected diversity at neutral sites and negatively selected sites, so that the total
expected per-site heterozygosity is p = pd + pn.
Distribution of steptimes and pd

where the ﬁrst sum runs from ℓ = 0 to the largest integer lessthan or equal to the smaller of k or pd/2. Note that inpractice we only have to evaluate the sum over k from0 to a multiple of Ud/s, since H(k, k9) will be negligible forlarger k.
These results for the distributions of genealogy lengths
and of pd involve several sums. However, all the terms inthese sums are straightforward and the numerical evaluations of their values are simple and fast. In Figure 4 we show
a representative example of the predicted distribution of the
per-site heterozygosity at negatively selected sites, r(pd),compared to simulation results. We explore the signiﬁcanceof the shape of the distribution r(pd), how this distribution
depends on the parameter values, and the source of thesmall but systematic deviations between the theoretical predictions and the simulation results in the Discussion.
The relationship between steptimes and timein generations

So far we have focused on the genealogies measured insteptimes, which allowed us to calculate the distribution ofheterozygosity among negatively selected sites. We would
now like to relate the steptimes to actual times in generations. To do this, we consider the probability thata coalescence event occurred at time t, given two individuals
sampled from classes k and k9 that coalesced in class k 2 ℓ,c(t|k, k9, ℓ). We compute this distribution in File S5, andﬁndcðt j k9; k; ℓÞ ¼

pX
d 21

spd ð21Þpd 2i21

i¼0

pd 2 1ik9 þ kpdB  2sBte
2 e2sAt ;A2B(25)

where we have deﬁned A [ k9+k 2i and B [ 2ðk2ℓÞþ

We begin by imagining that we sample two individuals at
random from the same ﬁtness class k. If they coalesce in classk 2 ℓ, they each acquired ℓ different deleterious mutations toreach class k. Thus the number of negatively selected sites atwhich they will be polymorphic is twice their coalescent
steptime, pd = 2ℓ. We therefore have

1=ðNshk2ℓ Þ.

Note that when Nhk2ts(k 2 ℓ) .. 1 (the same condition
required to neglect ﬂuctuations in hk, see Appendix B), thisexpression can be simpliﬁed; we ﬁndp k9 þ k:cðt j k9; k; ℓÞ ¼ sðpd þ 1Þe2sðk9þkÞt est 21 d
pd þ 1(22)(26)

where r(pd = 2ℓ) is the probability pd = 2ℓ.
More generally, if two individuals sampled from classes k
and k9 coalesce in class k 2 ℓ, we have pd = 2ℓ + k9 2 k. Thismeans we have

However, it is important to note that while this approximation may be valid in the bulk of the distribution, it will
always fail when coalescence occurs in the zero-class, wheres(k 2 ℓ) = 0. In this case, we must use the more complex
expression Equation (25) (or in the case when the coalescence time within the 0-class can be neglected compared tothe time taken to descend from the 0-class, the simpler expression described in Equation (39) below).
Averaging over the possible values of k, k9, and ℓ, we ﬁnd
the overall distribution of actual coalescent time between
two randomly chosen individuals,rðpd ¼ 2ℓÞ ¼ fkk ðℓÞ;rðpd ¼ 2ℓ þ k9 2 kjk; k9Þ ¼ fk9
k ðℓÞ:(23)

We can average this over the distributions of k and k9 to ﬁnd
the distribution of pd among individuals sampled at random
from the population. We ﬁndrðpd Þ ¼

N
XX
ℓ

d 22ℓ
Hðk; k9 ¼ k þ pd 2 2ℓÞfk9¼kþp
ðℓÞ;
k

k¼0

(24)

cðtÞ ¼

N X
k
X X
k9$k k¼0 ℓ¼0cðtjk; k9; ℓÞfk9k ðℓÞHðk; k9Þ;Coalescent Theory with Purifying Selection

(27)

761

In Figure 5, we compare this distribution of neutral heterozygosity to simulations. We ﬁnd good general agreement tothe shape of the distribution, though there are slight systematic errors (consistent with the effects of Muller’s ratchet,which we explore further in the Discussion). Note that, likeour results for the diversity at negatively selected sites, theseresults differ dramatically from the exponential distribution
a neutral model or effective population size approximation
would predict; we describe these comparisons further in the
Discussion.
We note that an alternative way to compute neutral heterozygosity is to further extend the sum of ancestral paths
approach which we used above to provide an alternative
derivation of the coalescence probabilities. In this formulation, we do not make any connection to real times. However,this approach provides an alternative way to compute the
distribution of neutral heterozygosity, r(pn). We carry out
this computation in File S6, and show that it leads to results
identical to our analysis above.
The total heterozygosity p

Figure 4 Characteristic examples of the distribution of pd. Here N ¼ 5 ·104, s ¼ 1023: (A) Ud/s ¼ 2; (B) Ud/s ¼ 4. Theoretical predictions areshown as a solid line, simulation results as a dashed line. Simulation
results are averaged across at least 300 independent simulations for eachparameter set; shaded regions show one standard error in the simulation
results. The ﬁt to simulations is good, but we tend to slightly underestimate pd, and this tendency is worse for larger Ud/s. This is consistent withthe effects of Muller’s ratchet, which becomes more problematic as we
increase Ud/s. This systematic underestimate becomes less severe (for all
values of Ud/s) as N increases, as expected, but comprehensive simulations for much larger N are computationally prohibitive.

where the distributions H(k, k9), fk9
k ðℓÞ; and c(t|k, k9, ℓ) are
as given above. However, as we will see below, in calculating
neutral diversity we will typically ﬁnd it easier to work directly with c(t|k, k9, ℓ) rather than this unconditional distribution for c(t).
The neutral heterozygosity pn

From the distributions of real times to a common ancestor
described above, we can calculate the distribution of pn, theneutral heterozygosity. Since the neutral mutations occur asa Poisson process with rate Un, and there are a total of 2t
generations in which these mutations can occur, pn follows
a Poisson distribution with mean 2 Unt, where t is drawn
from the distribution of coalescence times, Equation 27. We
haveðN
½2Un tpd 22Un tecðtÞdt:(28)rðpd Þ ¼
pn !
0

762

A. M. Walczak et al.

To calculate the distribution of total heterozygosity p =pn + pd, we must account for the fact that pd and pn are
not independent: large pd means a large coalescent steptime
and hence makes a large pn more likely. The distribution of
pd is given by r(pd) above. Above we found c(t|k, k9, ℓ),which implies thatðN½2Un tpn 22Un tecðt j k; k9; ℓÞdt:(29)rðpn jk; k9; ℓÞ ¼
pn !
0
We can compute this integral; we ﬁndrðpn j k9; k; ℓÞ

¼

pX
d 21
i¼0

0

pd ð21Þpd 2i21



pd 2 1ik9 þ k
pd





2Un
s

B B
@
A2B
2Un
s



pn þ1 2 pn þ1 A:
2Un
þA
s

þB

2Un
s

pn

1

pn

C

(30)

Since pd = 2ℓ + k 2 k9, this implies
Xrðpn jpd Þ ¼rðpn jk; k9; ℓÞ:(31)

pd ¼k9 2 kþ2ℓ

This describes the joint distribution of selected and neutral
variation, which is of interest in situations where we know
in advance which sites are likely to be neutral and which are
selected (e.g., when analyzing the joint distribution of synonymous and nonsynonymous variation). It implies a particular relationship between the observed diversity at selectedsites and the reduction in linked neutral variation.
In many situations, however, we will not know whichalleles are selected and which are neutral. In this case, wewant to understand the distribution of total heterozygosityp, which is given byrðpÞ ¼X
pn þpd ¼prðpd Þrðpn j pd Þ:(32)

Figure 5 Characteristic examples of the distributions of pn and the real coalescent times. (A) Theoretical predictions for the distribution of pn for Ud/s ¼
2, compared to simulation results. (B) Theoretical predictions for the distribution of pn for Ud/s ¼ 4, compared to simulation results. Simulation results
are averaged across at least 300 independent simulations for each parameter set; shaded regions show one standard error in the simulation results. (C)Theoretical predictions for the distribution of real coalescence times for Ud/s ¼ 2; note that these simply mirror the distribution of pn, as expected. (D)Theoretical predictions for the distribution of real coalescence times for Ud/s ¼ 4. In A–D, N ¼ 5 · 104 and s ¼ 1023. Our theory agrees well with thesimulations, but note that, as with pd, we tend to systematically underestimate pn, and this tendency is worse for larger Ud/s. This is consistent with
Muller’s ratchet and as expected becomes more problematic for larger Ud/s. This systematic underestimate becomes less severe (for all values of Ud/s) as
we increase N, as expected, but comprehensive simulations for much larger N are computationally prohibitive.

This is no more difﬁcult to calculate than r(pn), since it
involves analogous sums. In Figure 6, we compare this predicted distribution of total heterozygosity to simulations. As
with the other aspects of heterozygosity, we ﬁnd good general agreement to the simulations, with the slight systematic
errors that are consistent with the effects of Muller’s ratchet.
The mean pairwise heterozygosity

Above we have calculated the distribution of heterozygosity
for both neutral and deleterious mutations, as well as total
heterozygosity. It is straightforward to average these results
to calculate the mean pairwise heterozygosity for bothneutral and deleterious mutations; the mean total pairwise
heterozygosity is simply the sum of these. In Figure 7 andFigure 8 we show how this mean heterozygosity depends on
population size, mutation rate, and selection strength, forneutral and deleterious mutations respectively. We see thatthe dependence of hpdi on the population size is fairly weak.
While it increases roughly linearly with N in the weak selection regime, this quickly saturates and for Ns substantiallygreater than 1 the mean heterozygosity becomes almost independent of population size. The dependence on Ud/s, bycontrast, is much stronger. The dependence of hpni on theparameters is also interesting: this depends weakly on theparameters for small N or Ud/s, but for larger N becomesroughly linear. These results make intuitive sense, particularly in light of the “mutation-time” approximation that we
introduce in the Discussion, where we discuss these ﬁguresin more detail.
Statistics in larger samples

The distributions of pn and pd described above are very
different from the distributions of heterozygosity expectedin the absence of selection. We could certainly measure the
distribution of pairwise heterozygosity from a sample ofmany individuals from a population, and use this to infer
the action of selection. However, it may also be useful tounderstand the expected distribution of other statistics describing the variation in larger samples. One statistic often
used to describe variation in larger samples is the total

Coalescent Theory with Purifying Selection763

Figure 6 Characteristic examples of the distribution of total heterozygosity p. Here N ¼ 5 · 104, s ¼ 1023: (A) Ud/s ¼ 2; (B) Ud/s ¼ 4. Theoreticalpredictions are shown as a solid line, simulation results as a dashed line.
Simulation results are averaged across at least 300 independent simulations for each parameter set; shaded regions show one standard error in
the simulation results. The ﬁt to simulations is good, but we tend toslightly underestimate p, and this tendency is worse for larger Ud/s. Thisis for the same reasons as in the distributions of pn and pd.

number of segregating sites among a sample of n individuals, Sn. Here we describe how our framework allows us tocalculate the distribution of S3; similar methods can be used
to calculate the distribution of Sn for larger n. As we will see,
it is unwieldy to calculate closed form expressions for thesequantities in our framework, so here we merely lay out a prescription for calculating S3.
We ﬁrst consider the distribution of Sd3 ; the number ofsegregating negatively selected sites among three randomly
sampled individuals. In order to calculate the probability
a sample has a particular Sd3 ; we imagine picking three
individuals at random from the population and calculate theprobability of the coalescence events that lead to that Sd3 : We
illustrate such a situation where three individuals are
sampled from classes k, k9, and k$ in Figure 9. Two of thesethree lineages coalesced in class k1. We call the steptime atwhich two of the three lineages coalesced t3 (see Figure 9).
We next need to calculate the distribution of t2, the total
steptime to common ancestry of the three individuals.
This time of course cannot be smaller than t3. Given val-

764

A. M. Walczak et al.

Figure 7 Theoretical predictions for the mean pairwise heterozygosity at
negatively selected sites, hpdi, as a function of the parameters. (A) hpdi as
a function of Ud/s for several values of Ns. In the mutation–time approximation we expect this to be linear with a slope of 2, since on average
individuals are sampled from the mean class at k ¼ Ud/s and coalesce inthe 0-class, and hence we have pd ¼ 2Ud/s. We see that as expected thisapproximation becomes more and more accurate as Ns increases. For
smaller N, there is substantial probability of coalescence in the bulk of
the ﬁtness distribution, which is greater for larger Ud/s. Thus the slope of
hpdi as a function of Ud/s decreases as Ns decreases and has a downward
curvature. (B) hpdi as a function of Ns for several values of Ud/s. We seethat as Ns becomes large, hpdi approaches 2Ud/s, again consistent withthe mutation–time approximation. As Ns decreases, coalescence within
the bulk of the ﬁtness distribution becomes more likely, and hence
hpdi decreases.

ues of t3 and t2, it is clear from Figure 9 that the total
number of segregating negatively selected sites is
Sd3 ¼ 2t2 þ t3 2ðk$2kÞ2ðk$2k9Þ:Calculating the joint distribution of t2 and t3 is tedious,
because we must sum over all possible orderings of the coalescence events, but it can be computed using either our
lineage structure method or the sum of ancestral paths approach. The basic result is analogous to our results for the
coalescence steptime between a pair of individuals: coalescence probabilities within a given class are proportional tothe inverse size of that class times the number of real generations the ancestors of given individuals typically spend inthat class, times a factor that reﬂects the time that theFigure 9 The ﬁtness-class coalescence process for three individuals, A, B,
and C, where A and B coalesced t3 steptimes ago and C coalesced withthe other two t2 steptimes ago.

here. However, these issues will be important to explore in
future work aiming to use this framework for data analysis,
and our approach here can be used as the basis for genealogical simulations. Further, since our methods allow us toquickly compute the probability of a given genealogical history and to draw a particular genealogy from the appropriate
distribution, they may provide a useful basis for importancesampling or MCMC methods to infer selection pressures fromdata.

Numerical Simulations of the Genetic Diversity
Figure 8 Theoretical predictions for the mean real coalescence time hti.
In this ﬁgure we ﬁx s ¼ 1023 and show the dependence of the mean
pairwise coalescence time on N and on Ud/s. The mean pairwise heterozygosity at neutral sites, hpni, is simply hpni ¼ 2Unhti. (A) Mean coalescence time as a function of N for various values of Ud/s. We see that hti
increases slowly with N until for large enough N the EPS approximation
applies and hti becomes linear in N. (B) Mean coalescence time as a function of Ud/s for several values of N. For large N, the dependence is roughly
linear, consistent with the EPS approximation. For smaller N, coalescence
can occur in the bulk of the ﬁtness distribution, reducing the mean coalescence time.

ancestors of sampled individuals are present in each class atthe same time.
Given a particular value of Sd3 ; there is a relationshipbetween the steptimes and actual times (analogous toEquation 25), which we could use to ﬁnd the distributionof the total number of segregating neutral sites Sn3 : Morecomplex statistics involving even larger samples can becomputed using similar methods.
However, while this analysis provides a prescription for
calculating the distribution of Sd3 and Sn3 ; it is clear that thefull distributions are opaque. In the Discussion we provide
a simple approximation for Sn in a speciﬁc parameter regime
we refer to as the “mutation-time” regime, but the complexities of the general calculation are tangential to the ideas
behind our framework, so we do not pursue them further

We compare the predictions of our ﬁtness-class coalescenceanalysis to Monte Carlo simulations of the Wright-Fisher
model. In our simulations, we consider a population of
constant size N and we keep track of the frequencies of all
genotypes over successive, discrete generations. In each generation, N individuals are sampled with replacement fromthe preceding generation, according to the standard WrightFisher multinomial sampling procedure (Ewens 2004) inwhich the chance of sampling an individual is determinedby its ﬁtness relative to the population mean ﬁtness.
In our simulations, each genotype is characterized by
the set of sites at which it harbors deleterious mutations
and the set of sites at which it harbors neutral mutations.
In each generation, a Poisson number of deleterious mutations are introduced, with mean NUd, and a Poisson numberof neutral mutations are introduced, with mean NUn; eachnew mutation is ascribed to a novel site, indexed by a random number. The mutations are distributed randomly andindependently among the individuals in the population(so that a single individual might receive multiple mutations in a given generation). The simulations record the time(in generations) at which each distinct genotype was ﬁrst
introduced.
Starting from a monomorphic population, all simulations
were run for at least ð1=sÞlnðUd =sÞ or N generations (whichever was larger), to ensure relaxation both to the steady-

Coalescent Theory with Purifying Selection

765state mutation–selection equilibrium and to the PRF equilibrium of allelic frequencies within each ﬁtness class. The ﬁnal
state of the population—i.e., the frequencies of all surviving
genotypes—was recorded at the last generation. To produce
the empirical distributions of pd and pn shown in Figure 4
and Figure 5, we averaged across at least 300 independent
populations for each parameter set.
Our simulations allow for random ﬂuctuations in thefrequencies of each ﬁtness class, and for Muller’s ratchet. In
most of the parameter regimes we explored, the ratchet
proceeded during the simulation, so that the least-loadedclass at the end of each simulation typically contained anywhere from no deleterious mutations (typical for Ud/s = 2)
to 10 (typical for Ud/s = 4). We see that despite theseeffects, our theory agrees well with the simulations, although there are small systematic errors that are consistent
with effects of the ratchet. Generally speaking these errors
increase as we increase Ud/s, but become less severe forlarger N or s. We consider these effects of Muller’s ratchet
in more detail in the Discussion.

Discussion
In recent years, both experimental studies and sequence
data have pointed to the general importance of selective
forces among many linked variants in microbial and viral
populations and on short distance scales in the genomes ofsexual organisms (Hahn 2008). Our analysis provides
a framework for understanding how one particular typeof selection—pervasive purifying (i.e., negative) selection
against deleterious mutations—affects the structure of genetic variation at the negatively selected sites themselves
and at linked neutral loci. This type of selection is presumably widespread in many populations, in which there isa selective pressure to maintain existing genotypes, andmutations away from these genotypes at a variety of lociare deleterious.
A variety of earlier work has addressed aspects of thisproblem, as described in the Introduction. The key insight of
our approach is that instead of following the true ancestralprocess, we develop a ﬁtness-class genealogical approachthat focuses on how individuals “move” through the ﬁtness
distribution. Here each mutation plays the role of a reproductive event that moves individuals through the ﬁtness
distribution, and each ﬁtness class is a generation in which
coalescence can occur with some probability. We calculate
this probability using a simple approximation based on the
PRF model of Sawyer and Hartl (1992), rather than by considering the actual reproductive process within that class. By
extending formulas originally computed by Hudson and
Kaplan (1994), we showed that these coalescent probabilities can also be computed using a summation of ancestralpaths based on the structured coalescent described by
Kaplan et al. (1988). Hence the conclusions from our analysis also describe the simulations of Gordo et al. (2002) andare consistent with all other results based on this structured

766

A. M. Walczak et al.

coalescent approach. Our work is also closely related to recent work in a continuous-ﬁtness model by O’Fallon et al.
(2010), which uses a similar framework to analyze theweak-selection regime but not the Ns .. 1 situation we
study here. We explore the relationship between our analysis and earlier work in more detail in Appendix C.
Our approach leads to simple expressions for the coalescent probability at each step in our ﬁtness-class genealogical process. This makes it a complete effective coalescent
theory: using these probabilities, we can calculate theprobability that a sample of individuals has any particularancestral relationship. Our coalescent probabilities aredifferent from those in the standard Kingman coalescent(Kingman 1982), so the structure of genealogies has a different form.
Of course, since our process is an effective rather than an
actual coalescent, the relationship between a ﬁtness-classgenealogy and the expected statistics of genetic variation
given that genealogy is different than in the standard
neutral coalescent. Given a particular genealogy measured
in steptimes, the numbers of deleterious mutations are thecoalescent times, and to calculate the statistics of neutral
variation we have to make use of the relationship between
steptimes and actual coalescence times. This contrasts with
the Kingman coalescent, where numbers of neutral mutations are typically Poisson-distributed variables with means
proportional to coalescence times (Wakeley 2009). However, we can account for these differences by starting withthe distribution of ﬁtness-class genealogies and then converting these genealogies into actual coalescence times.
In this article, we have used this ﬁtness-class approach to
calculate simple statistics describing genetic variation, inparticular the distribution of pairwise heterozygosity. This
leads to analytic expressions for the quantities of interest,
although these expressions involve sums that are most easily
calculated numerically. These are easy to compute, do not
become harder to evaluate in larger populations, and hence
are more efﬁcient to evaluate than either simulations or
calculations within the ancestral selection graph.
An intuitive picture of the structure of genealogies

The most important aspect of our analysis is not the speciﬁc
results for heterozygosity, which match the conclusions ofearlier simulations. Rather, the ﬁtness-class coalescentapproach allows us to draw several important general
conclusions about how negative selection distorts the
structure of genealogies. For two individuals drawn fromparticular ﬁtness classes, the effect of negative selection issimilar to that of an effective population size that changes
as time recedes into the past. This is consistent with suggestions from earlier work (e.g., the simulation study ofWilliamson and Orive 2002 and the work of Seger et al.
2010). However, this is not a population size that decreases
in a simple way into the past. Our analysis shows the exact
form of this time-dependent population size. Further, it isclear from our analysis that this is not the only effect ofFigure 10 Relationship between
our results and an effective population-size approximation. (A) A
typical coalescent tree in a neutral
population of constant size. The
coalescent probability per generation between a random pair ofindividuals is the inverse population size. Time runs from the pastat the top to the present at thebottom. (B) An example of a neutral coalescent tree in a population that was smaller in the pastthan the present. The population
size is shown as the width ingreen. Coalescence events are
more likely to occur when the
population size is smaller. (C)The effective population-size history for an individual experiencing purifying selection according
to our model. The individual
spends on average 1=sk generations in class k, which has a total
size Nhk. Note that pairs of individuals are sampled from different classes k (i.e., they are not all
sampled from the bottom of thispicture). Further, the coalescenceprobabilities also include a factor of A/2, which reﬂects the probability that two lineages are in the same class at the same time. (D) The historicallyvarying effective population size Ne(t) for a pair of individuals sampled from classes k and k9, as deﬁned in the text, for several values of k and k9. The
Ne(t) for two individuals sampled at random from the whole population is also shown. Here N ¼ 5 · 104, Ud/s ¼ 6, and s ¼ 1023.

negative selection on genealogies. There are two key complications. First, the statistics of genetic variation (particularly at
the deleterious sites themselves) depend on the structure ofgenealogies differently in our ﬁtness-class coalescent than in
the standard neutral coalescent. Second, the time-varying
rate of coalescence between a pair of individuals depends
on the ﬁtness classes they were sampled from. In other words,
different pairs of individuals have a different time-varying
effective population size. This suggests that genetic diversity
cannot be represented by a single time-varying effective Ne(t)for the whole population, which means that it may be possible to develop statistical tests to distinguish negative selection
from population size. All of these general intuitive conclusions about the structure of genealogies in our ﬁtness-classcoalescent are illustrated in Figure 10.
We now pause to make this intuitive picture of the shapeof typical genealogies more precise. In general the probability that two individuals will coalesce within class k has theform Pc  ðA=2Þ1=ðnk skÞ; where nk is the population size ofthat class, sk is the effective selection pressure against individuals within that class, and A is a constant that depends onwhich classes the lineages began in, but not on any of the
population parameters. We have seen that each lineage
spends on average 1/(sk) generations in class k. Thus we
can think of each individual as seeing a historical effective
population size as shown in Figure 10C: it starts in someclass k with size nk and spends 1/sk generations in that classbefore moving to class k 2 1, and so on.

If we sample two individuals, however, they will notalways be in the same class at the same time. This effectreduces the coalescence probabilities in each class, as
captured by the factor A/2. This factor is the average fraction of the 1/(sk) generations each lineage spends in class kthat the two lineages spend there together. Alternatively, we
can think of this factor as consisting of two parts: A is theprobability that the two lineages are ever in the same classat the same time, and 1/(2sk) is the average amount of timethat they coexist in the class if they coexist at all (they eachspend on average 1/(sk) generations there, but on averageoverlap for only half this time if they overlap at all). While
the two lineages are in the class at the same time, the pergeneration coalescent probability is 1/nk.
This logic implies that genealogies in the presence of
purifying selection look like neutral genealogies witha speciﬁc type of historical population-size dependence.
Imagine, for example, that we picked two individuals fromthe same ﬁtness class k. They each spend on average 1/(sk)generations in class k, and during that time they have a probability ðA=2Þð1=nk Þ per (real) generation of coalescing (thisprobability includes the fact that on average they are both in
the class simultaneously for only a fraction of the mean timeeach spends there). So roughly speaking, they have aneffective population size of Ne  2nk =Ak;kℓ¼0 for the ﬁrst 1/(sk)generations. If they fail to coalesce, they then move to classk 2 1, where they spend 1=sðk21Þ generations and havea probability ðA=2Þð1=nk21 Þ per generation of coalescing,

Coalescent Theory with Purifying Selection767and hence an effective population size Ne  2nk21 =Ak;kℓ¼1 forthis time. If they again fail to coalesce, they move to classk 2 2, and so on.
So far, this picture of a time-dependent population size israther crude, but we can make it more precise. Speciﬁcally,
we can write the coalescence probability between two
individuals sampled from class k and k9 as a function of timein generations ascðtjk; k9Þ ¼

k
X
ℓ¼0

fk9
k ðℓÞcðt j k; k9; ℓÞ:(33)

We can then deﬁne the time-dependent effective population
size between these individuals, Ne(t), as the inverse probability of coalescence at time t given that coalescence has not
yet occurred,1
cðt j k; k9Þ
¼:
Ð
Ne ðtÞ 1 2 0t cðt9 j k; k9Þdt9(34)

In other words, the Ne(t) is deﬁned as usual as the inverse of
the probability that the two individuals will coalesce at time
t given that they have not yet done so.
We illustrate this precise time-dependent population size
Ne(t) in Figure 10D. We see that for two individuals sampledfrom the same ﬁtness class, Ne(t) typically increases into therecent past and then decreases into the more distant past.
This reﬂects the fact that the two individuals are becoming
less likely to be in the same ﬁtness class in the recent past,
but that as time recedes into the distant past they are likely
to be in the highly ﬁt classes that have smaller nk. For two
individuals sampled from classes near but not identical toeach other, Ne(t) starts high and then drops before exhibiting a pattern similar to that among individuals sampledfrom the same class. This reﬂects the fact that it takes atleast a short time before the two individuals have any
chance of being in the same class. Finally, for two individuals
sampled from more distant classes, Ne(t) simply declines
into the past, both because longer ago they were more likely
to be in the same class and more likely to be in the smallclasses near the high-ﬁtness tail.
Averaging over the whole population, Figure 10D shows
the precise time-dependent population size Ne(t) for two randomly sampled individuals. This average Ne(t) initially stays
roughly constant as time recedes into the past before decreasing thereafter. For these two randomly sampled individuals,
selection is indistinguishable from this particular historicallyvarying population size. The distribution of coalescence times
between this pair of individuals looks the same as neutral
coalescent histories with this speciﬁc population-size history.
The deleterious mutation rates and selection pressures matter
only in that they determine the form of this population size
history. We note that the average Ne(t) shown in Figure 10D
implies that recent branches of genealogies will typically belonger relative to ancient branches than we would expect
under neutrality. Thus background selection will lead to an

768

A. M. Walczak et al.

excess of low-frequency variants, and hence lead to negativevalues of Tajima’s D, consistent with expectations from previous work (Charlesworth et al. 1995; Fu 1997; Gordo et al.
2002).
However, a key difference from a neutral population oftime-varying size is that, as is clear in Figure 10D, pairs ofindividuals do not typically come from the same ﬁtnessclass. Rather, they come at random from different parts of
the ﬁtness distribution, and those that come from different
places have ancestries characterized by different historicallyvarying population sizes. The total distribution of ancestry is
the sum of all of these. In other words, the genetic variation
within the population is like that in a population in whichsome individuals had one type of historical population size
history, while others had another. If we restrict ourselves topairwise statistics such as p, the average Ne(t) across pairs ofindividuals will accurately describe the genetic diversity.
However, when we consider appropriately deﬁned statistics
in larger samples, the fact that there is no single Ne(t) for thewhole population could be important. It remains an interesting question for future work to explore how to exploit thisfact to develop statistical tests to distinguish the effects ofpurifying selection from that of a historically varying effective population size.
Approximations underlying our approach

Our analysis relies on several key approximations. First, both
our lineage structure and our sum of ancestral paths
methods assume that we can neglect ﬂuctuations in the
total frequency hk of each class. Related to this approximation, we have also implicitly assumed that the probability
a lineage in class k reaches a frequency close to hk can beneglected. In Appendix B, we analyze these approximations
in detail and show that they will hold in class k wheneverNhksk .. 1. In practice, this condition will often break downin the high- and low-ﬁtness tails of the ﬁtness distribution.
Fortunately, provided it holds in the bulk of the distributionin which most individuals will be sampled (which will typically be true provided Ns .. 1), our approach is still a good
approximation. We have also made several other more technical approximations in computing the ﬁtness-class coalescent probabilities. We discuss these in detail in File S1 andFile S4.
Our ﬁnal and most important approximation is that weassume that Muller’s ratchet can be neglected. The ratchet
occurs when h0 ﬂuctuates to 0, so we can think of this approximation as an extreme aspect of neglecting ﬂuctuationsin the sizes of each ﬁtness class. This approximation cansometimes be problematic; we discuss it in detail below.
Although we have focused primarily on situations inwhich selection is weak compared to total deleterious
mutation rates, our approach is also valid regardless ofwhether s is strong or weak compared to Ud. However, whenselection is sufﬁciently strong (Ns .. 1 and Ud/s , 1), thenan effective population size approximation accurately
describes the patterns of genetic variation, as we describebelow. Thus our methods are primarily useful for situations
in which selection is weak compared to mutation rates.
Relationship with an effective population
size approximation

Charlesworth et al. (1993) considered how selection against
many linked deleterious mutations affects linked neutral diversity in a model identical to ours. These authors foundthat when selection is sufﬁciently strong, the shape of genealogies and hence the statistics of variation at linked neutral
sites are identical to the neutral case, with a reduced effective population size. We refer to this as the effective population size (EPS) approximation.
The idea behind the EPS approximation is that whenselection is strong, deleterious mutations are quickly eliminated from the population by selection. Thus if we sample
individuals from the population, they must have very
recently descended from individuals within the class of
individuals that had no deleterious mutations (the 0-class).
The EPS approximation assumes that the time for this tohappen can be neglected and that individuals never coalesce
before it does. These individuals then coalesce within the 0class as a neutral process with effective population size
equal to the size of that 0-class, which is Ne2Ud =s . Thus thegenetic diversity within the population is identical to that ina neutral population of reduced size Ne ¼ Ne2Ud =s .
The EPS approximation is valid provided that the neutral
coalescence time within the 0-class, tneut, is large comparedto the time it takes for a typical individual to have descended
from the 0-class, tdesc. We know tneut  Ne2Ud =s, and sincea typical individual comes from ﬁtness class k  Ud/s, wehave that 
XUd =s 1
1
Ud

ln
:
tdesc j¼1 jsssThis means that the EPS approximation is valid provided U
2Ud =sNse..ln d :(35)s
Because of the exponential term on the left-hand side of thisexpression, it is clear that the EPS approximation is a strongselection, weak-mutation limit. It will tend to be valid
provided that Ns . 1 and Ud , s. However, whenever Ud
becomes much larger than s, it will typically break downeven in enormous populations, as has been suggested by
Nordborg et al. (1996) and Kaiser and Charlesworth (2009).
Our analysis describes the effects of background selection
beyond the EPS approximation. We do not assume that the
coalescence time through the ﬁtness distribution is small
compared to the coalescence times within the 0-class or thatcoalescence cannot occur among individuals carrying deleterious mutations. It is precisely these two effects that leadto distortions away from the neutral expectations, making it
impossible to describe genealogies using neutral theory witha revised effective population size. Although our analysis isa generalization of the EPS approximation, it is not inconsistent with it. However, we have focused primarily on
situations in which the EPS approximation breaks down,
and coalescence times through the ﬁtness distribution arelarge compared to those in the 0-class, because this is thesituation in which our approach is most useful.
Note also that in many situations it may be the case that
there are many linked weakly selected mutations and many
linked strongly selected mutations. In such circumstances,the process we consider and the EPS approximation can act
simultaneously, each for different classes of mutations.
Imagine we had one class of mutations with ﬁtness cost s1,
which occur with mutation rate U1, where U1 , s1 andNs1 .. 1 so that the EPS approximation applies. At the sametime, imagine another class of mutations with ﬁtness cost s2,
which occur with mutation rate U2, where U2 .. s2 so that
the EPS approximation breaks down for these mutations.
In this case, the genetic diversity we expect to see will becharacteristic of our ﬁtness-class coalescent theory (with Ud= U2 and s = s2), but with a reduced effective population
size Ne ¼ Ne2U1 =s1 . In other words, the strongly selected
mutations reduce the effective population size because allindividuals are very recently descended from an individualthat had no large-effect mutations, but the coalescence timethrough the distribution of weakly selected mutations cannot be neglected.
A “mutation–time” approximation

We have seen that our analysis accounts for two effects
missing from the EPS approximation: coalescence events
outside the 0-class and the time it takes for individuals tohave descended from the 0-class. Whenever Ud/s and N areboth sufﬁciently large, the former effect can be neglected
while the latter is still important, because the number oflineages in each ﬁtness class becomes large and hence coalescence events are very unlikely to occur outside of the 0class. This leads to an approximation that we can think of asa generalization of the EPS approximation. Rather than considering primarily the diversity generated within the most-ﬁt
background, we focus instead on the diversity that accumulates while lineages move between different less-ﬁt backgrounds. Hence we term this approach a mutation–timeapproximation (MTA) for short. In this approximation, we
assume that all individuals coalesce within the 0-class, aswith the EPS approximation. However, unlike the EPS approximation, we consider the time it took for individuals todescend from the 0-class in addition to the coalescence timewithin the 0-class. This approximation is valid for large N(when even Nh1 is enormous compared to 1/s) so that coalescence always occurs in the 0-class.
In this mutation–time approximation our results becomemuch simpler and provide a useful intuitive picture of thestructure of genealogies and genetic variation. Consider thedeleterious heterozygosity pd of two individuals sampledfrom ﬁtness classes k and k9. In this approximation, these
two individuals always coalesce in the 0-class so we always

Coalescent Theory with Purifying Selection769

have pd = k + k9. Since two individuals are sampled fromclasses k and k9 with probability H(k, k9), the distribution of
pd in the population as a whole is extremely simple: we haverðpd Þ ¼

X

Hðk; k9Þ ¼ e


2Ud pd:pd ! s

22Ud =s 1

k¼pd 2k9(36)This simple approximation makes it clear why the distribution of pd looks the way it does, and explains how it
varies with Ud/s and with N, both in this mutation–timeapproximation and more generally. For large N, when coalescence outside the 0-class can be neglected, two individuals from class k and k9 have pd = k + k9. Thus the
distribution of pd has roughly the same shape as the distribution of ﬁtness within the population. The mean pd is 2Ud/s,since the average individual comes from class k = Ud/s.
Smaller and larger pd are less likely; the distribution ofﬁtness in the population has variance equal to the mean,so the variance of the distribution of pd is also roughly equal
to its mean. As N gets smaller, there is sometimes coalescence outside of the 0-class. This reduces pd given k and k9.
Hence as we reduce N, the distribution of pd shifts somewhat leftward, with a peak somewhat below 2Ud/s, andhas slightly more variance relative to the mean since there
is a less deﬁnite correspondence between k, k9, and pd.
Since pn is determined by pd, this also explains why the
distribution of pn has the peaked form we observe andhow it depends on Ud/s and N (note that for pn the coalescence time within the 0-class, which increases linearlywith N, must also be included). All of these intuitive
expectations are reﬂected in our results, as shown in Figure
4, Figure 5, Figure 7, and Figure 8. Note for example thatin Figure 4, the peak of pd is slightly below 2Ud/s (reﬂecting the ﬁnite population size) and has variance about
equal to its mean; we have veriﬁed that as N increases,
the shape of the distribution remains roughly the same,but the mean increases toward 2Ud/s and the variancedecreases slightly.
More complex statistics of sequence variation are similarly straightforward to calculate in the mutation–time approximation. When considering larger samples, the genetic
diversity is determined by the ﬁtness classes these individuals come from, which is always simple since the probability
a given individual is sampled from ﬁtness class k is just the
Poisson-distributed hk. This approximation may therefore
prove useful in developing simple and intuitive expressions
for various statistics. For example, we can use this approximation to calculate a simple expression for the distributionof the total number of segregating negatively selected sites
in a sample of size n, Sdn ; which as we have seen above isotherwise rather involved. We haver Sdn ¼ x ¼

X

hk1 hk2 . . . hkn ;(37)

k1 ;k2 ;...kn

where the sum is over sets of the ki that sum to x. We ﬁnd

770

A. M. Walczak et al.

r Sdn ¼ x ¼ e2nUd =s
1 nUd x
:
x! s

(38)

This is a distribution that is peaked around a mean value of
nUd =s, for the same reasons the distribution of pd looks as it
does. We note, however, that as we increase the sample size n,
the population size N must be even larger for this MTA to hold.
We can also calculate the distributions of actual coalescence times and hence the distributions of statistics
describing neutral diversity in the mutation–time approximation. Consider the distribution of the real coalescence
time between two individuals chosen from classes k andk9. In the mutation–time approximation where the coalescence time within the 0-class can be neglected, the actual
coalescence time isk þ k92 1
:cðt j k; k9Þ ¼ sðk þ k9Þe2sðk þ k9Þt est 21(39)

Averaging over the values of k and k9, we have
2stcðtÞ ¼ 2Ud e2st22ðUd =sÞe :(40)

The distribution of coalescence times once within the 0-class
is c0 ðtÞ ¼ ð1=Nh0 Þe2t=ðNh0 Þ . From this distribution of real
coalescence times, we can ﬁnd the distribution of neutral
heterozygosity pn in the usual way,ðN
½2Un tpn 22Un trðpn Þ ¼ecðtÞdt:(41)
pn !
0
We can immediately see that the average coalescence
time in this MTA ist

X2Ud =s 1
i¼0

si

1
þ Nh0  lnð2Ud =sÞ þ Nh0 :s

We therefore expect that the neutral heterozygosity will on
average be


2Un
2Ud
ln
þ 2Nh0 Un :(42)
hpn i ss
The ﬁrst term in this expression comes from the time todescend through the ﬁtness distribution, while the second
term comes from the time to coalesce within the 0-class. If
this latter term is large compared to the former, the EPS
approximation applies. In the opposite case in which thetime to descend through the distribution dominates, we
can see from the MTA that, as with pd, the shape of this
distribution of pn is primarily determined by the shape ofH(k, k9). In this case, the peak in hk at k = Ud/s leads toa peak in the distribution of real times and hence a peak inthe distribution of pn. The width of the distribution of pn issomewhat wider, however, since even given individuals coming from ﬁtness classes near the mean, there is a broad
distribution of possible real times, and a broad distributionof pn even given a particular real time.

This average heterozygosity would correspond to aneffective population size of


1
2Ud
(43)
Ne  lnþ Nh0 ;ss
but as we have seen this effective population size cannot
correctly describe the full distribution of pn nor its relationship to other statistics describing the genetic diversity. For
smaller values of N where the mutation–time approximation
breaks down, the average pn would be somewhat lower thanthe MTA predicts and its distribution somewhat broader.
Muller’s ratchet

We have neglected Muller’s ratchet throughout our analysis
and assumed that the ﬁtness distribution hk is ﬁxed. Yet
Muller’s ratchet will certainly occur and in some circumstances could have a signiﬁcant impact on genetic diversity(Charlesworth and Charlesworth 1997; Gordo et al. 2002;
Seger et al. 2010). Thus this is a potentially important omission from our theory. In this section we discuss some of thecomplications associated with Muller’s ratchet that are important to keep in mind when considering our approach. We
discuss the parameter regimes where neglecting Muller’sratchet should be reasonable and those where it is likely
to cause more serious problems. We provide rough estimatesof how large we expect these problems to be and suggest
a few possible ways in which future work might incorporate
Muller’s ratchet into our general framework.
Muller’s ratchet causes several related problems within
our theoretical framework. First, it causes the values of hk
to change with time and means they may not always follow
a Poisson distribution. This changes the distribution of lineage frequencies within each class, and hence changes thecoalescence probabilities. After a “click” of the ratchet, thewhole distribution hk shifts in a complicated way, eventuallyreaching a new state where it is shifted left (so the class thatwas originally at frequency hk is now at frequency hk21, andso on). In a similarly complex way, the PRF distribution of
lineage frequencies in class k shifts from fk to fk21, and so on.
This naturally changes the coalescence probabilities in eachclass. Fortunately, since the coalescence probabilities in classk are generally very similar to those in classes k + 1 or k 2 1,
this effect is unlikely to lead to major inaccuracies provided the ratchet does not click many times within a coalescent time. This is true except when we start considering
coalescence in classes close to the 0-class, where the kdependence becomes signiﬁcant. This can be thought of as
an additional problem associated with Muller’s ratchet andis associated with the fact that the ratchet shifts the wholeﬁtness distribution. This effect is easiest to see with an example: imagine we sample two individuals within the k-classand that these individuals did not coalesce before theirancestors were both in the 0-class. At the time (in the past)when these individuals’ ancestors were in the 0-class, thiscurrent 0-class might have been the 1-class or 2-class (orhigher). Thus these two individuals within the 0-class might
not coalesce until, for example, their ancestors were in what
is currently the “22-class.” This clearly means that we mightin fact have pd . 2k, which our analysis assumes is impossible. In fact, we observe precisely this effect in simulations,and it is the reason why we commonly observe systematicdeviations where the simulated values of pd are larger thanour theory predicts.
From this discussion it is clear that the key factor in
determining whether Muller’s ratchet can reasonably be
neglected is how many times the ratchet clicks in a coalescence time. We have seen above that an average individual
coalesces through the ﬁtness distribution in a time at most oforder ð1=sÞlnðUd =sÞ generations. Once within the 0-class,coalescence times are of order Ne2Ud =s. We must compare
these times to the time it takes for the ratchet to click.
The rate of the ratchet is a complex issue that has been
analyzed by Gordo and Charlesworth (2000a,b) and Kim
and Stephan (2002) in the regime where Ne2Ud =s . 1 andby Gessler (1995) in the regime where Ne2Ud =s , 1. No general analytic expressions that are valid across all parameter
regimes exist. However, provided the ratchet does not
typically move a substantial fraction of the width of theﬁtness distribution in the coalescence time of two random
individuals, it will be a small correction to pd, and neglecting
it is a reasonable ﬁrst approximation. In practice we ﬁnd in
our simulations that for the parameter regimes we consider,
pd is at most of order 2 larger than our theoretical predictions, which would correspond roughly to the effect of a single
click of the ratchet during a typical coalescence time.
The discussion above suggests a way to incorporate Muller’sratchet within our theoretical framework, albeit in an ad hocway. The ratchet shifts the distribution hk underneath the ﬁtness-class coalescent process. The details of this shift are complicated, but on average every click of the ratchet shifts the
distribution one step to the left. We can deﬁne kmin to be thenumber of deleterious mutations (relative to the optimal genotype) in the most-ﬁt individual at any given time. For thecase where Ne2Ud =s . 1, the rest of the distribution will beapproximately a Poisson distribution, but with hk replaced by
hk2kmin . Muller’s ratchet can then be thought of as a process bywhich kmin increases over time. This increase is a random process, but has some average rate, leading to an average kmin(t).
As we look backward in time during the ﬁtness-class coalescentprocess, the value of kmin is decreasing due to Muller’s ratchet.
This suggests a simple approximation: we replace the actual
value of k with an “effective” value of k that accounts for thefact that kmin decreases as we look backward in time. For eachstep through the ﬁtness distribution, we imagine that kmin has
decreased by the appropriate amount, and hence the effective
value of k in the new ﬁtness class is decreased by ,1 compared to the old ﬁtness class. When Ne2Ud =s , 1 the ratchet isan almost deterministic process, so a similar approximation
may prove useful, but in this case the distribution hk is on
average shifted from the Poisson form (Gessler 1995). To
incorporate the ratchet into our analysis in this situation,

Coalescent Theory with Purifying Selection771

we ﬁrst must recalculate the relevant coalescence probabilities given the expected average form of hk, and then carry out
the above program. These and other methods to account for
Muller’s ratchet remain an interesting topic for future work.
Despite the potential relevance of Muller’s ratchet inpractical situations, we note that it does not affect our
results in the standard coalescent limit. As is apparent fromour general expressions for the coalescence probabilities, thestructure of our ﬁtness-class coalescent theory does not depend on all three parameters N, Ud, and s independently.
Rather, it depends only on the combinations NUd and Ns.
Thus our theory makes sense in the standard limit where
NUd and Ns are held constant while we take N / N. In thislimit, Muller’s ratchet does not occur. Whether this meansthat we can neglect the ratchet for large but ﬁnite N depends
on the convergence properties of the coalescent limit. This isa difﬁcult limit to explore with simulations, because it
requires large population sizes. However, we have used simulations to verify in a few cases in which, as expected, increasing N while keeping NUd and Ns constant does notchange the predicted structure of genealogies but decreases
some of the systematic differences between theoretical predictions and the simulations that are suggestive of the effectof the ratchet. Note that while this ratchet-free limit does
not change the structure of genealogies in our ﬁtness-classcoalescent, the distribution of real coalescent times doeschange, since all real time scales are proportional to s. Thus,
as might be expected, we must also take NUn constant as
N / N if we wish neutral diversity to also remain unaffected in this limit.
Note that this ratchet-free limit, while fairly standard incoalescent theory, is somewhat different from the mutation–
time approximation we discussed above. Of course, we caneasily imagine a population that is large enough that the
mutation–time approximation applies and then take thestandard coalescent limit.
Conclusion

Our ﬁtness-class coalescent approach provides a framework inwhich we can compute distributions of genealogical structuresin situations in which many linked negatively selected sites
distort patterns of genetic variation. We have used this framework to calculate the distributions of a few simple statistics
describing sequence variation. It remains for future work to
use this ﬁtness-class coalescent approach to compute a widearray of statistics to better understand the details of how
purifying selection on many linked sites distorts patterns ofgenetic variation. The eventual goal will be to use our results
to help interpret the increasing amounts of sequence data thatseem to point to the importance of negative selection on many
linked sites.

Acknowledgments
We thank Daniel Fisher and John Wakeley for many useful discussions, which inspired our ﬁtness-class coalescent approach.

772

A. M. Walczak et al.

M.M.D. acknowledges support from the James S. McDonnell
Foundation and the Harvard Milton Fund. A.M.W. thanks the
Princeton Center for Theoretical Science at Princeton University, where she was a fellow during some of her work on thisarticle. Many of the computations in this article were run on
the Odyssey cluster supported by the FAS Sciences Division
Research Computing Group at Harvard University. L.E.N. is
supported by the Department of Defense through the National
Defense Science and Engineering Graduate Fellowship Program, and also acknowledges support from an National
Science Foundation graduate research fellowship. J.B.P.
acknowledges support from the James S. McDonnell Foundation, the Alfred P. Sloan Foundation, the David and Lucille
Packard Foundation, the Burroughs Wellcome Fund, Defense
Advanced Research Projects Agency (HR0011-05-1-0057),
and the U.S National Institute of Allergy and Infectious
Diseases (2U54AI057168).

Literature Cited
Barton, N. H., and A. M. Etheridge, 2004 The effect of selection
on genealogies. Genetics 166: 1115–1131.
Charlesworth, B., 1994 The effect of background selection against
deleterious mutations on weakly selected, linked variants.
Genet. Res. 63: 213–227.
Charlesworth, B., and D. Charlesworth, 1997 Rapid ﬁxation ofdeleterious alleles can be caused by muller’s ratchet. Genet.
Res. 70: 63–73.
Charlesworth, B., M. T. Morgan, and D. Charlesworth, 1993 The
effect of deleterious mutations on neutral molecular variation.
Genetics 134: 1289–1303.
Charlesworth, D., B. Charlesworth, and M. T. Morgan, 1995 The
pattern of neutral molecular variation under the backgroundselection model. Genetics 141: 1619–1632.
Comeron, J. M., and M. Kreitman, 2002 Population, evolutionary
and genomic consequences of interference selection. Genetics
161: 389–410.
Comeron, J. M., A. Williford, and R. M. Kliman, 2008 The Hill–
Robertson effect: evolutionary consequences of weak selection
and linkage in ﬁnite populations. Heredity 100: 19–31.
Coop, G., and R. C. Grifﬁths, 2004 Ancestral inference on genetrees under selection. Theor. Popul. Biol. 66: 219–232.
Desai, M. M., and D. S. Fisher, 2007 Beneﬁcial mutation–selection
balance and the effect of linkage on positive selection. Genetics
176: 1759–1798.
Ewens, W. J., 1972 The sampling theory of selectively neutral
alleles. Theor. Popul. Biol. 3: 87–112.
Ewens, W. J., 2004 Mathematical Population Genetics: I. Theoretical Introduction. Springer, New York.
Fu, Y., 1997 Statistical tests of neutrality of mutations against
population growth, hitchhiking and background selection. Genetics 147: 915–925.
Gessler, D. D. G., 1995 The constraints of ﬁnite size in asexual
populations and the rate of the ratchet. Genet. Res. 66: 241–253.
Gordo, I., and B. Charlesworth, 2000a The degeneration of asexual haploid populations and the speed of Muller’s ratchet. Genetics 154: 1379–1387.
Gordo, I., and B. Charlesworth, 2000b On the speed of Muller’sratchet. Genetics 156: 2137–2140.
Gordo, I., A. Navarro, and B. Charlesworth, 2002 Muller’s ratchet
and the pattern of variation at a neutral locus. Genetics 161:
835–848.

Hahn, M. W., 2008 Toward a selection theory of molecular evolution. Evolution 62: 255–265.
Haigh, J., 1978 The accumulation of deleterious genes in a population-Muller’s ratchet. Theor. Popul. Biol. 14: 251–267.
Hermisson, J., O. Redner, H. Wagner, and E. Baake, 2002 Mutationselection balance: ancestry, load, and maximum principle.
Theor. Popul. Biol. 62: 9–46.
Hill, W., and A. Robertson, 1966 The effect of linkage on limits to
artiﬁcial selection. Genet. Res. 8: 269–294.
Hudson, R., and N. Kaplan, 1994 Gene trees with backgroundselection, pp. 140–153 in Non-neutral Evolution: Theories andMolecular Data, edited by B. Golding. Chapman & Hall, NewYork.
Hudson, R., and N. Kaplan, 1995 Deleterious background selection with recombination. Genetics 141: 1605–1617.
Kaiser, V. B., and B. Charlesworth, 2009 The effects of deleterious
mutations on evolution in non-recombining genomes. Trends
Genet. 25: 9–12.
Kaplan, N., T. Darden, and R. Hudson, 1988 The coalescent process in models with selection. Genetics 120: 819–829.
Kim, Y., and W. Stephan, 2002 Recent applications of diffusion
theory to population genetics, Modern Developments in Theoretical Population Genetics: The Legacy of Gustave Malecot, edited by
M. Slatkin, and M. Veuille. Oxford University Press, Oxford.
Kimura, M., 1955 Stochastic processes and distribution of genefrequencies under natural selection. Cold Spring Harb. Symp.
Quant. Biol. 20: 33–53.
Kingman, J. F. C., 1982 The coalescent. Stochastic Process. Appl.
13: 235–248.
Krone, S. M., and C. Neuhauser, 1997 Ancestral processes withselection. Theor. Popul. Biol. 51: 210–237.

McVean, G. A. T., and B. Charlesworth, 2000 The effects of Hill–
Robertson interference between weakly selected mutations on patterns of molecular evolution and variation. Genetics 155: 929–944.
Neuhauser, C., and S. M. Krone, 1997 The genealogy of samples
in models with selection. Genetics 145: 519–534.
Nordborg, M., B. Charlesworth, and D. Charlesworth, 1996 The
effect of recombination on background selection. Genet. Res.
67: 159–174.
O’Fallon, B. D., J. Seger, and F. R. Adler, 2010 A continuous-state
coalescent and the impact of weak selection on the structure ofgene genealogies. Mol. Biol. Evol. 27: 1162–1172.
Przeworski, M., B. Charlesworth, and J. Wall, 1999 Genealogies
and weak purifying selection. Mol. Biol. Evol. 16: 246–252.
Sawyer, S. A., and D. L. Hartl, 1992 Population genetics of polymorphism and divergence. Genetics 132: 1161–1176.
Seger, J., W. A. Smith, J. J. Perry, J. Hunn, Z. A. Kaliszewska et al.,2010 Gene genealogies strongly distorted by weakly interferingmutations in constant environments. Genetics 184: 529–545.
Tavare, S., 2004 Ancestral inference in population genetics, pp.
1–188 in Lectures on Probability Theory and Statistics, Vol. 1837,
edited by J. Picard. Springer, Berlin.
Wakeley, J., 2009 Coalescent Theory: An Introduction. Roberts,
Greenwood Village, CO.
Williamson, S., and M. E. Orive, 2002 The genealogy of a sequence subject to purifying selection at multiple sites. Mol. Biol.
Evol. 19: 1376–1384.
Zeng, K., and B. Charlesworth, 2011 The joint effects of background selection and genetic recombination on local genegenealogies. Genetics 189: 251–266.
Communicating editor: N. A Rosenberg

Appendix A: The Fitness-Class Coalescent Probabilities
PRF lineage-structure approach
In the main text, we used our PRF lineage-structure approach to write an integral expression for the probability Pk;k9/k2ℓthatctwo individuals sampled from ﬁtness classes k and k9 coalesce in class k 2 ℓ, Equation 13, above. In this appendix, we
evaluate this integral to calculate the coalescent probabilities.
Equation 13 depends on the transition probability for the change in the frequency of a lineage from x to y in a time
|t1 2 t2| in class k 2 ℓ, Gk2ℓ( y / x),|t2 2 t1|). This transition probability was calculated by Kimura (1955) and canbe expressed as an inﬁnite sum of Gegenbauer polynomials. Fortunately, it appears in the context of an integralð
IG ¼ y Gk2ℓ ð y/x; jt2 2 t1 jÞdy;(44)

which is simply the average of y over Gk2ℓ. Hence this integral is given by the deterministic result for the change in thefrequency of the lineage,
IG ¼ xe2sðk2ℓÞjt2 2t1 j :(45)

Note that this deterministic solution simply reﬂects the exponential decline in frequency of a rare deleterious allele.
Substituting Equation 45 into Equation 13, we ﬁndð
x2 fk2ℓ ðxÞ 2sðk2ℓÞjt2 2t1 j
Pk;k9/k2ℓ¼ dxdt1 dt2 Qk2ℓ
ðt
;
t
Þ
e
:
1
2
c
k;k9
h2k2ℓ

(46)

The x integral can be evaluated using standard asymptotic methods; we ﬁnd
ð10

dxx 2 fk2ℓ ðxÞ [ Ixk2ℓ ¼

1:
1 þ 2Nhk2ℓ sðk 2 ℓÞ

Coalescent Theory with Purifying Selection

(47)

773

Note that this and all further expressions for Ixk2ℓ incorporate the branching process correction for ﬂuctuations in hk describedin Appendix B. Plugging in this result, we ﬁndðk2ℓ2sðk2ℓÞjt2 2t1 j
dt1 dt2 Qk2ℓ¼
I:
Pk;k9/k2ℓcxk;k9 ðt1 ; t2 Þe(48)

k2ℓ ðt ; t Þ, the joint distribution of the times at which individuals sampled
To make further progress, we must understand Qk;k9
1 2
from ﬁtness classes k and k9 originally mutated from class k 2 ℓ to class k 2 ℓ + 1. In general, t1 and t2 are not independent,
since for the two lineages to have coalesced in class k 2 ℓ they must not have coalesced in any earlier classes, which makes
them less likely to have been in those classes at the same time. In File S1, we analyze these distortions and their effects on
the coalescence probabilities. Here we make use of a simpler approximation: since the coalescence probability in each step
will turn out to be small, conditioning on not coalescing in a particular class does not shift the distribution of mutation
timings much. We therefore neglect the complications associated with the probability distributions of the mutant timings
conditional on noncoalescence. We refer to this as the nonconditional approximation and discuss its validity further inFile S1.
k2ℓ
k2ℓ
In the nonconditional approximation, the times t1 and t2 are independent, Qk2ℓ
k;k9 ðt1 ; t2 Þ ¼ Qk ðt1 ÞQk9 ðt2 Þ. We calculate
k2ℓthese distributions of mutant timings Qk ðtÞ in File S2. Plugging these in, and evaluating the integrals as described in File S3,
we ﬁnd



k9
k
ð
k2ℓ
k2ℓ
2sðk2ℓÞjt2 2t1 j
 [ Ak;k9
dt1 dt2 Qk2ℓ
¼ 
(49)
ℓ :
k;k9 ðt1 ; t2 Þ ek þ k92ℓ þ k9 2 k

Plugging this result into Equation 48, we ﬁnd Pk;k9/k2ℓ¼ Ixk2ℓ Ak;k9cℓ , the result quoted in the main text. We note that2sðk2ℓÞjt
2 2t1 je
is the probability that the ancestor of the ﬁrst individual to mutate into class k 2 ℓ is still there when the ancestor
is the probability that the ancestors of the two individuals were inof the second individual mutated into that class. Thus Ak;k9ℓclass k 2 ℓ at the same time, while Ixk2ℓ is the probability that they coalesce if so, as described in the main text.

Sum of ancestral paths approach
In the main text, we considered the probability of any particular ancestral path in the history of a sample of two
individuals. In this section, we sum over the probabilities of all possible ancestral paths to compute the ﬁtness-classcoalescence probabilities. First, we consider sampling two individuals from the same ﬁtness class k. For these two
individuals to coalesce in class k, the ﬁrst event must be a coalescent event. Using the event probabilities computed in¼ Ixk , equivalent to our earlier lineage-based result. For these individuals to coalesce in classthe main text, we ﬁnd Pk;k/kc
k 2 1, the ﬁrst event must be a deleterious mutation event. Since both individuals’ ancestral lineages are currently in classk, the probability that the ﬁrst event is a deleterious mutation event is 12Ixk . After this event, there is now one ancestral
lineage in class k 2 1 and one in class k. The next event must be a deleterious mutation in the latter, which occurs withprobability k=ð2k21Þ. Finally, the third event must be a coalescent event. This implies
fkk ð1Þ ¼ 1 2 Ixk Ixk21

k:2k 2 1(50)Note that this logic has given us an expression for the probability that the coalescent steptime is 1, fkk ð1Þ, and not the, because we have alreadyprobability of coalescence in this class given that coalescence has not yet occurred, Pk;k/k2ℓcincluded the probability that the coalescence event does not happen in class ℓ.
We can continue to extend this logic to subsequent ﬁtness classes. For example, for coalescence to occur in class k 2 2,
there are six possible paths. We can label them as AABBc, BBAAc, ABABc, ABBAc, BABAc, and BAABc, where A corresponds to a mutation in the ﬁrst individuals’ ancestral lineage, B corresponds to a mutation in the second individuals’ ancestral lineage, and c corresponds to a coalescent event. We can calculate the probability of each path. For
example,PðAABBcÞ ¼

774

A. M. Walczak et al.

1 2 Ixk
2

!




k21
k
k 2 1 k22:
I
2k 2 1 2k 2 2 2k 2 3 x(51)

The probability of path BBAAc is identical, since it has the same probabilities at each step. However, the remaining four paths
have a different probability, because the ancestral lineages exist together in the k 2 1 class at the same time. This distorts theprobability of mutations at that step, since coalescence could also have occurred. For paths of this type, we havePðABABcÞ ¼

1 2 Ixk
2

!

!


k
1 2 Ixk21
k 2 1 k22
I:
2k 2 1
2k 2 3 x
2(52)

We add up each path to ﬁnd
fkk ð2Þ ¼ Ixk22¼ Ixk22

kðk 2 1Þ2 1 2 Ixk þ 4 1 2 Ixk
4ð2k 2 1Þð2k 2 3Þ

1 2 Ixk21

3kðk 2 1Þ
22
1 2 Ixk 2 Ixk21 þ Ixk Ixk21 :2ð2k 2 1Þð2k 2 3Þ
3
3

(53)

(54)

It is informative to consider the form of this result. The Ixk22 factor is the probability that the two ancestral lineages
coalesce in class k 2 2, given that they existed in class k 2 2 at the same time. The remaining factors represent theprobability that the two ancestral lineages existed at the same time in class k 2 2. This consists of a leading-order term
½kðk21Þ=½4ð2k21Þð2k23Þ (identical to our earlier result for Akℓ¼2 ), multiplied by a correction due to the distortion in paths
from the possibility of coalescence in previous steps.
We can continue on to consider the probability of coalescence in class k 2 3. There are now a total of ð 63 Þ possible paths.
These can be split into four types, depending upon whether the two ancestral lineages coexisted in both classes k 2 1 andk 2 2 (e.g., ABABABc), in class k 2 1 only (e.g., ABAABBc), in class k 2 2 only (e.g., AABBABc), or in neither (e.g.,
AAABBBc). The probability of each type of path is identical, except for a distortion factor ð12Ixk2i Þ for each class k 2 i inwhich the two ancestral lineages were together at the same time. The probabilities can be calculated as before and summed
to yield fkk ð3Þ. Using similar logic, we can extend this approach to the situation in which two individuals are sampled fromdifferent classes, k9 and k.
In File S4, we describe the details of carrying out this summation over all possible paths to determine the coalescentprobabilities. We ﬁnd

k2ℓ
fk9
k ðℓÞ ¼ Ix




2

k9
k9 2 k þ 2ik
2ℓ 2 2i
ℓX
21
k2ℓ
i
k2ℓ 6
ℓ2i

 6


1
2
Ixk2i
4
k9 þ k
k9
2
k
þ
2ℓ
i¼0
k9 2 k þ 2ℓ
ℓ

þ

ℓX
2 2 ℓX21

k9 2 k þ 2ii
i¼0 j.i

(55)





3
2j 2 2i
2ℓ 2 2j7j2iℓ2j
Ixk2i Ixk2j 2 . . . 75;
k9 2 k þ 2ℓℓ(56)

where as always we have assumed k # k9 by convention. The form of this solution is intuitive. The factor Ixk2ℓ is theprobability of coalescence in class k 2 ℓ, given that the two ancestral lineages existed in this class at the same time. The
remaining factors reﬂect the probability that the two lineages are together in class k 2 ℓ at some point. This consists ofcalculated previously, times a correction. The correction represents thea leading-order term, which is identical to the Ak;k9ℓdistortion in the paths due to the possibility that coalescence could have occurred at previous steps. There are a total of l + 1
terms in the correction, each of which is known and calculable.
Provided that 2Nhksk .. 1, we can neglect the higher-order terms in Equation 56. This is equivalent to calculating theprobability of coalescence in a given class, without considering the possibility that coalescence events could have occurred ink;k9/k2ℓ. Neglecting these terms alsoprevious classes. Thus it converts our expression for fk9
k ðℓÞ into an expression for Pc
implicitly makes the nonconditional approximation, as we did in the PRF method, because it assumes that the fact thatcoalescence did not occur in previous classes does not distort the likelihood of taking particular paths. Making thisapproximation, we ﬁnd¼
Pk;k9/k2ℓc1
1 þ 2Nhk2ℓ sðk 2 ℓÞ

Ak;k9
ℓ ;

(57)

Coalescent Theory with Purifying Selection775which exactly matches our expression for the coalescence probabilities in the nonconditional approximation in our PRF
approach, Equation 15.
The condition 2 Nhksk .. 1 is the condition we are already assuming in treating the frequencies of each class, hk asconstant (see Appendix B). Thus the results from the PRF method and the sum of ancestral paths are exactly equivalent in theregime where they are valid. We discuss the correspondence between approximations in the sum of ancestral paths method
as compared to the PRF method in more detail in File S4.

Appendix B: Fluctuations in hK
Throughout our analysis, we have neglected ﬂuctuations in the frequencies of each ﬁtness class hk. This approximation wasnecessary to write our PRF expressions for lineage structure, fk(x), which depend on hk. Similarly, it was necessary for us tocompute the probabilities of each possible ancestral event in our sum of ancestral paths method. In this appendix, we
examine this approximation in detail and analyze its regime of validity.
Fluctuations in the ﬁtness-class frequencies affect the coalescence probability within class k in three different ways. First,
ﬂuctuations in hk21 affect the rate at which mutations enter class k. When hk21 is larger than average, more mutations occur.
Within the PRF method, this means that there will be more small lineages than the steady-state fk(x) accounts for, which
reduces the coalescence probability. In the sum of ancestral paths method, this means that the probability of mutation events
increases relative to the probability of coalescence events, which similarly reduces the coalescence probability. When hk21 issmaller than average, fewer mutations occur, and the reverse is true.
Second, ﬂuctuations in hk affect the coalescence rates within this class. Consider the case in which hk is larger than
average. Within the PRF method, this means that the probability that two individuals randomly sampled from class k comefrom a given lineage of size x is less than our assumption of x2 =h2k . This reduces the coalescence probability. In the sum ofancestral paths method, this means that the probability of coalescence events decreases relative to mutation events, whichsimilarly reduces the coalescence probability. As before, when hk is smaller than average, the reverse is true.
The third effect of ﬂuctuations is speciﬁc to the PRF method, in which we assumed that the probability that two
individuals in class k come from a lineage of frequency x (given that the lineage exists) is x2 =h2k . This implicitly assumes thatthe fact that there exists a lineage of frequency x in ﬁtness class k does not affect the expected frequency of the class hk. This
is not strictly true: given that there exists a lineage at high frequency, it is likely that hk is larger than average, and vice versa.
In other words, there is a correlation between the size of a lineage and the frequency of the class, so the probability that two
individuals picked from a class come from a lineage of frequency x is not precisely x2 =h2k. When x is large, this expression
overestimates the probability that two individuals are from the same lineage, since given that those high-frequency lineages
exist, hk will be larger than average. Similarly (although less dramatically), when x is small our expression underestimates
the probability that two individuals are from the same lineage.
Note that this third effect of ﬂuctuations is distinct from the second effect above. The second effect describes ﬂuctuations
in hk that are uncorrelated to the frequency of a particular lineage. It thus applies to both the PRF and the sum of ancestralpaths methods; it reﬂects the general fact that when hk is larger, coalescence is less likely. The third effect, on the other hand,reﬂects the fact that if we assume that we sample an individual from a lineage of size x, this biases the value of hk. Since oursum of ancestral paths method never makes any references to lineages, this third effect of ﬂuctuations applies only to the PRF
method.
These three effects all depend on the size of the ﬂuctuations relative to the average size of the each ﬁtness class. Thus
neglecting ﬂuctuations will be a good approximation provided that the ﬂuctuations in hk are small compared to hk. To
determine when this will hold, we note that each lineage in class k can reach, at most, a maximum size of order 1=sk
individuals (selection prevents any individual lineage from becoming more common than this). The total number ofindividuals in the class is on average Nhk. This means that, provided that Nhk ..1=sk , each ﬁtness class is made up of many
individual lineages. Thus we would expect that the ﬂuctuations in the sizes of each one would tend to cancel, and the overall
ﬂuctuations in hk should be negligible provided that this condition holds.
To make this intuition more precise, we must calculate the variance in hk and compare it to hk. In principle this information iscontained in our PRF expressions, but it is much simpler to compute using a continuous-time branching process method. That is,rather than use a diffusion approximation to describe the dynamics of each lineage, we use a continuous-time branching process.
As before, we imagine that new lineages in class k are created at a rate uk/2. In steady state there will be some time-independent
probability that there are n total individuals across all the lineages in the class, P(n). Note that on average we must have n/N =hk, and that P(n) contains information on the ﬂuctuations in the hk. We ﬁrst compute the generating function for P(n),HðzÞ [

N
X
n¼0

776

A. M. Walczak et al.

PðnÞzn :(58)

To do so, we start by computing the generating function for the probability distribution of the number of individuals fromeach lineage, as described by Desai and Fisher (2007, Equations 7–9). We substitute this expression into Desai and Fisher(2007, Equation 24) and integrate. We ﬁndNXu2ð12sÞ
s
Pðn; tÞz [ hz i ¼;HðzÞ [
12zð12sÞ
n¼0

n

n

(59)

where angle brackets denote expectation values, and we have suppressed the k subscripts. Note that this calculation is basedon a continuous-time branching process, in which individuals have a different distribution of offspring number than ina Wright–Fisher process, leading to a transient distribution of the frequencies of individual lineages that is half as large as in
the Wright–Fisher model for lineages of substantial frequency. Thus to make comparisons with the Wright–Fisher process, wehave to take u / 2u (as we would in comparing Wright–Fisher to Moran models), as described by Desai and Fisher (2007).
Equation 59 describes the ﬂuctuations in the size of an individual ﬁtness class: the mean, variance, and higher moments ofn can be easily computed by taking derivatives of H(z). Thus we can immediately compute Var(hk)/hk using standard
generating function methods. We ﬁnd that in fact the ﬂuctuations in hk are indeed negligible provided thatNhk sk ..1:

(60)

In practice, this condition will often break down in the high- and low-ﬁtness tails of the ﬁtness distribution. Fortunately,
provided it holds in the bulk of the distribution in which most individuals will be sampled, which will typically be true
provided Ns .. 1, our approach will still be a good approximation.
All three effects of ﬂuctuations in hk described above are negligible in the same parameter regime, Nhksk .. 1. However, the factthat the third effect applies only to our PRF result obscures the precise relationship between our two approaches and the relationshipto earlier work. Further, relaxing this approximation provides a useful comparison of the subtle differences between the assumptions
underlying the approaches. Thus we describe here an alternative approach to understanding the lineage structure in a ﬁtness classthat allows us to account for these correlations between the size of a lineage, x, and the frequency of the ﬁtness class, hk.
We ﬁrst note that, in his original calculation of the neutral ESF, Ewens (1972) used a diffusion result, f(x), roughly
analogous to our PRF expression to describe the probability that there exists a lineage with frequency x in the population at
a given time. However, Ewens’ f(x) was derived as the solution to the diffusion approximation to the K-allele Wright–Fisher
process, in the limit of inﬁnite alleles. This process explicitly imposes the constraint that the sum of all lineages in the
population at a given time must add to 1. This means that there is no correlation between the size of a lineage and the total
number of individuals in the population.
The PRF calculation of the lineage structure does not involve this explicit constraint. This is what makes it possible tocompute a simple analytical expression for fk(x). This lack of constraint means that the PRF result admits ﬂuctuations in hk,
which lead to corresponding correlations between x and hk. We could partially avoid this by deﬁning gk = Nhksk, rather thanNhk, as we have so far. This would effectively mean that each lineage is assumed to be diffusing between 0 and hk rather thanbetween 0 and 1 and forbid any lineage from reaching a frequency larger than hk. Thus it reduces the discrepancies
associated with the correlations between x and hk. However, even with this redeﬁnition, there is no constraint that thelineages in a given class all add to precisely hk, and so correlations still exist.
To correct exactly for the effects of correlations between x and hk, we extend the continuous-time branching processmodel introduced above. We now imagine that there are B sites in the genome, each of which can mutate to create a newlineage in class k. In the large-B limit, each distinct lineage in class k arose from a mutation at a different site in the genome(and we later make the inﬁnite-sites assumption B / N, which makes this exactly true). The rate at which new mutations
found lineages in class k due to mutations at a speciﬁc one of these B sites is uk =2B. This means that, analogous to Equation
59, the generating function for the probability that there are n mutations at a particular site i in class k is
Hi ðzÞ ¼

u=Bð12sÞs;
12zð12sÞ

(61)

where again we have suppressed the k subscripts and we have taken u / 2u to match to the Wright–Fisher model as
described above.
If we deﬁne ni,k to be the total number of mutants at site i in class k, we have that
sk [

B
X

ni;k

(62)

i¼1

is the total number of individuals in the class (note that on average we expect sk = Nhk). We now imagine that we sample
some number m individuals from class k. The probability that they are all from the same lineage is

Coalescent Theory with Purifying Selection

777

*
ðkÞ
Jm ¼

B nm
X
i;k
i¼1

+

sm
k

*


¼

nm
1;k

nm
2;k

n1;k þ . . . n1;B

nm
B;k

+

m þ m þ ::: þ m :n1;k þ . . . n1;B
n1;k þ ::n1;B

(63)

m
m
Note that this has the same form as our PRF expression, except that we average over nm
i =s rather than averaging over ni
and then dividing by the average sm. In other words, we explicitly account for the correlations between x and hk.
We can rewrite Equation 63 using the identity

ð N m211x¼
e2xsk dx:
sm
ðm
2
1Þ!
0
k

(64)

This identity can easily be veriﬁed by integrating the RHS by parts. Using this, and noting that lineages at each of the B sites
are independent, we ﬁnd+BÐ N xm21 2xs
P
ðkÞ
m
k
e
ni 0
dx
Jm ¼
ðm 2 1Þ!
i¼1
Ð N xm21  m 2xs (65)
n e k dx
¼B 0
ðm 2 1Þ! 1Ð N x m21
2xn1
¼B 0
dx:
he2xni iB21 nm
1e
ðm 2 1Þ!
The ﬁrst expectation value inside the integral can be computed by noting that
he2xni i ¼ Hðz ¼ 1 2 xÞ ¼ 1 þ x

12s u=½Bð12sÞ
:
s

(66)

2xn1 i. Plugging these results in and
Differentiating this result m times with respect to x results in an expression for hnm
1e
integrating, taking the limit B / N, and neglecting higher-order terms in s, we ﬁndðkÞJm ¼ um
21
X

ð21Þ j

j¼0



m21
j



1
ðm 2 1Þ!
1
:¼
¼ Qm21uþjuþm21
j¼1 ðu þ jÞu(67)

If we were to use the original PRF result to calculate the probability that two individuals sampled simultaneously fromÐðkÞclass k are from the same lineage, we would ﬁnd 01 ðhxk Þ2 fk ðxÞdx ¼ 1=u. Using our branching process result for J2 , we see that
correcting the PRF result for the third effect of ﬂuctuations in hk yields the modiﬁed probability 1=ð1 þ uk Þ. As expected, thebranching process result precisely matches the sum of ancestral paths approach, which is also unaffected by this third effectof ﬂuctuations in the hk. All of the formulae quoted in the main text and shown in the ﬁgures incorporate this correction,
which appropriately handles the correlations between the frequency of an individual lineage and the size of the ﬁtness class.

Appendix C: Relation to Previous WorkIn this appendix we compare our analysis to related work and summarize the key approximations that we and others haveused. We have presented two main approaches to calculating coalescence probabilities in this article. The ﬁrst approach isbased on the lineage structure within each ﬁtness class, described using a PRF-based method. The second approach involves
summing over all possible ancestral paths, based on the structured coalescent framework introduced by Kaplan et al. (1988)
and Hudson and Kaplan (1994, 1995). We show in this article that both approaches involve closely related approximations
and yield equivalent expressions for the coalescence probabilities.
Historically, attempts to describe the coalescent process in the presence of selection go back to the structured coalescentintroduced by Kaplan et al. (1988). These authors considered a sample of individuals from given ﬁtness classes andcomputed the relative probabilities that the next event to occur backward in time would involve a mutation or coalescentevent, without explicitly describing lineage structure. In their original work, Kaplan et al. (1988) used a full stochastic
description of the frequencies of each ﬁtness class, in which one keeps track of the probability distribution of these frequencies to account for selection. They derived diffusion equations for the transition probabilities between states. This approach isvery general, but as a result is complex and requires numerical evaluation. Barton and Etheridge (2004) developed this
diffusion approach to compute the effect of selection on genealogies in a system in which selection acts only on a single locus.

778

A. M. Walczak et al.

TABLE C1 A summary of related approaches to the coalescence process in the presence of purifying selection
This
work
Analytical expressions forgenealogy structure
Accounts for frequency classﬂuctuations (valid for Ns 1)
Valid for Nse2U/s ,, ln [U/s]
Valid for Ns .. 1
Valid for many classes
Accounts for Muller’s ratchet
Discrete ﬁtness classesabKaplan et al.
(1988)

Hudson and
Kaplan(1994, 1995)

Gordo
et al.
(2002)xCharlesworth
et al. (1993)

Barton andEtheridge(2004)

x

x

x
x
x
x
x

x
x
x
x
x

x
x
x
x

x
x
x
x

x
x
x

Seger
et al.
(2010)O’Fallon
et al.
(2010)
x

x

x

xa

x
x

x
x
x
x
x

x

xb
x

x

Addresses Ns 1 situation, but assumes deterministic ﬁtness distribution.
Within a two-class framework.

Hudson and Kaplan (1994) later simpliﬁed their original structured coalescent approach to describe the case in which
ﬂuctuations in the frequencies of ﬁtness classes can be neglected. In this deterministic approximation, they showed that one
can compute very simple expressions for the relative probabilities of the next event to occur backward in time in the historyof a sample. In this manner, Hudson and Kaplan (1994) were able to generate a simple recursion relation for the mean time
to a common ancestor, their Equation 12. Gordo et al. (2002) used this equation as the basis for a coalescent simulation, andZeng and Charlesworth (2011) recently extended this method to describe the joint effects of recombination and backgroundselection.
Recursion relations of the Hudson and Kaplan (1994) form can be solved numerically and have been used to generate
data describing coalescent statistics, but have not yet led to an analytic description of the structure of genealogies in thepresence of negative selection at many linked sites. In this article we have shown that one can sum over ancestral paths
within this framework, to derive analytical formulas for the coalescence probabilities that are equivalent to those computedfrom our lineage-based formalism. This equivalence means that our analytical results in this article match earlier numericaland simulation results based on the Hudson and Kaplan (1994) formulation. However, like the Hudson and Kaplan (1994)
framework, neither of our approaches in this article account for ﬂuctuations in the frequencies of ﬁtness classes.
In reality, the frequency of each ﬁtness class will ﬂuctuate due to genetic drift. As we have described in Appendix B, these
ﬂuctuations are substantial in classes whose deterministic size is small compared to the inverse of the effective selection
pressure against individuals in that class, Nhksk , 1. This leads to important effects on the structure of genealogies if mostﬁtness classes through the bulk of the ﬁtness distribution ﬂuctuate substantially. This will occur whenever Ns ≲ 1, so
ﬂuctuations must therefore be taken into account for small Ns. While the diffusion approach of Kaplan et al. (1988) inprinciple provides a complete solution to this problem for all values of Ns, this formalism and the related results of Barton
and Etheridge (2004) are computationally strenuous. A need for further work on accurate but more analytically tractableapproaches that are able to account for the frequency ﬂuctuations remains.
We note that the work of O’Fallon et al. (2010) and of Hermisson et al. (2002) introduced analytical approaches valid forthe case of Ns 1, although these methods are not based on a model related to the ideas of Kaplan et al. (1988). We also notethat the problem of ﬂuctuating ﬁtness-class sizes has been considered in the case of other problems (for example, forward
selection; Coop and Grifﬁths 2004), but a detailed discussion is outside the scope of this work (Table C1).
Neglecting the ﬂuctuations in ﬁtness-class frequencies is in principle reasonable when Ns .. 1. However, we note thateven when Ns .. 1, the sizes of the smallest ﬁtness classes near the tails of the distribution may still ﬂuctuate substantially.
Muller’s ratchet is one aspect of this general effect. Recently Seger et al. (2010) extended the simulation scheme of Gordo
et al. (2002) to address this problem by ﬁrst doing a forward-time simulation, recording the ﬂuctuations in the classes(including Muller’s ratchet) from this simulation, and then putting these ﬂuctuations into a backward simulation by hand.
Our methods do not account for these effects. They are therefore less general than the work of Seger et al. (2010), and break
down due to ﬂuctuation effects more quickly as Ns decreases. On the other hand, our analysis does not rely on forward
simulations and is able to compute simple analytic expressions for coalescence probabilities.
We also note that although we consider the large Ns approximation, our approach has a broader range of applicability
than the effective population-size approximation, which assumes that the coalescence time is dominated by the time tocoalescence within the most-ﬁt class. For the EPS approximation to be valid requires that this latter time ð Ne2Ud =s Þ is small
compared to the time average individuals took to descend from the most-ﬁt class ð ð1=sÞ lnNsÞ. Thus for the EPS
approximation to hold, we require Ne2Ud =s ..ð1=sÞln ½Ud =s, not just Ns .. 1. Thus we can easily have Ns .. 1, yet
Nse2Ud =s ,,ln½Ud =s, in which case the EPS approximation breaks down and yet our approach is still valid.

Coalescent Theory with Purifying Selection

779

GENETICS
Supporting Information
http://www.genetics.org/content/suppl/2011/12/01/genetics.111.134544.DC1

The Structure of Genealogies in the Presenceof Purifying Selection: A Fitness-Class Coalescent
Aleksandra M. Walczak, Lauren E. Nicolaisen, Joshua B. Plotkin, and Michael M. Desai

Copyright © 2012 by the Genetics Society of America
DOI: 10.1534/genetics.111.134544

FILE S1: THE FULL CONDITIONAL CALCULATIONIn the main text, we focused primarily on the non-conditional approximation to the coalescence probabilities,
which led to our simple expression for the coalescence probabilities, Eq. (15). In this Supplementary File,we show how this approximation can be relaxed in our lineage-structure framework by carrying out the fullconditional calculation for some of the simplest possible cases. We use this to understand the structureof the conditional results and discuss the validity of the non-conditional approximation. We note that thefull conditional result can also be obtained from the sum of ancestral paths approach by keeping the higher
order terms in Eq. (56) of Appendix A, as described in File S4, and the validity of the non-conditional
approximation can be directly assessed with that approach.
We begin by considering the full conditional result for the probability that two individuals both sampledfrom class k coalesce in class k − 2. From Appendix A of the main text, we have∫Pck,k→k−2 = Ixk−2

Qk−2
k,k (t1 , t2 ) exp [−s(k − 2)|t1 − t2 |] dt1 dt2 .

(S.1)

In order to evaluate this integral, we need to determine the probability distribution of mutant timings
Qk−2
k,k (t1 , t2 ). The time t1 is the sum of the time for one individual to have mutated from class k − 2 to classk −1 plus the time for it to have mutated from class k −1 to class k, and analogously for t2 . However, in orderfor the two lineages to coalesce in class k − 2, they must not have coalesced in class k − 1. To illustrate themain point, we neglect the distortion in the mutant timings due to the fact that individuals did not coalesce
in class k and focus only on the distortions due to the fact that coalescence did not occur in class k − 1;
if desired, the former distortion can also be included using analogous methods. We refer to the probability
distribution of the times when these individuals mutated from class k − 1 to class k conditional on them nothaving coalesced in class k − 1 as Qk−1
k,k (t1 , t2 |nc). The distribution of the times for these individuals to thenhave mutated from class k − 2 to class k − 1 is then given by
2 −s(k−1)(t1 +t2 )Qk−2.
1step (t1 , t2 ) = [s(k − 1)] e(S.2)

Thus the distribution of t1 and t2 is given by
k−1
k−2
Qk−2
k,k (t1 , t2 ) = Qk,k (t1 , t2 |nc) ⋆ Q1step (t1 , t2 ),(S.3)

where ⋆ indicates a convolution. Note that much of the time when the individuals did coalesce in class k − 1,they did so because t1 happened to be close to t2 (since this increases the chance the two individuals mutatedk−1from the same lineage). Thus in Qk−1
k,k (t1 , t2 |nc), t1 and t2 are on average further apart than in Qk,k (t1 , t2 ),and t1 and t2 are no longer independent random variables.
2 SI

A. M. Walczak et al.

We now need to calculate Qk−1
k,k (t1 , t2 |nc). We havek−1
k,k→k−1
Qk−1
k,k (t1 , t2 ) − Qk,k (t1 , t2 |c)Pc

Qk−1
k,k (t1 , t2 |nc) =1 − Pck,k→k−1

,

(S.4)

where Qk−1
k,k (t1 , t2 |c) is the distribution of timings of mutations from class k − 1 to k given that the lineages
do coalesce in class k − 1. Applying the general probability identity P (t1 , t2 |c) = P 1(c) P (c|t1 , t2 )P (t1 , t2 ),and reading oﬀ the coalescence probability given t1 and t2 from Eq. (13), we ﬁnd that
Ixk−1

Qk−1
k,k (t1 , t2 |c) =Pck,k→k−1−s(k−1)|t1 −t2 |Qk−1.
k,k (t1 , t2 )e(S.5)

We therefore ﬁnd
Qk−1
k,k (t1 , t2 |nc) =[]2 −sk(t1 +t2 )k−1
2 −2k(t1 +t2 ) −s(k−1)|t1 −t2 |
(sk)
(sk)
e
−
I
e
e
.
x
k,k→k−11

1 − Pc(S.6)

Plugging this into our convolution formula for Qk−2
k,k (t1 , t2 ) and evaluating the integrals by separating out
the possible time orderings, we ﬁnd2

Qk−2
k,k (t1 , t2 ) =k 2 [s(k − 1)]

1 − Pck,k→k−1e−s(k−1)(t1 +t2 )

[

(

1−e

−st1

)(

1−e

−2t2

)

]
Ixk−1
−
B ,
k−2

(S.7)

where we have deﬁned
B

=

[
)
1
2(
1 − e−2s min(t1 ,t2 ) −1 − e−sk min(t1 ,t2 )(k − 2)
k
(
)
(
)]
1
+
1 − e−2k|t1 −t2 | e−2s min(t1 ,t2 ) − e−sk min(t1 ,t2 ) .
k

(S.8)

We can now use this expression in Eq. (S.1) to calculate the coalescence probability Pck,k→k−2 . Since theresult is tedious and does not further illuminate the structure of the full conditional calculation, we do not
do so explicitly here, but the integrals are straightforward to evaluate with the methods we have used above.
To motivate the validity of the non-conditional approximation, we need to consider the full calculation
going back one additional step. Thus we consider the probability that two individuals both sampled fromclass k coalesce in class k − 3, Pck,k→k−3 . This will be given by∫
Pck,k→k−3 =

x2
Qk−3
(t
,
t
)
fk−3 (x)e−s(k−3)|t1 −t2 | dt1 dt2 dx,
1
2
k,k
h2k−3

(S.9)

where here Qk−3
k,k (t1 , t2 ) is the distribution of the time at which the ancestors of the two sampled individuals
originally mutated from class k − 3 to class k − 2, conditional on them not coalescing in classes k − 2 or k − 1.
k−2
We can calculate Qk−3
k,k (t1 , t2 ) in the same way we calculated Qk,k (t1 , t2 ). Explicitly,
k−2
k−3
Qk−3
k,k (t1 , t2 ) = Qk,k (t1 , t2 |nc) ⋆ Q1step (t1 , t2 ),(S.10)A. M. Walczak et al.

3 SIwhere analogously to the expression in the previous stepQk−2
k,k (t1 , t2 |nc) =1
1 − Pck,k→k−2[]
k−2
k,k→k−2Qk−2.
k,k (t1 , t2 ) − Qk,k (t1 , t2 |c)Pc(S.11)

We note that Qk−2
k,k (t1 , t2 ) is the expression in Eq. (S.7) we calculated above. As before, we havek,k→k−2−s(k−2)|t1 −t2 |Qk−2= Ixk−2 Qk−2,
k,k (t1 , t2 |c)Pc
k,k (t1 , t2 )e(S.12)[]Qk−2
k,k (t1 , t2 )k−2 −s(k−2)|t1 −t2 |
1
−
I
e
.
x
1 − Pck,k→k−2(S.13)

hence we can write
Qk−2k,k (t1 , t2 |nc) =Plugging the above expression back into Eq. (S.10), we obtainQk−3
k,k (t1 , t2 )=s2 (k − 1)2 k 2 s2 (k − 2)2−s(k−2)(t1 +t2 )∫ t2 ∫ t1ees(k−2)(y+z) es(k−1)(y+z)(1 − Pck,k→k−1 )(1 − Pck,k→k−2 )
0
0
]
[
][
I k−1
× 1 − Ixk−2 e−s(k−z)|y−z| (1 − e−sy )(1 − e−sz ) − x B .
(S.14)
k−2

We could evaluate the integrals in the above expression for Qk−3
k,k (t1 , t2 ) in the same way that we did
k−3
in our calculation for Qk−2
k,k (t1 , t2 ). We would then substitute this result for Qk,k (t1 , t2 ) into an analogous

calculation of Qk−4
k,k (t1 , t2 ), and so on. In this way we can build up the full conditional results. The mostuseful way to go about this is to separate the results into powers of Ix , which is a small parameter related tothe coalescent probability in each step. We see from the expression for Qk−3
k,k (t1 , t2 ) that there is a term in(Ix )0 , which is exactly the non-conditional approximation. There are two terms involving (Ix )1 , and a single
0
term involving (Ix )2 . In general, in the expression for Qk−ℓk,k (t1 , t2 ), we will have one (Ix ) term (which equals()
the result in the non-conditional approximation) plus ℓ terms proportional to Ix , 2ℓ terms proportional to(Ix )2 , and so on. Fortunately, the dependence on the population parameters is entirely contained withinthese powers of Ix . That is, the coeﬃcients of these various powers of Ix depend only on k and ℓ, and not
at all on the population parameters N , s, and Ud . Thus we could simply calculate a table of coeﬃcients
once, and then would be able to understand all the distributions of mutant timings (and from this all the
coalescent probabilities).
In practice, it is easier to make these full conditional calculations within the sum of ancestral paths
approach. As we show in File S4, that approach leads naturally to a power series in Ix of exactly the form
described above, in which the leading order term is the non-conditional approximation and the additionalterms represent the conditional corrections. This calculation shows that provided Ix ≪ 1, which is true
provided our usual condition that N hk sk ≫ 1 holds, these higher order terms are all small, and our nonconditional approximation is valid.
These full conditional results are, however, very complex and unilluminating. Therefore we focus here on
4 SI

A. M. Walczak et al.

understanding the general structure of these results, and on showing why the non-conditional approximationis good description of the distribution of mutation timings. We can see that at each step back through theﬁtness distribution, the probability distribution of times shifts from the non-conditional results by a factor
which is roughly proportional to the coalescence probability at that step. That is, in general we haveQk−ℓk,k (t1 , t2 ) =1
1 − Pck,k→k−ℓ][k,k→k−ℓ k−2
Qk−ℓ
Qk,k (t1 , t2 |c) .
k,k (t1 , t2 ) − Pc

(S.15)

The ﬁrst term in square brackets reﬂects the fact that the probability distribution at a given step conditional
on non-coalescence at that step is almost equal to the unconditional probability distribution at that step.
The second term represents the correction: note that it is proportional to the coalescence probability in thatstep, Pck,k→k−ℓ . The nature of the correction can be seen by plugging in the distribution of times conditionalon coalescence, giving
Qk−ℓk,k (t1 , t2 ) =

[
]
Qk−ℓ
k,k (t1 , t2 )k−ℓ −s(k−ℓ)|t1 −t2 |
e
1
−
I
.
x
1 − Pck,k→k−ℓ(S.16)We see that the correction acts to reduce the probability that |t1 − t2 | is small — that is, it makes it more
likely that t1 and t2 are further apart, because this is more likely to be the case given that coalescence did
not occur.
Since at each step the shift in the distribution of mutant timings is proportional to the coalescenceprobability, and the coalescence probability at each step is small, it seems clear that the non-conditional
approximation where we simply ignore this shift in mutant timings is reasonable. However there is one
potential caveat we must consider: although the shift in the distribution of mutation timings due to conditioning on non-coalescence is small in each step, we typically take many steps before the lineages coalesce.
In fact, since the shift in mutation timings is proportional to the coalescence probability, and we typicallygo back a number of steps of order one over the coalescence probability, in principle the shifts in mutation
timings could add up to a substantial shift.
Fortunately, there are three factors which prevent this from happening. First, the shift in mutation
1
timings at each step is always to reduce the probability of times t1 and t2 where |t1 − t2 | . (k−ℓ)s. Sinceat each step ℓ is increasing, and the range of separations between mutation timings at which coalescencecan happen is also increasing, the shifts in mutation timings from many steps ago are not a huge factor in
determining coalescence probabilities in a particular step. That is, though the shifts in mutation timings add
up over many steps, the shifts most relevant to the coalescent probability in a given step do not. Second,the coalescence probabilities at each step are diﬀerent. This reduces the chance that we take enough steps
to shift the overall mutation timings substantially by the time we coalesce. Finally, and most importantly,
we will see that the there is a substantial probability that the ancestors of the two individuals sampled doA. M. Walczak et al.

5 SInot coalesce until they are in the most-ﬁt class. This means that the total sum of coalescence probabilities(and hence the total possible weight in the shift of mutation timings) remains small even in the worst casewhere the two lineages do not coalesce for the maximum possible number of steps. The non-conditional
approximation will always be good in the regime where this is true. All of these heuristic conclusions are
reﬂected in the fact that the full conditional result we calculate in the sum of ancestral paths approach isequal to the non-conditional result plus corrections that are small provided Ix ≫ 1.

6 SI

A. M. Walczak et al.

FILE S2: THE NON-CONDITIONAL DISTRIBUTIONS OF MUTANT TIMINGS

Within the non-conditional approximation we need to calculate the distribution of mutant timings, as usedin Eq. (48). Speciﬁcally, we need to calculate
k−3
k−ℓ
k−1
Qk−ℓ
(t) ⋆ Qk−2
k (t) = Qk
k−1 (t) ⋆ Qk−2 (t) ⋆ . . . ⋆ Qk−ℓ+1 (t),(S.17)

where ⋆ refers to a convolution and−s(k−ℓ+1)tQk−ℓ,k−ℓ+1 (t) = s(k − ℓ + 1)e(S.18)

as given by Eq. (6). In general, the convolution of n exponential distributions with parameters λ1 . . . λn is
given by

n−1
∑

n−1
∏

λi e−λi t

i=0

j=0,̸=iλj.
λj − λi(S.19)

Applying this identity with λi = s(k − i), we ﬁnd

 ℓ−1
∏
k−j 

ℓ−1


∑

k−ℓ
−s(k−i)t  j=0
Qk (t) =
se

 ℓ−1

 ∏
i=0

i−j

(S.20)

j=0,̸=iWe can simplify this expression by noting thatℓ−1∏(k − j) =

j=0

k!
,
(k − ℓ)!

(S.21)and similarly thatℓ−1∏

(i − j) = i!(ℓ − 1 − i)!(−1)ℓ−1−i .

(S.22)

j=0,̸=iThis means we haveQk−ℓk (t) =ℓ−1∑

−s(k−i)t

sℓe

(
(−1)

ℓ−i−1

i=0

ℓ−1
i

)(

)
k
.
k−ℓ

(S.23)

We can evaluate this sum by recognizing the binomial expansion formula( )n,(1 + x) =xi
i=0(S.24)( ))ℓ−1
k −skt ( st
.
Qk−ℓ
(t)
=
sℓ
e
e −1
k
ℓ

(S.25)A. M. Walczak et al.

7 SI

n

n
∑

i

where we identify x = −est . We ﬁndMore generally, we have( ))a−b−1a −sat ( st
e
e −1
.
b

Qba (t) = s(a − b)

8 SI

A. M. Walczak et al.

(S.26)

FILE S3: GENERAL COALESCENCE PROBABILITIES IN THE NON-CONDITIONALAPPROXIMATION

The probability of coalescence for two individuals originally in two diﬀerent classes k and k ′ , as deﬁned in
Eq. (48) can be rewritten as′′

Pck,k →k −ℓ =1[I1 + I2 ] ,1 + 2N hk−ℓ s(k − ℓ)

(S.27)

where we have deﬁned
∫ ∞
I1

=
0

∫ ∞
I2

=
0

−s(k−ℓ)t1
Qk−ℓk′ (t1 )e−s(k−ℓ)t2
Qk−ℓk (t2 )e∫ t10s(k−ℓ)t2Qk−ℓ
dt2 dt1
k (t2 )e(S.28)s(k−ℓ)t1
Qk−ℓ
dt1 dt2 .
k′ (t1 )e(S.29)∫ t2
0

Note that both I1 and I2 involve integrals of the form
∫ t
Ia =

′

Qba (t′ )esbt dt′ .

(S.30)

0Plugging in the results for the non-conditional distributions of mutant timings, Eq. (S.26), and making useof the binomial expansion formula for (1 + x)n noted in File S2, we ﬁnd this integral becomesIa( )∫ t( ′)a−b−1′aes(b−a)t est − 1
dt′
b
0
(
)∫ t
( ) a−b−1
′
a−b−1
a ∑
(−1)a−b−1+is(a − b)es(b−a+i)t dt′
i
b
0
i=0
( )
(
)
a−b−1
)
∑ (−1)i a − b (a(a − b)(−1)a−bes(b−a+i)t − 1
b
a−b
i
i=0
( )
(
)
a−b
)
∑
a
a − b ( s(b−a+i)t(−1)a−b
(−1)i
e
−1
b
i
i=0
(
)
( )
a−b
∑
( st )i a − ba(−1)a−b es(b−a)t
−e
i
b
i=0
( )
)a−b
a s(b−a)t ( st
e
e −1
.
b

= s(a − b)

(S.31)=(S.32)====(S.33)(S.34)(S.35)(S.36)

We now substitute this result for Ia into our expressions for I1 and I2 . We note that both have terms of
the form

∫ ∞
Ib =
0

( )
)c−b
c −sct ( stdt.
Qba (t)ee −1b(S.37)A. M. Walczak et al.

9 SI

Using similar manipulations to those above, we ﬁndIb( )( ) ∫ ∞()a+c−2b−1a ce−s(a+c)t est − 1
dtbb 0( )( )∫ ∞a+c−2b−1∑ (a + c − 2b − 1)a c= s(a − b)(−1)a+c−1(−1)ie−s(a+c−i)t dt
b
b
i
0
i=0
(
)
( )( )
a+c−2b−1∑
a + c − 2b − 11a c.
= (a − b)(−1)i(−1)a+c−1
i
a+c−i
b
b
i=0
= (a − b)(S.38)(S.39)(S.40)

Using the partial fraction decomposition
( )
n
∑
1
n
i
(n+x) =,(−1)i−1i x+in
i=1
we ﬁnd( )( )

Ib =(S.41)( )( )

a c
a c
a−b
a−b
a+c
2b
a+c−2b b b (−1)
a+c−2b b b (−1)()=.
( −2b−1 )a+ca+c−2ba+c−2b

(S.42)′′

We can now use this result for Ib to determine I1 and I2 , and hence compute Pck,k →k −ℓ . We ﬁnd′′
Pck,k →k −ℓ =1
1 + 2N hk−ℓ s(k − ℓ)( k′ )( k )(k−ℓk+kk−ℓ′ ) .

(S.43)

2ℓ+k′ −kAs we noted in the main text, this is just′Pck,k →k−ℓ =′1
Ak,k ,
1 + 2N hk−ℓ s(k − ℓ) ℓ′

(S.44)

with Ak,kas deﬁned in Eq. (16). Note that when k = k ′ , this result simpliﬁes to Pck,k→k−ℓ as deﬁned inℓthe main text, as expected.

10 SI

A. M. Walczak et al.

FILE S4: COMPUTING SUMS OF ANCESTRAL PATHS′

In this Supplementary File, we describe the calculation of ϕkk (ℓ) using the sum of ancestral paths approach.
Calculation of ϕkk (3): We begin by considering a simpler speciﬁc case, where k = k ′ and ℓ = 3. There()
are a total of 63 = 20 possible ancestral paths by which two individuals sampled from class k can coalesce
in class k − 3. These can be separated into four types, according to whether the two ancestral lineages wereever together in classes k − 1 or k − 2. We can list all paths of each type, using the notation that A is amutation event in the ﬁrst lineage, and B is a mutation event in the second lineage. We have
ABABAB ABABBA  ABBAAB 






 ABBABA  BAABAB  BAABBA 






 BABAAB BABABA{z}|(21)(21)(21)=8 ways


ABAABB ABBBAA  BAAABB BABBAA|{z}(21)((42)−(21)(21))=4 ways
AABBAB AABBBA  BBAAAB BBAABA|{z}(21)((42)−(21)(21))=4 ways
AAABBB AABABB .
 BBBAAA 
BBABAA|{z}(63)−others=4ways

The probabilities of all paths of a particular type are identical. We can calculate the probability of eachof the four types of paths using the same logic as outlined in the main text. We ﬁnd
P (AAABBBc)=

P (AABBABc)=

P (ABAABBc)=

P (ABABABc)=()
k(k − 1)(k − 2)
1 − Ixk ,8(2k − 1)(2k − 3)(2k − 5)()()
k(k − 1)(k − 2)
Ixk−3
1 − Ixk 1 − Ixk−1 ,8(2k − 1)(2k − 3)(2k − 5)()()
k(k − 1)(k − 2)
Ixk−3
1 − Ixk 1 − Ixk−2 ,8(2k − 1)(2k − 3)(2k − 5)
(
)(
)(
)
k(k − 1)(k − 2)
Ixk−3
1 − Ixk 1 − Ixk−1 1 − Ixk−2 .
8(2k − 1)(2k − 3)(2k − 5)

Ixk−3

Summing over all the possible paths, we ﬁnd( k )( k ) [(2)(4)ϕkk (3) = Ik−3 k−3(2k)k−3
6(2)(4)

1 − (6) Ik−1 − (6) Ik−2 +
1

2

3

1

2

3

(2)(2)(2)](61) 1
3

Ik−1 Ik−2 .

1(S.45)(S.46)(S.47)(S.48)(S.49)

We now pause to consider the form of the probabilities of each type of ancestral path. These probabilities
diﬀer only by factors of (1 − Ixk−i ). One such factor arises each time the two ancestral lineages are together
in class k − i. In other words, we can rewrite the probability of each path as the probability of an undistorted
path (deﬁned to be a path in which the contributions due to the possibility of coalescence in previous classes
A. M. Walczak et al.

11 SI

are neglected), times a correction for each class in which the two lineages are together:
P (AAABBBc)=

P (AABBABc)=

P (ABAABBc)=

P (ABABABc)=()
P (Undistorted Path) 1 − Ixk()()
P (Undistorted Path) 1 − Ixk 1 − Ixk−1()()
P (Undistorted Path) 1 − Ixk 1 − Ixk−2
(
)(
)(
)
P (Undistorted Path) 1 − Ixk 1 − Ixk−1 1 − Ixk−2 .

(S.50)(S.51)(S.52)(S.53)

By deﬁnition, the “undistorted path” probability is the probability neglecting the contributions due to thepossibility of coalescence in previous steps, and is therefore the same for all paths. We have
P (Undistorted Path) ==k(k − 1)(k − 2)k(k − 1)(k − 2)I k−ℓ2k(2k − 1)(2k − 2)(2k − 3)(2k − 4)(2k − 5) x

(S.54)k!
k!
(k−3)! (k−3)! k−ℓIx .
2k!
(2k−6)!

(S.55)

Using these results, we can write ϕkk (3) as
ϕkk (3) =[[# of Paths] P (Undistorted Path) Fk (1 − Ixk ) + Fk,k−1 (1 − Ixk )(1 − Ixk−1 )]+Fk,k−2 (1 − Ixk )(1 − Ixk−2 ) + Fk,k−1,k−2 (1 − Ixk )(1 − Ixk−1 )(1 − Ixk−2 ) ,(S.56)

where we have deﬁned F{a} to be the fraction of paths that are together in the set of classes {a} (and are
not together in any other class).
Calculation of ϕkk′ (ℓ): We now use this approach to calculate the coalescence probability in the generalcase. The probability of any particular ancestral path from k and k ′ to k − ℓ is the product of the individual
probabilities of each mutational step that makes up this path. Each such individual probability consistsof three parts: a numerator, which depends only on the current class of the lineage that mutates, divided
by a denominator, which depends only on the sum of the current set of classes for both lineages, times a
correction factor of (1 − Ixk−i ) if the two lineages are in the same class at that step.
Although in each ancestral path the mutations will occur in a diﬀerent order, all paths will ultimately
consist of the same set of mutations (k ′ → k ′ − 1 → . . . → k − ℓ and k → k − 1 → . . . → k − ℓ). Therefore,
regardless of the path taken, the product of the numerators from each step will be identical. Similarly, thesum of the current set of classes will begin at k ′ + k, and decrement by one each time a deleterious mutation
occurs, until both lineages are in the ﬁnal class (k ′ + k → k ′ + k − 1 → . . . → 2k − 2ℓ). Therefore, regardlessof the path taken, the product of the denominators from each step will also be identical. Therefore, thepaths will diﬀer only by the correction factor (1 − Ixk−i ) for each class in which the two ancestral lineages
are together. This means that, analogous to the case of ϕkk (3) we described above, the probability of each
12 SI

A. M. Walczak et al.

path is the probability of an “undistorted path” times the appropriate correction factor. The probability of
the undistorted path is
P (Undistorted Path) =k ′ (k ′ − 1) . . . (k − ℓ + 1)k(k − 1) . . . (k − ℓ + 1) k−ℓ
Ix .
(k ′ + k)(k ′ + k − 1) . . . (2k − 2ℓ + 1)

(S.57)

We can now sum up all possible paths to obtain[ϕkk′ (ℓ)=[# of Paths] P (Undistorted Path) F∅ +ℓ∑

Fk−i (1 − Ixk−i )

i=0

+

ℓ−1 ∑
ℓ
∑

Fk−i,k−j (1 − Ixk−i )(1 − Ixk−j )(S.58)

i=0 j>i

+

ℓ−2 ∑
ℓ−1 ∑
ℓ
∑


Fk−i,k−j,k−m (1 − Ixk−i )(1 − Ixk−j )(1 − Ixk−m ) + . . . ,i=0 j>i m>jwhere as before F{a} is the fraction of paths that are together in the set of classes {a} (and are not together
in any other class). Note that there are a total of ℓ + 1 terms in this equation, representing the possibility
that the two lineages can be together in anywhere from 0 to ℓ of the classes. We can rearrange these terms
to write[ϕkk′ (ℓ)=[# of Paths] P (Undistorted Path) 1 −ℓ∑

Gk−i Ixk−i

i=0

+

ℓ−1 ∑
ℓ
∑

Gk−i,k−j Ixk−i Ixk−j(S.59)

i=0 j>i

−

ℓ−2 ∑
ℓ−1 ∑
ℓ
∑


Gk−i,k−j,k−m Ixk−i Ixk−j Ixk−m + . . . ,i=0 j>i m>jwhere we have deﬁned G{a} to be the fraction of paths that are together in at least the set of classes {a}.
We can evaluate each of these factors of G. For example, the fraction of paths that are together in classk − i equals the number of ways for the two lineages to descend from classes k ′ and k to be together in class( ′)k − i, k −k+2i, times the number of ways for the two lineages to descend from class k − i to be together ini(2i−2ℓ)class k − ℓ, i−ℓ , divided by the total number of ways for the two lineages to descend from classes k ′ and( ′)k to be together in k − ℓ, k −k+2ℓ. Using this logic, we ﬁndℓϕkk′ (ℓ)=[# of Paths] P (Undistorted Path)(S.60)


)
(
)(
)(
(
)(
)
′
′
ℓ−1 k −k+2i 2ℓ−2iℓ−2 ∑ℓ−1 k −k+2i 2j−2i 2ℓ−2j
∑
∑
i
j−i
ℓ−j
× 1 −
Ixk−i +
Ixk−i Ixk−j . . . .
(ik′ −k+2ℓℓ−i)(k′ −k+2ℓ)ℓ

i=0

The total number of paths is(k′ −k+2ℓ)ℓ

i=0 j>iℓ, so we ﬁnally ﬁnd that the full probability of coalescence in class
A. M. Walczak et al.

13 SI

k − ℓ is(′ϕkk (ℓ)

=

)(

)[

k′
k
Ixk−ℓ (k−ℓk′ +kk−ℓ)
k′ −k+2ℓ

1−

ℓ−1
∑

(k′ −k+2i)(2ℓ−2i)i=0

Ixk−i +
(ik′ −k+2ℓℓ−i
)
ℓ


( ′
)(
)(
)
ℓ−2 ∑
ℓ−1 k −k+2i 2j−2i 2ℓ−2j
∑
i
j−i
ℓ−j
Ixk−i Ixk−j − . . . .
(k′ −k+2ℓ)
i=0 j>i(S.61)ℓ′

This is Eq. (56) from the main text. Note that it equals our non-conditional result for Pck,k →ℓ times a
correction factor. There are a total of ℓ + 1 terms in this correction factor. This full correction factor can bearbitrarily complex for large ℓ, so we do not write out a general form here. However, it is straightforward to
calculate for any values of k, k ′ , and ℓ; a Mathematica script to do so is available on request.

14 SI

A. M. Walczak et al.

FILE S5: THE CORRESPONDENCE BETWEEN STEPTIMES AND REAL TIMESIn this Supplementary File, we calculate the correspondence between steptimes and the actual times measured
in generations. Our goal is to calculate the probability distribution of real coalescence times, ψ(t|k, k ′ , ℓ),given that individuals were initially in classes k and k ′ and coalesced in class k − ℓ.
To begin, we neglect the coalescence time within class k − ℓ, and consider the distribution of the time atwhich an ancestor of one of the two sampled individuals ﬁrst mutated from class k − ℓ to class k − ℓ + 1. We
refer to this as ψ1 (t|k, k ′ , ℓ). We ﬁrst calculate the joint distribution of the times at which both ancestorsk−ℓ
k−ℓ
mutated out of the class, Rk,k
′ (t1 , t2 ). Conditional on coalescence in class k − ℓ, Rk,k ′ (t1 , t2 ), is given by theprobability of t1 and t2 and coalescence divided by the total probability of coalescence. That is,R(t1 , t2 ) =P (coal|t1 , t2 )P (t1 , t2 ).
P (coal)(S.62)

Substituting in the relevant expressions from the main text, this gives
k−ℓ
Rk,k′ (t1 , t2 ) =

1
′
Ak,k
ℓ

−s(k−ℓ)|t1 −t2 |Qk−ℓ.
k,k′ (t1 , t2 )e(S.63)The time at which the ﬁrst ancestor mutated out of class k − ℓ is the longer of the two times t1 and t2 ,ψ(t|k, k ′ , ℓ) =

[∫ t
0

∫ t
k−ℓ
Rk,k
′ (t1 , t)dt1 +

0

]
k−ℓ
Rk,k
(t,
t
)dt
′
2
2 .

(S.64)

k−ℓ
Substituting in our expression for Rk,k
′ (t1 , t2 ) and carrying out the integrals as in File S3, we ﬁnd′ψ1 (t|k, k ′ , ℓ) = sπd e−s(k +k)t (est − 1)πd −1

( ′
)
k +k
,
πd

(S.65)

where we have used πd = k ′ − k + 2ℓ.
We can alternatively calculate ψ1 (t|k, k ′ , ℓ) using our sum of ancestral paths approach. As before, we
imagine two individuals sampled from classes k and k ′ and condition on them coalescing in class k − ℓ.
Consider a case where k ̸= k ′ . Then the ﬁrst event in the history of these two individuals must be adeleterious mutation. Since these mutations happen at rate sk and sk ′ in each lineage, the distribution oftimes since this mutation occurred in one of the two ancestral lineages is′

P (t) = s(k + k ′ )e−s(k+k )t .

(S.66)′k′With probability k+k′ , this mutation is in the lineage sampled from class k , in which case the two lineages

are now in classes k and k ′ − 1. Alternatively, the mutaion occurred in the lineage sampled from k and thelineages are in classes k − 1 and k ′ .
We can now consider the time to the next event backwards in time. If the two lineages are in the same
A. M. Walczak et al.

15 SI

class (but not yet in class k − ℓ), the distribution of times to the next deleterious mutation event is somewhat
shorter, because we are conditioning on coalescence not occuring. However, provided that 2sk1 ≫ N1hk (thecondition we are already making elsewhere), this shortening of the time will be a small correction andneglecting it is a good approximation.
Making this approximation, the rate at which the next deleterious mutation event occurs when the twolineages are in classes k1 and k2 is just s(k1 + k2 ). Regardless of the order in which these mutations happen
between the two lineages, this sum is simply decreased by s at each step. This will continue until the bothancestral lineages are in class k − ℓ. Therefore, the distribution of times until the original mutation out ofclass k − ℓ is given by:′′

ψ1 (t|k ′ , k, ℓ) = s(k ′ + k)e−s(k +k)t ⋆ s(k ′ + k − 1)e−s(k +k−1)t ⋆ . . . ⋆ s(2k − 2ℓ + 1)e−s(2k−2ℓ+1)t .

(S.67)

This can be written as
ψ1 (t|k ′ , k, ℓ) = λ0 e−λ0 t ⋆ λ1 e−λ1 t ⋆ . . . ⋆ λk′ −k+2ℓ−1 e−λk′ −k+2ℓ−1 t ,(S.68)

where we have deﬁned:
λi = s(k ′ + k − i).

(S.69)

We can compute this convolution as in File S2 (compare to Eq. (S.17) for Q2k−2ℓk+k′ (t)). We ﬁnd′ψ1 (t|k, k , ℓ) = sπd e−s(k′ +k)t(e − 1)
st

πd −1

( ′
)
k +k
,
πd

(S.70)

identical to the result of our lineage structure calculation above.
Distribution of Coalescence Times: To calculate the correspondence between steptimes and realtimes, we now need to add the time it takes two individuals two coalesce in class k − ℓ, which we refer to as
ψ2 (t|k, k ′ , ℓ), to the time it took them both to get to that class, ψ1 (t|k, k ′ , k − ℓ). The rate of coalescenceonce in class k − ℓ is N h1k−ℓ , so we have
ψ2 (t|k ′ , k, ℓ) = (2s(k − ℓ) + 1/N hk−ℓ ) e−[2s(k−ℓ)+1/N hk−l ]t .

(S.71)

Putting this together, the full distribution of times since coalescence isψ(t|k ′ , k, ℓ) = ψ1 (t|k ′ , k, ℓ) ⋆ ψ2 (t|k ′ , k, ℓ).

(S.72)

Carrying out this convolution (and expanding the binomial factor (est − 1)πd −1 in ψ1 ), we ﬁnd′ψ(t|k , k, ℓ) =

π∑
d −1
i=0

16 SI

sπd (−1)πd −i−1()())
πd − 1 k ′ + k
B ( −sBte− e−sAt ,iπdA−B

A. M. Walczak et al.

(S.73)

where we have deﬁned A ≡ k ′ + k − i and B ≡ 2 (k − ℓ) + N sh1k−ℓ .

A. M. Walczak et al.

17 SI

FILE S6: AN ALTERNATIVE APPROACH TO NEUTRAL DIVERSITY
Instead of calculating the distribution of neutral heterozygosity by ﬁrst computing the distribution of realtimes, we could alternatively incorporate neutral mutations directly into the sum of ancestral paths framework. This completely bypasses the correspondence with real coalescence times. To do this, we characterize
ancestral paths not only by the ordering of deleterious mutation and coalescence events, but also by theordering of neutral mutations. This means that if we sample two individuals A and B, there are ﬁve typesof events that can happen in their ancestral paths: a deleterious mutation (DM) in A or in B, a neutral
mutation (NM) in either A or in B, and or a coalescence (C) event (if A and B are currently in the sameclass).
We now imagine that we sample two individuals from classes k and k ′ , and that they coalesce in classk − ℓ. Our goal is to calculate the probability distribution of πn given k, k ′ , and ℓ, ρ(πn |k, k ′ , ℓ). We willﬁnd it helpful to divide the ﬁve types of events that can occur into two classes: neutral mutations on the one
hand, and deleterious mutations or coalescence (which we call “steps”) on the other. We begin by computing
the probability that a given number of NMs occur before the next DM or C events (i.e. the number of neutral
mutations that occur at this “step”). We have
(
′

′

′

P (a NMs, then DM in k or k |k , k) =

)a

2Un
s

k ′ + k + 2Us n

k + k′,
k ′ + k + 2Us n(S.74)

where we have made our usual assumption that N hk sk ≫ 1, allowing us to neglect the rates of coalescence
events (when k = k ′ ) in writing this expression.
This probability only depends on the sum of the current classes the individulas are in. At each subsequent
step, regardless of the path taken, this sum of the classes will decrease by one. Therefore, the probability
that ai neutral mutations occur at step i is independent of the path taken. This observation allows us tocalculate the probability that a given total number of neutral mutations have occurred since coalescence.
We ﬁrst calculate the probability that a given number of neutral mutations have occurred since the ﬁrst
deleterious mutation out of the k − ℓ class. We will add in the additional neutral mutations once in the k − ℓclass at the end.
In order for πn neutral mutations to have occurred since the ﬁrst deleterious mutation out of class k − ℓ,we require that a0 mutations occurred at the ﬁrst step, a1 mutations occurred at the second step, and soon, such that a0 + a1 + . . . + ak′ −k+2ℓ−1 = πn . This gives′ρ(πn = X|k , k, ℓ) =18 SI

(k′ +k)!
(2k−2ℓ)!

∑ (( 2Usn +k′ +k)!
a|=X
( 2Usn +2k−2ℓ)! |⃗2Un /s2Un /s + k + k ′

)a0

(
...

A. M. Walczak et al.

2Un /s2Un /s + 2k − 2l + 1)ak′ −k+2l−1. (S.75)

We can deﬁne x ≡ 2Un /s + k + k ′ , recognize πd = k ′ − k + 2ℓ, and relabel the ai as
a0 → X − b0 ,a1 → b0 − b1 ,aπd −2 → bπd −3 − bπd −2 ,...

aπd −1 → bπd −2 .

(S.76)

This gives(k′ +k)′ρ(πn = X|k , k, ℓ) =(πd

( 2Un +k′ +k)
s

πd

2Un
s

)b
b0 (
∑
x−1 1
x−2

b1 =0

)X ( )X ∑
)b0
X (
1
x
x
x−1

(S.77)b0 =0

bπd −3∑...

bπd −2 =0(x − πd + 2
x − πd + 1)bπd −2.

To simnplify this expression, it is helpful to deﬁne a function f such that:f (A, B)≡( )X ∑
)b0
X (
x
1
x
x−1
b0 =0

)b
b0 (
∑
x−1 1
b1 =0

x−2

(S.78))b bA−1 ()b(X∑
x−A+1 0 ∑ x−A A
...
x−A
x−B
bA−1 =0bA =0In other words, f (A, B) is a set of A nested sums, each of the same form, except for the ﬁnal sum, which
can have a diﬀerent denominator. Using this deﬁnition, we have(k′ +k)′P (πn = X|k , k, ℓ) = ( 2Unsπd+k′ +k
πd

(
)

2Un
s

)X
f (πd − 2, πd − 1) .

(S.79)The virtue of this deﬁnition is that this sum can be solved recursively. We havebA−1 (

∑

bA =0

x−A
x−B

)bA

x−B
x−A
=
−
A−B
A−B

(

x−A
x−B

)bA−1
.

(S.80)x−A
x−B
f (A − 1, B) −f (A − 1, A) .
B−A
B−A

(S.81)

Therefore we havef (A, B) =Repeatedly inserting this result yields:f (A, A + 1) →f (A, A + 1) →..
.
f (A, A + 1) →()(x − A)(x − A − 1) f (A − 1, A + 1) f (A − 1, A)
−
1
x−A−1
x−A
[
]
(x − A + 1)(x − A)(x − A − 1) f (A − 2, A + 1) 2f (A − 2, A) f (A − 2, A − 1)
−
+
2
x−A−1
x−A
x−A+1
( )
) m
m
x − A − 1 + m ∑ (−1)i+m
f (A − m, A + 1 − i) .
(m + 1)
x
−
A
−
1
+
i
i
m+1
i=0

(S.82)A. M. Walczak et al.

19 SI

(

Note that f (−1, B) = 1/B X , since there are no more sums to compute. Thus, for m = A + 1 we have
(

) A+1
(
)
∑
x
(−1)i+A+1
A+1
f (A, A + 1) = (A + 2).
A + 2 i=0 (x − A − 1 + i)X+1
i(S.83)

Relabeling the sum and taking A = πd − 2, we have(f (πd − 2, πd − 1) = πd

x
πd

(
)
πd − 1(−1)i.
i(x − i)X+1
i=0

) π∑
d −1

(S.84)

We can now substitute these results into our expression for πn , to ﬁndρ1 (πn = X|k ′ , k, ℓ) = πd

)X π∑
(
( ′
)
)(
d −1
(−1)i
πd − 1k +k
2Uns(2Un /s + k + k ′ − i)X+1
i
πd
i=0

(S.85)

Note, however, that this is only the distribution of neutral mutations since the ﬁrst deleterious mutation out
of class k − l. It is also possible for neutral mutations to occur prior to the coalescence event. Adding in thisfactor, we ﬁnd) πd −1()( ′k +k ∑
i πd − 1
πd
(−1)
πd
i
i=0

′

ρ(πn = X|k , k, ℓ) =

×

πn
∑
X=0

X

(2Un /s)(2Un /s + k + k ′ − i)X+1(S.86)(

2Nk−l Un
1 + 2Nk−l Un + 2Nk−l s(k − l))πn −X.

Rearranging this expression gives′ρ(πn |k , k, ℓ) =

π∑
d −1

πd −i−1πd (−1)

i=0

(
)
(
)(
)
( 2Us n )πn( 2Us n )πnπd − 1 k ′ + k
B, (S.87)− 2Un
A − B ( 2Us n + B)πn +1i
πd( s + A)πn +1

where we have deﬁnd
A = k ′ + k − i,

B = 2 (k − ℓ) +identical to our earlier result.

20 SI

A. M. Walczak et al.

1,
N shk−l(S.88)