Theoretical Population Biology 81 (2012) 144–157

Contents lists available at SciVerse ScienceDirect

Theoretical Population Biology
journal homepage: www.elsevier.com/locate/tpb

The structure of allelic diversity in the presence of purifying selection
Michael M. Desai a,b,c,1 , Lauren E. Nicolaisen a,b,c,1 , Aleksandra M. Walczak d , Joshua B. Plotkin e,∗a

Department of Organismic and Evolutionary Biology, Harvard University, United Statesb

Department of Physics, Harvard University, United Statesc

FAS Center for Systems Biology, Harvard University, United Statesd

CNRS-Laboratoire de Physique Théorique de l’École Normale Supérieure, 24 rue Lhomond, 75005 Paris, FranceeDepartment of Biology, University of Pennsylvania, United Statesarticleinfo

Article history:
Received 13 September 2011
Available online 16 December 2011
Keywords:
Allelic diversity
Purifying selection
Ewens sampling formula
Linkageabstract
In the absence of selection, the structure of equilibrium allelic diversity is described by the elegant
sampling formula of Ewens. This formula has helped to shape our expectations of empirical patternsof molecular variation. Along with coalescent theory, it provides statistical techniques for rejecting thenull model of neutrality. However, we still do not fully understand the statistics of the allelic diversity
expected in the presence of natural selection. Earlier work has described the effects of strongly deleterious
mutations linked to many neutral sites, and allelic variation in models where offspring fitness is unrelated
to parental fitness, but it has proven difficult to understand allelic diversity in the presence of purifying
selection at many linked sites. Here, we study the population genetics of infinitely many perfectly linked
sites, some neutral and some deleterious. Our approach is based on studying the lineage structure withineach class of individuals of similar fitness in the deleterious mutation-selection balance. Consistent withprevious observations, we find that for moderate and weak selection pressures, the patterns of allelicdiversity cannot be described by a neutral model for any choice of the effective population site. We
compute precisely how purifying selection at many linked sites distorts the patterns of allelic diversity,
by developing expressions for the likelihood of any configuration of allelic types in a sample analogous to
the Ewens sampling formula.
© 2011 Elsevier Inc. All rights reserved.

1. Introduction
In any evolving population, new clonal lineages are constantly
being created and destroyed. The balance between the creationof lineages by new mutations and their destruction by natural
selection and genetic drift determines the statistics of theclonal structure of the population. In the absence of natural
selection, Ewens (1972) computed an elegant sampling formula
describing the clonal structure of a neutral population, andexplained how the allelic (i.e. lineage) configuration in a sample ofindividuals from the population provides a window into this clonal
structure.
Natural selection distorts the clonal structure of a population
away from this neutral expectation. Of particular interest is purifying (negative) selection against many linked deleterious mutations(‘‘background selection’’). Recent evidence has suggested that this

∗ Correspondence to: Department of Biology, University of Pennsylvania, 219
Lynch Labs, 433 S. University Avenue, Philadelphia, PA 19104, United States.
E-mail address: jplotkin@sas.upenn.edu (J.B. Plotkin).
1 These authors contributed equally to this manuscript.
0040-5809/$ – see front matter © 2011 Elsevier Inc. All rights reserved.
doi:10.1016/j.tpb.2011.12.002

may be generally important in a wide range of populations (seeHahn (2008) for a recent review). In this paper, we explore how thistype of selection alters the clonal (i.e. allelic) structure of a population. Our analysis leads to a generalization of the Ewens sampling
formula to situations involving background selection.
Over the past few decades, numerous authors have studied
allelic diversity in infinite-alleles frameworks that incorporate
selection. Li (1977) and Watterson (1978) introduced models inwhich alleles may have a few different selective effects. Li (1978)and others (Li, 1979; Ewens and Li, 1980; Griffiths, 1983) analyzed
the structure of allelic diversity in these models. More recent
work has analyzed a very general model of selection introducedby Ethier and Kurtz (1987), which allows for diverse types of
selection pressures (Ethier and Kurtz, 1994; Joyce and Tavare,
1995; Grote and Speed, 2002; Joyce, 1995). This work has helped
us understand the general effects of selection in distorting thefrequency spectrum of sampled alleles. However, the models theseauthors have analyzed cannot be directly connected to a concrete
description of mutations and selection occurring at specific sites.
Rather, they assume that each new mutation creates a new allele
whose fitness is completely independent of the fitness of its
parent. In other words, there is no sense of relatedness among

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

alleles, or of a correlation in fitness between closely related alleles.
Etheridge and Griffiths (2009) and Etheridge et al. (2010) havemore recently derived a coalescent dual of the Moran process withan arbitrary number of types, mutation rates between types, andgenic selection coefficients, but it is not clear how this corresponds
to selection acting on some fraction of an infinite number of
specific sites.
In this paper we take a different approach, based on the specificmodel of linked sites described by Charlesworth et al. (1993) andHudson and Kaplan (1994). That is, we imagine that each individual
has a genome comprised of many neutral and many negatively
selected sites. The fitness of each individual is determined by the
number of mutations it carries at the negatively selected sites. We
make the infinite-sites assumption that no two mutations at thesame site ever segregate simultaneously. This is also an infinitealleles model, but it is based on a specific model of mutations at
individual sites, and the fitness of each new allele depends on thefitness of its parent.
Earlier studies have investigated the effects of purifying
selection in models identical or closely related to the one we
consider here. Charlesworth et al. (1993) introduced a model
essentially identical to the one we analyze here, and Kaplan
et al. (1988) and Hudson and Kaplan (1994) developed a simple
algorithm which can be used to recursively compute how
purifying selection alters the structure of genealogies. Hudson
and Kaplan (1995) and Gordo et al. (2002) further developed thisidea, resulting in a simple computational method for sampling
genealogical relationships in the presence of background selection.
Related simulation and analytical work has further characterized
the structure of genealogies and the statistics of genetic diversity
at the level of individual sites in this or closely related models(McVean and Charlesworth, 2000; Seger et al., 2010; Charlesworth
et al., 1993; Comeron and Kreitman, 2002; Comeron et al., 2008;
Barton and Etheridge, 2004). However, this earlier work does notprovide an analytic description of lineage structure, or sampling
formulas for allelic diversity in the presence of purifying selection
on many linked sites.
In this paper, we explicitly analyze the lineage structure, andwe derive a selected version of the Ewens sampling formula. We
begin by noting that the balance between mutations at deleterious
sites and selection against them leads to a steady state mutationselection balance (Haigh, 1978). Our approach is to study the
structure of lineages within this steady state, using the Poisson
Random Field (PRF) method developed by Sawyer and Hartl(1992). We show that this lineage structure can alternatively
be derived using a retrospective approach, by considering theprobabilities of mutation and coalescence events in the ancestryof each individual; these probabilities are calculated by Hudson
and Kaplan (1994) and Gordo et al. (2002) (and implicitly in a
related context by Barton and Etheridge (2004)). Our descriptionof lineage structure is thus precisely consistent with the analysis
of genealogical structures in this earlier work. Finally, we use ourdescription of lineage structure to calculate sampling formulas
for allelic diversity, and compare our predictions to the results ofMonte Carlo simulations.
Provided that selection is strong and deleterious mutation
rates are sufficiently small, our results show that the effect ofbackground selection on allelic diversity is to reduce the effective
population size without otherwise distorting the lineage structure.
Our results are thus consistent with the effective population size
approximation to background selection proposed by Charlesworth
et al. (1993). For weaker selection, however, or higher mutation
rates, the effective population size approximation breaks down,
and the effects of background selection become more complex.
We show that in this case the allelic diversity cannot be describedby neutral theory with some appropriately chosen effective

145

population size. This is consistent with earlier observations thatbackground selection leads to distortions in the structure ofgenealogies (McVean and Charlesworth, 2000; Seger et al., 2010;
O’Fallon et al., 2010; Comeron and Kreitman, 2002; Comeron et al.,2008; Barton and Etheridge, 2004; Gordo et al., 2002; Hermisson
et al., 2002; Williamson and Orive, 2002). Our analysis here allows
us to compute precisely how these distortions due to purifying
selection at many linked sites alter patterns of allelic diversity,and hence provides an analytical framework for exploring where
statistical power may lie to distinguish purifying selection fromneutrality.
Our approach relies on the assumption that we can describe
the distribution of fitnesses within the population with thesteady state mutation-selection balance. In particular, we neglect
fluctuations within this balance. We note that the PRF andretrospective approaches depend somewhat differently on thiskey approximation, which offers some insight into the roleof fluctuations in our model. We analyze the validity of thisapproximation in more detail below, and describe a correction
for some aspects of the effects of fluctuations in the PRF
formalism, which allows us to make a precise correspondencewith the retrospective approach. Related to this approximation,
we also neglect the effects of Muller’s ratchet. We discuss thisapproximation in detail in the Discussion. We further test thevalidity of our analysis via Monte Carlo simulations; we find thatthese approximations are reasonable across a broad parameter
regime spanning weak and strong selective pressures.
Our analysis in this paper is limited to allelic diversity, andit does not address the degree of relatedness among sampled
alleles. In other words, our analysis only tells us the probability
that individuals are genetically identical, not the distribution of thenumber of specific sites at which individuals may differ. Our results
are thus not directly comparable to the work described above,
which makes predictions about expected diversity at the level ofindividual sites. However, while our allele-based results provide
an incomplete picture of genetic diversity within the population,
they do provide a useful perspective on how purifying selection
distorts patterns of molecular evolution. Most importantly, we areable to make precise analytical predictions about how purifying
selection distorts allelic diversity, in ways that cannot be describedby a single reduced effective population size.
2. ModelWe imagine a finite haploid population of constant size N.
Each haploid genome has a large number of sites, which begin insome ancestral state and mutate at a constant rate. Each mutationis either neutral or confers some fitness disadvantage s (whereby convention s > 0). We assume an infinite-sites framework,so there is negligible probability that two mutations segregate
simultaneously at the same site.
We assume that there is no epistasis for fitness, and that eachdeleterious mutation carries fitness cost s, so that the fitness of anindividual with k deleterious mutations is wk = (1 − s)k . Since weassume that s ≪ 1, we will often approximate wk by 1 − sk. Laterwe comment briefly on extensions to our method to consider thecase when the selection coefficient of a deleterious mutations isdrawn from some fixed distribution.
The population dynamics are assumed to follow the diffusion
limit of the standard Wright–Fisher model. That is, 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 define θd /2 ≡ NUd , theper-genome scaled deleterious mutation rate. Similarly, neutral
mutations occur at a rate Un per individual per generation, andwe analogously define θn /2 ≡ NUn . We assume that each newly

146

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

arising mutation occurs at a site at which there are no other
segregating polymorphisms in the population (the infinite-sites
assumption). Since in this paper we focus only on allelic diversity,
this infinite-sites approximation simply means that each new
mutation creates a unique allele. Throughout the analysis we
assume that Muller’s ratchet can be neglected; we discuss thevalidity of this approximation in the Discussion.
We study the case of perfect linkage. In other words, we
imagine that all the sites we are considering are in an asexual
genome or within a short enough distance in a sexual genome thatrecombination can be entirely neglected. Although our model isdefined for haploids, this assumption means 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 fitness cost as being heterozygous).
We believe that this is the simplest possible model based on aconcrete picture of mutations at individual sites that can describe
the effects of a large number of linked negatively selected sites
on patterns of genetic variation. It is essentially equivalent to themodel described by Charlesworth et al. (1993) and Hudson and
Kaplan (1994), which has formed the basis for much of the analysis
of background selection (Gordo et al., 2002; Seger et al., 2010).
3. Analysis
The balance between mutations and selection leads to asteady state distribution of fitnesses within the population; thisis the well-known ‘mutation-selection balance’. However, theindividuals of a given fitness are not all genetically homogeneous,
but rather comprise a number of different alleles. The number andfrequency distribution of these alleles depends on how quickly
new alleles are created by deleterious mutations from more-fit
individuals, and hence on the overall fitness distribution.
We begin by describing the relevant aspects of the mutationselection balance that leads to a steady state distribution of
fitnesses within the population. Our description of this steady state
fitness distribution is entirely deterministic. Of course, in a finite
population, there will be random fluctuations in the values of hk ,the fraction of the population harboring k deleterious mutations. In
the most extreme case, these fluctuations lead to Muller’s ratchet.
In our analysis below, we will neglect these fluctuations in hk ,assuming that these frequencies are always at their deterministic
steady state. Consistent with this approximation, we will also
neglect the effects of Muller’s ratchet. We will then return in a later
section to use our results to determine when these approximations
are valid.
If we assume for a moment that these approximations arereasonable, we can already guess the form of our result for theallelic diversity. New alleles are constantly being generated withinfitness class k due to deleterious mutations from class k − 1and neutral mutations from class k. Within class k, all alleles
drift neutrally with respect to each other. Therefore, conditional
on mutations and selection keeping the frequency of the class at
hk , the allelic diversity within this class will be the same as in aneutral population of size Nhk in which new alleles are created bymutations at the appropriate rate. Thus for example the probability
two individuals are of the same allelic type is the probability thatthey are both in the same class k times the appropriate neutral
result for the homozygosity within that class, summed over all
possible classes. Sampling formulas for larger samples can becalculated in the analogous way.
The remainder of our analysis in this paper is, essentially,
devoted to making this simple intuition precise and showing when
it is accurate. We start by summarizing earlier results for the steadystate mutation-selection balance hk , and then compute the allelicdiversity in detail, neglecting all fluctuations in hk . This allows us
to see precisely when this approximation is reasonable, and hence
prove when the simple intuition described above holds.

baFig. 1. Schematic of the allelic diversity in the mutation-selection balance. (a)U
Sketch of the mutation-selection balance in the case sd = 5. The steady state
distribution of fitness within the population is maintained by a balance between
mutations moving individuals toward lower fitness and selection favoring thoseclasses more fit than average at the expense of those less fit than average. (b) The
inset shows the processes maintaining a class of individuals with k deleterious
mutations. Deleterious mutations from class k − 1 found new lineages within classk at rate Nhk−1 Ud . Neutral mutations found new lineages in the class at a rate
Nhk Un . Selection favors or disfavors individuals from each lineage at a per capita
rate −(k − k̄)s, and deleterous mutations eliminate individuals from each lineage
at a per capita rate Ud + Un .

3.1. The steady state fitness distribution
In our model, all deleterious mutations have the same fitness
cost s, so we can characterize individuals by their Hamming class,k, relative to the wildtype (which by definition has k = 0). Thatis, individuals in class k have k deleterious mutations more thanthe most-fit individuals in the population. Here k refers only tothe number of deleterious mutations an individual has; individuals
with the same k can have different numbers of neutral mutations.
We normalize fitness such that by definition all individuals in classk = 0 have fitness 1. Individuals in class k then have fitness 1 − ks.
Imagine that at a given time a fraction hk (t ) of the population isin class k. This class is acquiring new individuals due to deleterious
mutations arising in class k − 1, and it is losing individuals due todeleterious mutations away to class k + 1. It also gains or loses
individuals at a rate −(k − k̄)s due to selection, where k̄ is themean k within the population, k̄ ≡khk . This is illustrated inFig. 1. Note that the term involving k̄ simply normalizes the effectof selection (selection favors a class if it is more fit than the average
individual, and vice versa). This means that on average hk (t ) will
evolve according to the equationdhk (t )dt= Ud hk−1 − Ud hk − (k − k̄)hk s.

(1)

Note this is a system of k equations for all the hk (t ). Of
course random genetic drift will also affect the hk (t ), and these
deterministic equations are only true on average. We return to thispoint below, but for now we neglect drift and focus on the steadystate distribution.
The steady state fitness distribution (the mutation-selection
balance) is given by the values of hk (t ) after a long time. We
can find this mutation-selection balance by setting the right handside of Eq. (1) equal to 0 for all values of k. This calculation was
originally carried out by Kimura and Maruyama (1966) and Haigh(1978); they found that the steady state, ĥk , is given by a Poisson
U
distribution with mean sd ,ĥk =e−Ud /s
k!



Ud
s

k

.

(2)

Note that this means that the average fitness in the population isU
1 − Ud , and that k̄ = sd .

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

3.2. Allelic diversity within a given fitness class
We now look more closely at individuals within a given fitnessclass, as illustrated in Fig. 1b. For the moment we neglect neutral
mutations; we consider their effects further below.
All lineages in class k originally arose from a deleterious
mutation to an individual in class k − 1. Each of these deleterious
mutations founds a new lineage within class k. Such lineages arefounded at a rate θk /2, where we defineθk = 2Nhk−1 Ud .

(3)

Note this is true whether or not the hk are at their steady-statevalues, though for the purposes of our analysis we will always
assume the steady state.
In our infinite-alleles approximation, each new lineage is anallele that is unique within the population. The fate of this lineage(allele) is then determined by the forces of random drift, selection,
and additional mutations. Additional mutations that occur within
this lineage go on to found new alleles. Thus from the point of viewof this particular lineage, additional mutations cause individuals to
be lost from the lineage. This means that individuals are removed
from a lineage in class k at a per capita ratesk ≡ −Ud − s(k − k̄).

(4)

(5)

Note that in steady state, when the fitness distribution hk takes
the mutation-selection balance form ĥk derived above, k̄ = Ud /sand the effective selection coefficient sk is negative for all fitnessclasses with k > 0. This makes intuitive sense: each fitness class(except k = 0) is constantly receiving new individuals due tomutations. Thus older individuals must on average die out, if thefitness class is to stay at a constant steady state size. The only
exception is the k = 0 class, for which sk = 0. This class drifts
effectively neutrally, with its actual selective advantage relative
to the mean exactly balanced by the loss of individuals due todeleterious mutations. For k = 1 we have s1 = −s, and in generalsk = −ks. On the other hand, θk /hk increases with k, reflecting
the fact that the stronger selection against the larger-k classes isbalanced by a larger influx of new deleterious mutations into theseclasses.
We can now incorporate the effect of neutral mutations. Eachneutral mutation within an individual in class k creates a newlineage in class k. Thus we may simply redefine the rate at which
new lineages are founded, givingθk ≡ 2Nhk−1 Ud + 2Nhk Un .

(6)

When the hk ’s are in steady state this definition simplifies to θk =2Nhk (sk + Un ). Each neutral mutation also causes an individual to
be lost from the lineage it was in before the mutation, so we also
redefine the effective selection coefficientsk ≡ −Ud − Un + s(k − k̄).

lineage is independent of all the others. This assumption is valid
provided that no lineage ever becomes a substantial fraction of
the overall population, which will be true whenever N |sk | ≫ 1(i.e. all lineages are selected against strongly enough). A sufficientcondition for this to hold in the bulk of the fitness distribution is
simply N (Un + Ud ) ≫ 1, and in fact our approximation will also
hold even in some circumstances when this condition breaks down
(we describe this further below).
3.3. Poisson Random Field description of a lineage structure
Using the independence assumption, we have reduced the
problem of describing a lineage within a given fitness class toexactly the situation addressed by the Poisson Random Field modelof Sawyer and Hartl (1992). Thus the frequency distribution of
lineages (alleles) in fitness class k is a Poisson Random Field (PRF)with parameters θk and γk (where as before γk ≡ Nsk ). That is, thenumber of distinct lineages in class k segregating at a frequency
between a and b in the entire population is Poisson distributed with
mean

 b

fk (x)dx,(8)a

We refer to sk as the effective selection coefficient against an allele
in class k, because it is the rate at which any particular lineage
in class k loses individuals (note we have defined signs such thatsk < 0). Note that sk depends implicitly on the hk through theterm involving k̄ (recall k̄ is the average value of k, k̄ ≡khk ).
For convenience we will define the scaled effective selection
coefficient γk by

γk = Nsk .

147

where
fk (x) =(9)This is equivalent to saying that the probability that there exists a
lineage in class k with frequency between x and x + dx is fk (x)dx, for
infinitesimal dx. Note that this PRF result implicitly assumes that θkand γk are constant (which requires constant hk ), and hence only
describes the diversity in steady state.
This PRF description offers a convenient and well-established
way to describe the lineage structure. It is similar in spirit
to the diffusion result used by Ewens (1972) in his originalcomputation of the neutral ESF. However, there is an importantdifference: Ewens’ f (x) was derived as the solution to the diffusion
approximation to the K -allele Wright–Fisher process, in the limitof infinite alleles. This explicitly constrains all lineages to add toa total frequency of 1. The PRF does not impose this constraint.
This makes it possible to compute a simple analytical expression
for fk (x) in the presence of selection. However, it does involvean implicit approximation. In the Supplementary Appendix B, we
describe this approximation along with a way to relax it using analternative branching process model to describe lineage structure.
3.4. The self-consistency condition
It is clear from our PRF formulation above that the allelicdiversity within each fitness class depends on the θk and γk , which
in turn depend on the hk . Yet the sum of the frequencies of all
the alleles within fitness class k is, by definition, hk . In steadystate, these two quantities must be equal. Verifying under what
conditions these quantities are equal allows us to determine inwhat parameter regime the PRF formulation is self-consistent.
More specifically, we have derived the steady state value of hk
in Eq. (2),hk =e−Ud /sk!

(7)These neutral mutations are also reflected in Fig. 1b. Note thatfor all k, neutral mutations tend to increase θk , and make sk morenegative. In the presence of neutral mutations, even s0 is negative.
We have seen that new lineages are founded within fitnessclass k at rate θk /2, and then drift randomly subject to an effective
selective pressure sk . We now make the key assumption that eachθk
1 − e−2γk (1−x).
x(1 − x) 1 − e−2γk



Ud

k

s

.

When we plug these hk into our PRF result, the summed allele
frequencies according to the PRF must agree with steady-statevalue we used for hk , for consistency. According to our PRF result,the sum of the frequencies of all the alleles in fitness class k is

 1
hk =
0

xfk (x) dx.

(10)

148

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

Because Eq. (2) is equivalent to requiring θk /2 = |γk |hk for all k(i.e. in steady state the net influx of individuals into a class must
equal the average rate at which individuals within that class arelost), we can rewrite the self-consistency equation asθk=
2|γk |

 1
x·
0

θk
1 − e−2γk (1−x)
dx.
x(1 − x) 1 − e−2γk(11)Some algebra reduces this to the condition

 10

1 − e−2γk xxdx =1 − e−2γk
2|γk |.

(12)

The analysis in Appendix A shows that this condition holds to thelevel of approximation considered whenever |γk | ≫ 1. When thisis true, the steady state mutation-selection balance of Eq. (2) isalso the distribution hk that makes our PRF analysis of the allelicdiversity within each fitness class self-consistent.
The condition |γk | ≫ 1 corresponds to saying that the effective
selection coefficient in each class is large compared to 1/N. This
will be true for all k whenever NUn ≫ 1. In practice, even when
this condition fails in some fitness classes, it is still valid for allclasses in which |γk | ≫ 1. Thus our results still give a good
approximation to the population allelic diversity provided |γk | ≫1 for the classes around k̄ that make up the bulk of the population.
This will hold whenever γk̄ = N (Ud + Un ) ≫ 1. When thiscondition does not apply, our PRF result for the allelic diversity
within each fitness class is inaccurate. This is because, when |γk | ≫1, the growth of some mutant lineages is limited by the size ofthe population, which is ignored by the PRF approximation. Thus
the PRF approximation overestimates the probability that lineages
become common, and the self-consistency breaks down.
It is important to note that we also require an additional,
stronger condition for other aspects of our analysis to be valid.
The self-consistency condition ensures that the average size of thefitness class implied by the PRF analysis equals the steady state
hk . However, even when this holds, there could be substantial
fluctuations in hk around its average value. The PRF result for fk (x)tells us the probability that a set of lineages exists at any given
frequencies. Therefore it contains detailed information about thesefluctuations. However, we have neglected these fluctuations in
substituting the hk into our expressions for θk and sk , and will also
neglect these fluctuations below in calculating sampling formulas.
We return to consider this additional approximation in a later
section.
3.5. An alternative, retrospective approach
It is possible to derive the neutral Ewens sampling formulain two quite different ways. Ewens (1972) imagined new alleles
being created continuously by new mutations, and considered
the frequency distribution of lineages set up by the balance
between the continual creation of new alleles and the extinctionof older alleles. This leads to expressions analogous to those in
our PRF calculation of the lineage structure. We can calculate
sampling formulas from this lineage structure, as Ewens did in
the neutral case. First, however, we note that in a companion
paper to Ewens (1972), Karlin and McGregor (1972) showedthat the Ewens sampling formula could also be derived using a
retrospective analysis, by considering the ancestral history of asample of individuals. This same type of retrospective approachis also possible in our model; in this section we describe thisalternative derivation of the allelic diversity as relevant to the caseof purifying selection.
In order to calculate the probability of a particular allelicconfiguration, we consider the ancestral history of a sampled setof individuals. In particular, we are interested in the most recent

event to occur in the history of a sample, backwards in time. We
classify these possible events into one of three possible types:
coalescence events (i.e. identity by descent), neutral mutations,and deleterious mutations.
This method is easiest to understand if we begin by consideringa sample of size two. In order for two individuals to have thesame genotype, they of course must be in the same fitness classk. Furthermore, if we look at the ancestral history of each of these
two individuals, the most recent event to occur, backwards in time,must be a coalescent event. In contrast, for them to have a different
genotype, the most recent event to occur must be a mutation event.
Therefore, to calculate the probability of either configuration, we
need only calculate the probability that the most recent event is acoalescent event.
In order to calculate the probabilities of each possible mostrecent event, we must know the distribution of times until eachtype of event. In general, neutral mutations are exponentiallydistributed with rate Un per generation. Assuming the steadystate values for hk , deleterious mutations are also exponentiallydistributed with rate sk per generation (Hudson and Kaplan, 1994).
Finally, within  each class, coalescence occurs as a neutral processwith ratei2

per Nhk generations. Therefore, for a sample of size

2, each of which are sampled from class k, we have that:
P (1st Event : Coal.) ∞=dtP (Coal at t )P (No Neut. Mut by t )0× P (No Del. Mut. by t ) ∞dte−t e−2Nhk Un t e−2Nhk skt

=
0

=

1
1 + 2Nhk (Un + sk)=

1
1 + θk,(13)

where we have defined θk ≡ 2Nhk (sk + Un ). Of course, this result
agrees with the standard neutral result, replacing θ by θk (seebelow).
This same logic can be easily extended to larger sample sizes.
For example, if we consider i individuals within the same class, theprobability that the first event is a coalescence event is
P (1st Event: Coal.) = ∞dtP (Coal at t )0× P (No Neut. Mut by t )P (No Del. Mut. by t ) ∞   i
=

dt
0

i

2

e

−

2 t

e−iNhk Un t e−iNhk skt

 
=  
i
2

i
2

+ iNhk (Un + sk)=

i−1
i − 1 + θk

.

(14)

If the first event is a coalescence event, that means two of theindividuals are of the same allelic type. This leaves us with i − 1
individuals in the class which may or may not be identical; we
can now use the identical method to ask whether any of theseremaining individuals are of the same allelic type. Similarly, if thefirst event is a mutation event, the remaining i − 1 individuals could
still coalesce with each other before they also experience mutation
events.
We note that our analysis in this section is very similar in
spirit to that of Hudson and Kaplan (1994), Barton and Etheridge(2004), and particularly to Gordo et al. (2002). These earlier authors
considered the relative probabilities of mutations and coalescencein the ancestry of each individual, leading to expressions that
implicitly contain results analogous to those in this section. They
did not however consider the implications of these results for theoverall patterns of allelic diversity in the population, which we now
turn to.

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

3.6. Sampling formulas
We can now calculate the probability of sampled configurationsof allelic types. Our goal is to calculate the probability that asample of n individuals will have some distribution of allelic types(e.g. n1 individuals with allele 1, n2 individuals with allele 2,etc). Specifically, we aim to calculate a negative selection versionof the neutral Ewens sampling formula (ESF). As we will see,this calculation proceeds exactly analogously whether we use thelineage structure (PRF) or retrospective analysis.
We begin with the simplest case, a sample of n = 2 individuals from the population. What is the chance that these individuals
are the same genotype? In other words, what is the allelic homozygosity, Q2 , in the population? In order to be the same genotype, thetwo individuals must carry the same number of deleterious mutations – i.e. they must fall in the same Hamming class, k. In addition,
they must also be of the same mutant lineage within class k. This
must equal the probability that the most recent event in the historyof these 2 individuals is a coalescence event; from Eq. (14) this is1. Alternatively, we could calculate the probability the two in1+θk
dividuals are in the same lineage directly from our PRF result; it isthe expected value of x2 , where x is integrated over the distributionof lineage frequencies in class k:

 1

x

0

h2k

2

fk (x)dx =1
1 + θk,(15)

where we have evaluate the integral as described in Appendix A(see also the corrections in the Supplementary Appendix B).
We therefore find that the full probability that two sampledindividuals have the same genotype, which we denote Q2 , is given
by

Q2 =

∞


h2k





1
1 + θk

k=0

.

(16)

Note that, in the case Ud = 0, all individuals are in the zero class,such that hk̸=0 → 0 and h0 → 1. Therefore:
Q2Neutral →1
1 + 2NUn,(17)

in agreement with the neutral Ewens sampling formula.
In order for two individuals to have a different genotype, there
are two possibilities: either the two individuals could be sampledfrom different classes (in which case they must have a different
genotype), or they could be sampled from the same class, and be ofdifferent allelic types (cf. the first event in their ancestral history is
a mutation event). Therefore:
Q1,1 =

hk hk′ +k,k′ ̸=k

= 1−

h2k



k



h2k



k

1
1 + θkθk1 + θk= 1 − Q2 .

(18)

Note that:
Q1Neutral,1→

3.6.1. Relationship with the neutral result
At this point, it is informative to consider the form of thisresult. The presence of selection serves to subdivide the population
into classes, as given by the mutation-selection balance result.
Thus, in order for a sample of individuals to have a particularallelic configuration, they must be sampled from a set of classes
consistent with that configuration. However, within each class,the population behaves identically to that of a neutral population,with a different population size (N → Nhk ) and mutation rate
(Un → Un + sk). We can see this explicitly by defining:
Q{ESF
Configuration},k ≡ ESF Result for {Configuration} with

θ → 2Nhk (Un + sk).

1 + 2NUn,in agreement with the neutral Ewens sampling formula.

(19)

(20)

For example, we have that:
Q{ESF2},k =

1
1 + θk,Q{ESF1,1},k =θk.
1 + θk(21)

We can then rewrite our results as:
Q2 =

h2k Q{ESF2},k ,(22)k

Q1,1 =

h2k Q{ESF1,1},k +k

hk hk′ .

(23)k,k′ ̸=kThus we see that, within each class, the probability of a particularconfiguration is effectively neutral with parameter θ = 2Nhk (Un +sk), consistent with our initial intuitive guess for the form of ourresult. The overall probability of a given allelic configuration is thenthe probability that a specific configuration is achieved within eachclass, summed over all possible sets of class configurations that areconsistent with the allelic configuration.
3.6.2. Sample size n = 3This logic can be extended to larger sample sizes. In order for
three randomly-selected individuals to have the same genotype,
all three individuals must be sampled from the same class and theymust all be from the same lineage (i.e. both of the first two events
must be coalescence). This can be computed by considering the
1average of x3 over the PRF, 0 x3 fk (x)dx, or by using the results from
Eq. (14). We find:
Q3 =



h3k

k





2
2 + θk1
1 + θk.

(24)

Note that, for Ud = 0, hk̸=0 → 0 and h0 → 1, such that:
Q3Neutral →2(2 + θ )(1 + θ )

,

(25)

in agreement with the neutral Ewens sampling formula.
In order for two individuals to have the same genotype and thethird individual to have a different genotype – a configuration weterm bizygotic – there are two possibilities. First, two individuals
could have been selected from the same class and the third
individual could have been selected from a different class. In thiscase, the two individuals in the same class must be from the same
lineage (i.e. coalesce prior to a mutation event). Alternatively, allthree individuals could have been selected from the same class. In
this case, two must be from the same lineage and the third from adifferent lineage, which occurs with probability

 1
2NUn

149

3x2 (1 − x)fk (x)dx.

(26)

0

Thinking about this retrospectively, this is equivalent to the sum of
two possibilities: either the first event could be a mutation event,in which case the next event among the other two lineages must be150

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

a coalescent event, or the first event could be a coalescent event, inwhich case the next event among the third lineage and the merged
lineage must be a mutation event. We findQ2,1 =

3h2k hk′

1
1 + θkk,k′ ̸=kθk+h3k
2 + θk1 + θk
k


 3h2
2hk
k
1−
.
=
1 + θk2 + θk
k






2




+

θk
2 + θk

1
1 + θk(27)

Note that:3θ
,(1 + θ )(2 + θ )Q2Neutral
→
,1



Qn =

+

h3k



k

= 1−





θk
2 + θk


3h2k

h3k

θk
1 + θk




1



θk
1 + θk(1 + θk )(2 + θk )= 1 − Q3 − Q2,1 ,(29)

as expected. Note that
Q1Neutral,1,1 =

θ2,(1 + θ )(2 + θ )

hnk



=

hnk



k

n−1



n − 1 + θkkθk +n−1θk

QnNeutral → 





4



n−2
n − 2 + θk




···

1



1 + θk

.

(35)

,

(36)

Note that

1 + θkk

3h2k hk′k,k′ ̸=kk+

hk hk′ hk′′ +k,k′ ̸=k,k′′ ̸=k′ ,k

3.6.4. Sampling formulas for arbitrary sample size
We can extend this method to arbitrary sample size. For
example, in order for a sample of n individuals to each have thesame genotype, all individuals must be sampled from the sameclass. They must all be of the same allelic type, which occurs with
1
probability 0 xn fk (x)dx. Or equivalently, the first event among then lineages must be a coalescent event, the next event among theremaining n − 1 lineages must also be a coalescent event, and soon. We find(28)

in agreement with the neutral Ewens sampling formula for thisconfiguration, which we call bizygotic.
Analogous considerations lead to the probability that all three
individuals are of different allelic types,
Q1,1,1 =θ → 2Nhk (Un + sk). The overall probability of a given allelicconfiguration is then the probability that a specific configurationis achieved within each class, summed over all possible classconfigurations that are consistent with the allelic configuration.

(30)

in agreement with the neutral Ewens sampling formula.
In principle, this method can be extended to calculate theprobability of any allelic configuration. Alternatively, we can use
the relationship between these results and the neutral Ewens
sampling formula to infer the probabilities. We found that, for thecases n = 2 and n = 3, we can write the probability of a given
allelic configuration as the probability that, within each class, aparticular configuration is achieved, summed over all sets of classconfigurations that are consistent with the allelic configuration.
Similarly, we see that for Qn :Qn =in agreement with the neutral Ewens sampling formula.

1
θ+n−1
θ



hnk Q{ESFn},k ,(37)k

3.6.3. Relationship with the neutral result
As before, we define a class-specific version of the neutral
Ewens sampling formula with θ → 2Nhk (Un + sk):
Q{ESF
Configuration},k ≡ ESF Result for {Configuration} with

θ → 2Nhk (Un + sk).

(31)

In particular, we have that:
2
3θk
,
Q{ESF
,
2,1},k =(1 + θk )(2 + θk )(1 + θk )(2 + θk )
θk2
.
Q{ESF
1,1,1},k =(1 + θk )(2 + θk )Q{ESF3},k =

Using these formulas, we can rewrite our results:
Q3 =

h3k Q{ESF3},k ,(32)k

Q2,1 =

h3k Q{ESF2,1},k +kQ1,1,1 =

h3k Q{ESF
1,1,1},kk+3h2k hk′ Q{ESF2},k ,(33)

k,k′ ̸=k+k,k′ ̸=k

Q{ESF
Configuration},k ≡ ESF Result for {Configuration} with

θ → 2Nhk (Un + sk).

(38)

Using this logic, we can infer the probability of additional
configurations. For example, in order to sample n individuals of one
genotype and n − m of another, there are two possibilities: First,
m individuals could be sampled from class k and n − m individuals′
could be sampled from another  class k . The probability of samplingn −m nin this manner is hmh.
Withinclass k, the probability of thek k′mm individuals having the same genotype is given by the neutral′
result Q{ESFm},k with θ → 2Nhk (sk + Un ). Similarly, within class k ,the probability of the n − m individuals having the same genotype
is Q{ESF. Alternatively, all n individuals could be sampled fromn−m},k′the same class k. This occurs with probability hnk . The probabilityof m individuals having the same genotype and n − m individuals
having another is then given by Q{ESFm,n−m},k . Combining these results
and summing over all sets of k and k′ , we have that:
Qm,n−m =3h2k hk′ Q{ESF
1,1},khk hk′ hk′′ .

where we have defined:

hnk Q{ESFm,n−m},kk(34)k,k′ ̸=k,k′′ ̸=k′ ,kTherefore, we again see that, within each class, the probabilities ofa particular configuration are effectively neutral with parameter+k,k′ ̸=k

n−m
hm
k hk′

n
m

ESF
Q{ESF
m},k Q{n−m},k′ .

(39)Note, however, that if m = n − m we must divide by two in thesecond term in the above expression, to avoid double-counting.

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

Extending this logic, we have that:
Qn−m−p,m,p =mean), our approach will still give a good approximation to the
population allelic diversity.
We note that in addition to assuming hk are in their steadystate values in defining θk and sk for both the PRF and retrospective
approaches, the PRF contains an additional implicit approximation.
In writing the PRF sampling formulas, we assumed that, forexample, the probability two individuals in class k come from a

hnk Q{ESFn−m−p,m,p},k

k

+



n−m−p m+p
hk
hk′



n



m+p

k,k′ ̸=k

ESF
× Q{ESFn−m−p},k Q{m,p},k′ p n −p  n ESF
hk hk′+
Q{ESFp},k Q{n−m−p,m},k′

+



n −m
hm
k hk′

n
m

k,k′ ̸=k



+

2

p

k,k′ ̸=k

ESF
Q{ESFm},k Q{n−m−p,p},k′n−m−p m p
hk′ hk′′hk

k,k′ ̸=k,k′′ ̸=k,k′


×
×

n



n − m − p, m, p
ESF
ESF
Q{ESF
n−m−p},k Q{m},k′ Q{p},k′′

.

(40)

Note, however, that we must correct the above expression for overcounting if two or more classes require identical configurations(e.g. if n − m − p = m = p we must divide the second throughfourth terms in the above expression by 3 and the last term by 6).
In general, the probability of any allelic configuration can be written as the sum over all possible class combinations that are consistent with a given allelic configuration, where the probability ofeach configuration within a class is given by the neutral result with
θ → 2Nhk (sk + Un ). In the Supplement we provide a computer
algorithm that performs this sum symbolically, for any allelic configuration Qi,j,k,... .
Note that, in the case Ud = 0, all individuals are sampled fromthe zero-class, such that hk̸=0 → 0 and h0 → 1. In this case,only the leading-order term will be non-zero in the above results.
Therefore, the results reduce exactly to the neutral Ewens sampling
formula.
3.7. Fluctuations in the steady state hk

lineage of frequency x (given that lineage exists) is hx . This assumeskthat hk and x are independent quantities. That is, we assume
that all the lineages in the class always add up to a frequency
hk (i.e., we neglect fluctuations in hk ). However, the existence ofa high-frequency lineage naturally implies that hk is likely to belarger than average, and vice versa.
These correlations between the frequency of an individual
lineage and the hk do not pose a problem to our retrospective
analysis, which never makes reference to lineages, but it does leadto small errors in the PRF results. We show in the Supplementary
Appendix B that these errors are negligible provided thatfluctuations in hk can be neglected (i.e. provided Nhk sk ≫ 1).
However, they do lead to small discrepancies between the PRF
and retrospective results (and between the PRF results and theneutral ESF in the Ud → 0 limit, since the neutral ESF is derived
assuming a strict constraint on the total population size). Thus in
the Supplementary Appendix B we describe a method to correct
for these effects, making the lineage-based and retrospective
approaches to allelic diversity exactly equivalent. All of the above
sampling formulas include this correction, as do all our figures.
As a result of fluctuations in the values of hk , there will also
be fluctuations in the value of the average class, k̄. But these are
negligible in the same situations that fluctuations in hk are.
There is one additional extreme effect of fluctuations in hk : afluctuation in h0 can lead to loss of this most-fit class, a process
referred to as Muller’s ratchet. We expect that, provided the ratchet
does not click many times over the timescale in which individual
lineages exist, this will not significantly affect the allelic diversity.
Thus we have neglected the ratchet in our analysis. We return toconsider this in more detail in the Discussion, and test the validityof our approximation with numerical simulations.
3.8. A distribution of fitness effects of deleterious mutations

Even when the self-consistency condition holds, the frequencies hk will fluctuate about their steady state frequencies. However,both our PRF description of the lineage structure and our retrospective analysis assume that the fitness distribution is always in thesteady state, hk . We have previously studied this approximation in
Walczak et al. (2011). Here we summarize our analysis of the validity of this approximation, as relevant for the present paper.
Each allele in class k can at most contain s1 individuals;kselection prevents any individual allele from becoming morecommon than this. The total number of individuals in the classis on average Nhk . Thus when Nhk ≫ s1 , each fitness classkcontains many individual alleles. Thus we expect that the overall
fluctuations in hk should be negligible provided that this condition
holds. This intuition can be made precise: we can calculate thevariance in hk in steady state from our PRF approach, or moreeasily from a branching process approximation described in the
Supplementary Appendix B. By computing Var(hk )/hk , we showthat in fact the fluctuations in hk are indeed negligible provided
that
Nhk sk ≫ 1.

151

(41)

In practice, this condition will often not hold in the high-fitness(and low-fitness) tails of the distribution. However, provided it
holds in the center of the fitness distribution from which mostindividuals will be sampled (i.e. for those fitness classes near theWe have analyzed a model in which all deleterious mutations
have the same fitness cost, s. However, in most real populations
it is likely that deleterious mutations have a range of possible
fitness effects. We could model this by assuming that the overall
deleterious mutation rate is still Ud , but that deleterious mutations
have a fitness cost between s and s+ds with probability ρ(s)ds. Thatis, ρ(s) is the distribution of fitness effects of deleterious mutations.
In this more general situation, there is still a steady state
distribution of fitness within the population. Generalizing ourearlier notation, we can write this distribution as h(k), where Nh(k)is the steady state number of individuals with a fitness between sk
and (s + ds)k, where s is the average fitness cost of a deleterious
mutation and k is no longer constrained to be an integer. For certain
ρ(s) (e.g. an exponential distribution) it is possible to calculate h(k)analytically, but even when this is not possible there does exist
some steady state h(k).
The basic ideas behind our analysis still apply in this moregeneral situation. The rate at which new lineages within fitness‘‘class’’ h(k) are created is now

θ (k)/2 = Nh(k)Un + N kh(k′ )ρ((k − k′ )/s)dk′ .

(42)

0

The effective selection pressure against individuals in this class is
s(k) = Un + Ud − (k − k̄)s.

(43)

152

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

a

b

c

d

Fig. 2. A comparison between simulation results (dots) and the predictions of our theory (black lines), for the case where some mutations are deleterious and others areneutral. For comparison we also show the predictions of NS interpretation of the neutral Ewens Sampling formula (gray lines; the NM interpretation gives a worse fit to thedata). (a) Homozygosity Q2 as a function of Ud /s for N = 5 × 104 . (b) Q2,1 as a function of Ud /s for N = 5 × 104 . (c) Homozygosity Q2 as a function of N for Ud /s = 6. (d) Q2,1
as a function of N for Ud /s = 6. In all plots Un = 3.2 × 10−4 , s = 10−3 .

Using these modified parameters, we can now apply our analysis
as before; the distribution of lineage frequencies in class k is givenby the PRF formula f (k; x) with appropriate θ (k) and s(k). We canthen find sampling formulas as before—the only difference is thatinstead of summing over a discrete set of fitness classes, we must
integrate over a continuous set of possible fitnesses. For example,
∞1
we have Q2 = 0 0 x2 f (k, x)dxdk.
This extension of our model allows us to calculate the effects
of more general forms of purifying selection on allelic diversity.
However, there is a wide array of possible distributions ρ(s), andusing this more general form obscures the basic effects of selection.
Thus in analyzing our results and comparing to simulations wefocus on the simpler case in which all deleterious mutations havethe same fitness cost s. This focus has the advantage of simplicity,and it allows us to explore more clearly how the strength of
selection affects the patterns of allelic diversity.

mutations are introduce, with mean NUn . The mutations are
distributed randomly and independently among the individuals in
the population (so that a single individual might receive multiple
mutations in a given generation). Each new mutation is ascribed toa novel site, so that each mutation results in a new genotype.

3.9. Simulations

Starting from a monomorphic population, all simulations wererun for at least 1s ln(Ud /s) generations (or for at least several
times N generations when Ud /s < 1), to ensure relaxation bothto the steady-state mutation-selection equilibrium and to the
PRF equilibrium of allelic frequencies within each fitness class.
Appropriate relaxation to steady state was checked by extending
the simulations and ensuring our results did not change. The
final state of the population – i.e. the frequencies of all surviving
genotypes – was recorded at the last generation, and Q2 and Q2,1
were calculated from these frequencies. This was repeated andaveraged over 250 replicate simulations to produce the points
shown in the figures.

In order to check the validity of our analysis, we have performed
simulations of a Wright–Fisher population. In our simulations,
we consider a population of constant size N and keep trackof the frequencies of all genotypes over successive, discrete
generations. In each generation, N individuals are sampled withreplacement from the preceding generation, according to thestandard Wright–Fisher process (Ewens, 2004) in which the chanceof sampling an individual is determined by its fitness relative to the
population mean fitness.
In each generation, a Poisson number of deleterious mutations
are introduced, with mean NUd , and a Poisson number of neutral

Our simulations allowed for random fluctuations in thefrequencies of each fitness class, as well as for Muller’s ratchet. The
ratchet did not proceed substantially for the simulations relevant
for Fig. 3, except for the highest Ud point shown in that figure.
However, it did proceed substantially in the simulations shownin Fig. 2, such that the most-fit individuals at the end of eachsimulation contained typically a few (for small Ud /s) to more thana dozen (for larger Ud /s ∼ 10) deleterious mutations. We can see
that, despite the effects of Muller’s ratchet and fluctuations in the
hk , our simulations are generally in excellent agreement with ourtheoretical predictions.

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

a153

small samples. Our aim is to highlight the essential differences
between neutral diversity and the diversity in situations involving
linked deleterious mutations.
Aside from likelihoods of configurations, and associated statistics, our approach could also be used to calculate the full distribution of branch lengths, following the generating function approach
used by Lohse et al. (2011).
4.1. Relationship to the neutral Ewens sampling formulab

Fig. 3. Allelic diversity as a function of ln Ud , for Un = 10−4 , s = 10−3 , and N =5 × 104 . Our predictions are shown as a solid line, compared to the predictions of the
NS-ESF (dotted line) and NM-ESF (dash-dotted line). We also compare our results
to the predictions of a neutral ESF using the effective population size that would be
predicted by background selection (EPS, dashed line), though we emphasize this is
not the situation the BGS approximation was developed to address. These analyticalpredictions can be compared to simulation results (dots). (a) Homozygosity Q2 . (b)Q2,1 . Note that Q3 ≈ 0 everywhere for these parameters, so for these predictions
Q1,1,1 ≈ 1 − Q2,1 .

4. Results and discussion
Using the approach we have described, we can calculate theprobability of any allelic configuration within a sample of nindividuals from a population experiencing negative selection
at many linked sites. From this, we can calculate the expecteddistribution of any statistic describing allelic diversity. To doso we must first determine which allelic configurations lead towhat values of the statistic. The probability of each possible
value of the statistic is then the sum of the probabilities of all
allelic configurations leading to that value. This is identical to thecalculation we would do in the neutral case—the only differenceis that to calculate the probability of each allelic configuration, we
use our sampling formula rather than the neutral Ewens sampling
formula.
In practice, some statistics are easier to calculate than others.
While we can easily calculate the distribution of statistics
describing diversity in a small sample, and we could in principle
calculate certain statistics in larger samples (e.g. the total numberof alleles in a sample of size n, Kn ), further work is needed todevelop efficient methods of calculating arbitrary statistics inlarge samples. This is clearly important for applications of our
method to analysis of sequence data, but the combinatoric andcomputational issues involved are an extensive topic which istangential to the ideas underlying our method. Instead, we focus
here on describing the distributions of simple statistics involving

Although it may seem counterintuitive, our analysis applies
even when Ud = 0 (that is, in the case where all mutations areneutral). In this case, our model is the same as that studied by
Ewens (1972). If we apply our methods to this Ud = 0 case, allgenotypes are in the fitness class k = 0, and we have h0 = 1,
γ0 = −NUn and θ0 = θ = 2NUn . Provided that |γ0 | ≫ 1, theconditions for our PRF analysis to be valid are met, and all of ourprevious results still apply, but are greatly simplified. And from ouranalysis of sampling formulas above we can immediately see that,as expected, setting Ud = 0 always causes our results to exactly
reduce to the neutral Ewens sampling formula. Note that we musttake the limit Ud → 0 rather than s → 0 to recover the neutral
result, because taking s → 0 with finite Ud causes the steady state
mutation-selection balance to break down (i.e. we have hk → 0
and fluctuations in the frequencies of each class become crucial).
For nonzero Ud , we expect that our results will differ from thepredictions of the neutral ESF. To illustrate these differences inmore detail, we study the allelic configurations in samples of sizen = 2 and n = 3. Consider first the homozygosity Q2 in a sampleof size n = 2. In Fig. 2a and c we show how Q2 depends on Ud andthe population size N, both under our theory and in Monte Carlosimulations. We compare these results with the predictions of theneutral ESF. We make the same comparisons for the heterozygosity
Q2,1 in Fig. 2b and d. We note that the simulation results agree wellwith our predictions and differ from those of the ESF.
In making this comparison, there is some ambiguity about how
to interpret the ESF, which depends only on θ , for Ud > 0. In
one interpretation, we neglect selection against the deleterious
mutations and set θ = 2N (Un + Ud ); we refer to this as the NSESF case. Alternatively, we could neglect the deleterious mutations
entirely and set θ = 2NUn ; we refer to this as the NM-ESF case.
In Fig. 3 we explore the ambiguity in the interpretation of
the ESF, and compare the predictions of our theory to the twodifferent interpretations of the ESF. For small Ud , our prediction
is equivalent to both interpretations of the neutral ESF. As Ud
increases, our predicted homozygosity decreases slowly until it
experiences a sharp transition at Ud ≈ s. This transition makes
intuitive sense: when Ud < s, most individuals in the population
have no deleterious mutations, and hence the allelic diversity issimilar to the neutral case. As Ud increases past s, most individuals
have deleterious mutations, so these mutations decrease theexpected homozygosity. These deleterious mutations decreasehomozygosity by less than they would if they were neutral, so
our predicted homozygosity is higher than the NS-ESF (neglecting
selection against deleterious mutations) but lower than the NMESF (neglecting deleterious mutations entirely).
We can gain further insight into this behavior by comparing
our predictions to those of the NS-ESF and the NM-ESF in moredetail (Fig. 3). We see that even when Ud = Un , our predictedhomozygosity is only slightly lower than when Ud = 0, despite
the fact that there are twice as many mutations occurring (andhence the NS-ESF prediction for Q2 has declined by a factor of two).
Here the NM-ESF prediction is fairly accurate, reflecting the factthat selection is still strong (with Ud ≪ s) so that most individuals
have no deleterious mutations at all. However, as Ud increases pasts, most individuals now have one or more deleterious mutations

154

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

and hence these mutations decrease our prediction for the allelichomozygosity. In this regime, the NM-ESF becomes inaccurate,
because the deleterious mutations are sufficiently weakly selected(Ud & s) that their presence is important to the diversity. However,despite this being weak selection, the fact that selection eliminates
deleterious mutations from the population more rapidly than ifthey were neutral means that the allelic homozygosity is higher
than the NS-ESF, even as Ud becomes very large. As Ud increases,
our predictions become more similar to the NS-ESF, and in thelimit of infinite Ud will equal the NS-ESF. In Fig. 3b we show the
bizygosity Q2,1 as a function of Ud . Through this parameter range Q3
is small, and so Q1,1,1 ≈ 1 − Q2,1 . As Fig. 3b shows, the dependenceof bizygosity on Ud is similar to the behavior of heterozygosity, foressentially the same reasons.
This shift in our results from being approximately equal to the
NM-ESF for small Ud to the NS-ESF for large Ud has an intuitive
explanation from the form of our results for θk . For Ud ≪ s, h0
is close to 1, since most individuals have no deleterious mutations.
In this class, we have θ0 = 2Nh0 s0 ≈ 2NUn , the same as the θ forthe NM-ESF. Since diversity within each class is neutral with theappropriate θ , in this Ud ≪ s regime the diversity is approximatelythat predicted by the NM-ESF. On the other hand, in the limit ofvery large Ud , hk becomes sharply peaked about k = Ud /s, so
almost all individuals have approximately the same fitness, andindividual deleterious mutations change fitness by a negligible
amount. Thus the diversity is approximately that predicted by the
NS-ESF. This behavior is exactly as reflected in Fig. 3, with thetransition between the two regimes occurring at Ud ∼ s, as thisanalysis would predict.
Our analysis above makes it clear that the difference between
weak and strong selection for the purpose of allelic diversity is set
by whether s is small or large compared to Ud . We have potentiallythree regimes of selection strength. For Ns < 1, selection isineffective relative to drift, and we always have nearly neutral
diversity. For Ns > 1, we can have weak, moderate, or strong
selection. When s ≪ Ud , we have weak selection as describedabove; the NS-ESF is accurate. When s . Ud , we have a ‘‘moderate
selection’’ regime where the diversity generated by the deleterious
mutations themselves can be important, and hence the NM-ESF
is inaccurate. However selection is not so weak that the NS-ESF
is accurate either; the selection against the deleterious mutations
does reduce the amount of diversity they contribute. In this regime,
neither interpretation of the Ewens neutral sampling formulaprovides an accurate prediction for allelic diversity. Finally, fors ≫ Ud , we have a ‘‘strong selection’’ regime, where deleterious
mutations are eliminated quickly from the population and hence
do not contribute to diversity, and the NM-ESF is accurate. The
NS-ESF is also accurate in this regime when Ud ≪ Un but it will
underestimate homozygosity when Ud & Un . Note that in Fig. 3 weshow a case where s > Un , so there is a regime where s ≫ Ud but
Ud & Un and hence the NM-ESF is accurate but the NS-ESF is not.
Such a regime does not exist in the case s < Un , but otherwise thesame qualitative patterns exist for the same reasons.
4.2. Comparison to the effective population size approximation
The background selection model we have studied has been
the subject of much earlier work, although this has largely been
focused on the structure of genealogies in the presence of purifying
selection, rather than allelic diversity (Hudson and Kaplan, 1994,
1995; Gordo et al., 2002; Seger et al., 2010). A particularly simple
and useful approximation to the effects of background selection
was developed by Charlesworth et al. (1993); Charlesworth (1994);
Charlesworth et al. (1995). This approximation is widely used tosummarize the effects of background selection (Hartl, 1988). We
refer to it here as the effective population size approximation(EPS). The EPS analysis makes predictions about the structureof genealogies and hence about genetic diversity at the levelof individual sites, not just the allelic diversity we consider
here. Further, it focuses on the genetic diversity among neutral
mutations only. Thus it is not directly comparable to our resultsin this paper. Despite this, we find it instructive to briefly examine
how EPS compares to our results, if we apply it to predict allelicdiversity. We stress that this is not the interpretation intended
by Charlesworth et al. (1993) and does not provide a fair pictureof its accuracy in general. Since EPS describes the structure ofgenealogies, we defer a detailed discussion of the accuracy of the
EPS approximation and its relationship to our results to Walczak
et al. (2011), where we calculate the structure of genealogies under
our model.
The EPS approximation assumes that deleterious mutations
are eliminated by selection quickly compared to the coalescence
time between two individuals who do not have any suchmutations. When this is true, almost all neutral mutations weobserve occurred in individuals that did not have any deleterious
mutations, because they have little time to occur in individuals
that do have deleterious mutations before these individuals are
eliminated by selection. Thus, according to the EPS approximation,
the genetic diversity among neutral sites linked to negatively
selected sites is exactly the same as the entirely neutral case, butwith the population size N replaced by the size of the least-loaded(i.e. most-fit) class. That is, N is replaced by the effective population
size
Ne = Nh0 = Ne−Ud /s .

(44)

Given this Ne , EPS predicts that any properties of neutral diversity
are identical to those of coalescent theory with the appropriate
Ne . Applying this to the allelic diversity, this predicts that the
sampling properties of neutral alleles will be given by the classicalEwens’ sampling formula, using θ = 2NUn h0 = 2NUn e−Ud /|s| . Note
this is effectively a NM-EPS case, which seems most natural. An
alternative NS-EPS case can be defined using θ = 2N (Un + Ud )h0 ;this leads to similar conclusions.
In the strong selection regime where Ud ≪ s, most individuals
are in the 0-class. Thus our analysis predicts that this class will
dominate allelic diversity, which will be neutral with θ0 =2Nh0 s0 = 2Ne−Ud /s Un . Thus our analysis reduces exactly to thepredictions of the NM-EPS in this regime. This is the regime inwhich the EPS approximation is expected to hold (Walczak et al.,2011), so our analysis reduces to the EPS in the regime in which it
should.
However, for the moderate and weak selection regimes, Ud &s, the EPS prediction breaks down dramatically, consistent withthe earlier observations of Nordborg et al. (1996) and Kaiser andCharlesworth (2009). We graph this prediction in Fig. 3 (using
the NS interpretation of the EPS, which provides a slightly better
prediction than the NM interpretation). In this regime the EPS
predicts that the neutral homozygosity increases dramatically,
since the least-loaded class becomes negligible in size. However,
the homozygosity is not so large in reality, as our predictions
demonstrate. Rather, both neutral and deleterious variation amongindividuals that harbor one or more deleterious mutations isimportant. Our theory accounts for this effect, while EPS fails
because the approximation that the coalescence time between
individuals is dominated by the time in the least-loaded classbreaks down.
We note that, contrary to the intuition one might be tempted
to draw from EPS, having more deleterious mutations can neverdecrease allelic diversity. That is, if we fix all other parameters,
simply having more deleterious mutations (i.e. increasing Ud )does not reduce heterozygosity. Certainly it reduces neutral
heterozygosity, but accounting for all variation a population
with a larger deleterious mutation rate will have more allelicheterozygosity.

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

4.3. Distortions in allelic diversity
The above discussion makes clear that for given population
sizes, mutation rates, and selection strengths, purifying selectionchanges the probabilities of particular allelic configurations in asample. However, this does not necessarily imply that selection
leads to distortions in the patterns of genetic variation compared to
the neutral case. In the neutral case, the probabilities of all allelicconfigurations in a sample are determined by a single parameter θ .
This means that we can infer θ from a statistic which depends on
the probabilities of one set of allelic configurations, and this θ thenpredicts the expected distribution of all other statistics describing
genetic variation within the population, provided it is evolving
neutrally.
Our discussion of the EPS approximation above makes clear thatfor sufficiently strong selection, genetic diversity is not distorted
relative to the neutral case. In this section, we show that formoderate to weaker selection (relative to mutation rates), there
is no effective population size Ne which can describe genetic
diversity in our model. As we noted in the Introduction, this isconsistent with earlier observations that background selection
leads to distortions in the structure of genealogies (McVean andCharlesworth, 2000; Seger et al., 2010; O’Fallon et al., 2010;
Comeron and Kreitman, 2002; Comeron et al., 2008; Barton andEtheridge, 2004; Gordo et al., 2002; Hermisson et al., 2002;
Williamson and Orive, 2002). Here we compute precisely how
these distortions alter particular aspects of the patterns of allelicdiversity. Our analysis in this section demonstrates one place inwhich statistical power exists to distinguish purifying selection
from neutral processes at a reduced effective population size.
Our framework can in principle be used to explore where suchstatistical power lies more generally, but we leave this moregeneral question for future work.
In this section, we simply show that there is no effective neutral
population size Ne to describe diversity in our model. To do this, it
is sufficient to show that the effective θ that one would infer fromone statistic predicts the incorrect values of other statistics. The
simplest way to do this is to begin with the Q2 we would predict
given some set of parameters. We calculate the effective θe one
would infer from this Q2 using the neutral ESF (i.e. we choose θe
1such that Q2 = 1+θ
). We then calculate the neutral prediction fore
Q2,1 (or Q3 ) based on this θe . We compare this with our predictions
for Q2,1 (or Q3 ) given the real parameters. The difference between
these two predictions is a measure of the deviation from neutrality.
We show this deviation from neutrality, expressed as the ratio of
the neutral effective population size prediction to the actual result,for Q2,1 in Fig. 4a and for Q3 in Fig. 4b.
We see from Fig. 4 that negative selection distorts the allelicdiversity away from high-frequency polymorphisms and toward
lower-frequency polymorphisms, for a given level of overall
heterozygosity. The effects are strongest when Ud is of order (orslightly larger than) s, and the distortion is stronger for smaller Un
and N.
These two simple statistics measuring deviations from neutrality demonstrate that there is no effective population size describing allelic diversity. These particular comparisons are presumably
not the most statistically powerful way to detect this type of negative selection, but they do show that statistical power exists. Using the framework developed in this paper, it is now possible tosystematically investigate exactly how linked negatively selectedsites generate different patterns of allelic diversity from the neutral case, and to determine which statistics provide the most power
detect this type of selection. Note for example that the deviationfrom neutrality is much stronger in Fig. 4b than in Fig. 4a. This
reflects the fact that we are inferring θ from Q2 , which in ourtheory is more closely related to Q2,1 than it is to Q3 . Even more

155

a

b

Fig. 4. The deviation from neutrality. We take Q2 as predicted by our theory, and1
use the neutral ESF to find the effective θ that this implies by setting Q2 = 1+θ
.
e
We then use this effective θe in the neutral ESF to predict the values of Q2,1 and Q3 itcorresponds to. We compare this to the Q2,1 and Q3 predicted by our theory. This isa measure of the deviation from neutrality, the skew in the frequency spectrum ofallelic diversity away from neutral results with some modified effective population
size. (a) The ratio of Q2,1 from the effective population size description to the Q2,1
from our theory, as a function of ln(Ud ), for s = 10−3 and three different values ofUn and N. (b) The ratio of Q3 from the effective population size description to theQ3 from our theory as a function of ln(Ud ), for s = 10−3 and three different valuesof Un and N.

powerful tests for selection are presumably possible. While muchearlier work has anticipated that purifying selection distorts the
structure of genealogies (McVean and Charlesworth, 2000; Gordo
et al., 2002; Hahn, 2008; Comeron et al., 2008; Seger et al., 2010;
Betancourt et al., 2009; Comeron and Kreitman, 2002; Williamson
and Orive, 2002), no analytic formalism has previously provided a
way to determine precisely how selection alters patterns of allelicdiversity (and hence, where statistical power may lie).
While we have shown that there is no neutral effective
population size describing allelic diversity, this allelic diversity isa summary statistic of the full per-site diversity. Thus our result
also implies that genetic diversity at a per-site level also cannot be
described by a neutral effective population size, and that additional
power to distinguish neutrality from negative selection can befound in data on site-based variation, consistent with the earlier
work described above.
4.4. Muller’s ratchet
Throughout our analysis, we have assumed that Muller’s ratchet
can be neglected. This is clearly not true in general. The problem
Muller’s ratchet creates is that hk can change with time, and thischanges the distribution of allele frequencies within each class.

156

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

After a ‘‘click’’ of the ratchet, the distribution of hk shifts, eventuallyreaching a new state shifted left by one class (so the class that was
originally at frequency hk is now at frequency hk−1 , and so on). The
PRF distribution of lineage frequencies in class k correspondingly
shifts from fk to fk−1 , and so on, which changes the allelic diversity.
Fortunately, since fk (x) is similar to fk+1 and fk−1 , this effect isunlikely to cause major inaccuracies, provided the ratchet does
not click many times over the timescale on which the lineage
frequency spectrum turns over. We expect that this is generally
true within the bulk of the fitness distribution. At the tails of the
distribution, where hk is small, the allele frequency distribution
can sometimes be substantially different than expected due to theratchet. However, by definition these classes represent a small
fraction of the overall population and hence we do not expect themto contribute substantially to allelic diversity.
We tested the accuracy of our approximation neglecting
Muller’s ratchet using the simulations described above, all ofwhich included the possibility of the ratchet. Our predictions
remain very accurate, even in simulations in which the ratchet wasobserved to operate. Note, however, that the ratchet is potentiallymore problematic in considering the genetic diversity at the levelof individual sites, because the high-fitness tail of the fitness
distribution can be important for the structure of genealogies even
if it does not contribute substantially to allelic diversity at any time.

Science and Engineering Graduate Fellowship Program, and also
acknowledges support from an NSF graduate research fellowship.
Many of the computations in this paper were run on the
Odyssey cluster supported by the FAS Sciences Division Research
Computing Group at Harvard University. AMW thanks the
Princeton Center for Theoretical Science at Princeton University,where she was a fellow during some of her work on thispaper. JBP 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), andthe US National Institute of Allergy and Infectious Diseases(2U54AI057168).

4.5. Conclusionf ( x) =We have introduced a formalism to calculate the statistics ofallelic diversity in the presence of purifying selection at many
linked selected sites. We have done so by calculating the structureof the individual lineages that maintain the deleterious mutationselection balance. This analysis is based on the PRF frameworkof Sawyer and Hartl (1992), which was originally developed todescribe the frequency of mutations at completely unlinked sites.
We have adapted this framework to our problem with a shift
in perspective: rather than treating new mutations at individualsites as the basic and independently fluctuating quantities, weconsider the lineages founded by new mutations as the basic
independent quantities. This allows us to describe aspects of thegenetic diversity despite the fact that selection is acting on many
linked non-independent sites. We showed that this approach isexactly equivalent to a retrospective perspective, which studied
the probability individuals are in the same lineage by considering
the probability that coalescence events preceded mutations.
Of course, each lineage we describe contains many different
mutations, and the fluctuations in lineage frequency described by
the PRF framework represent correlated fluctuations in all of these
individual mutations. If we could also describe how lineages are
related to each other, and hence the statistics of which mutations
they share, we could combine this with the results in this paper
to describe the full per-site patterns of genetic diversity despite
the correlations between sites introduced by linkage and selection.
In this paper, however, we have focused on describing allelicdiversity, leading to a negatively selected version of the neutral
Ewens sampling formula. This analytical framework allows us tocompute precisely how patterns of allelic diversity are distorted
by negative selection at many linked sites, and hence understandexactly where statistical power may lie to distinguish purifying
selection from neutrality.

where we have suppressed the subscripts and used the notationa ≡ −2γ .
Whenever n and m are both ≥ 1, these integrals are easy toevaluate analytically. When either n or m equals zero, the integrals
can be separated into an exactly solvable analytical part and a partthat involves the integral

Acknowledgments
We thank Warren Ewens and Isabel Gordo for helpful
discussions and comments on the manuscript. MMD acknowledges
support from the James S. McDonnell Foundation. LEN is supportedby the Department of Defense through the National Defense

Appendix A. Integrals involving fk (x)Our expressions for the probabilities of various allelic configurations involve integrals of the form

 1
I =

A(x)f (x),

(45)

0

where A(x) is a polynomial function of the form A(x) = xn (1 − x)m(with n and m integers). Here f (x) is the expression from Eq. (9),

I′ =

1

ah
ea

− 1 x(1 − x)

 1eay − 10y

 a(1−x)
e


−1 ,

dy.

(46)

(47)

This integral I ′ is a known special function Ein(−a); see p. 228 of
Abramowitz and Stegun (1965).
Consider for example the integral

 1
I2 =

x2 f (x)dx.

(48)

0

Substituting in for f (x) and substituting y = 1 − x in the integral
gives
I2 =

ah
ea

−1

 1
0


1 − y  aye − 1 dy.
y

(49)

1−y

We now simply write y = 1y − 1 and evaluate the analyticallysolvable parts of this integral to get
I2 =

ah
ea

−1

I′ − h +

ah
ea

−1

.

(50)

Fortunately, we can calculate a simple analytic approximation
for I ′ in the limit a ≫ 1 (i.e. |γ | ≫ 1), which is the limit we arealways working in. This is a standard asymptotic expansion of theEin function; we have
I′ ≈





1
1 a
e 1+
.
a
a

(51)

We can now plug our approximation for I ′ into our result for I2 to
get
I2 =

h
a

.

(52)

M.M. Desai et al. / Theoretical Population Biology 81 (2012) 144–157

For more complex integrals, we need to keep higher order terms
in the asymptotic expansion of I ′ . In general, we find

 1
In =

xn f (x) =0(n − 1)!han − 1.

(53)

Similar calculations can be used to find an analogous approxima1
tion for Im = 0 (1 − x)m f (x)dx, but this integral is not necessary
for our purposes in this paper.
These calculations allow us to give simple analytic expressions
for any integrals of the form xn (1 − x)m f (x)dx. Whenever m andn are both ≥ 1, the integrals can be evaluated exactly in termsof elementary functions, and when either m or n are 0 we canuse the above results to provide simple analytic approximations
to whatever precision we require.
Appendix B. Supplementary data
Supplementary material related to this article can be foundonline at doi:10.1016/j.tpb.2011.12.002.
References
Abramowitz, M., Stegun, I.A., 1965. Handbook of Mathematical Functions. Dover,
New York.
Barton, N.H., Etheridge, A.M., 2004. The effect of selection on genealogies. Genetics
166, 1115–1131.
Betancourt, A.J., Welch, J.J., Charlesworth, B., 2009. Reduced effectiveness of
selection caused by a lack of recombination. Current Biology 19, 655–660.
Charlesworth, B., 1994. The effect of background selection against deleterious
mutations on weakly selected, linked variants. Genetical Research 63, 213–227.
Charlesworth, B., Morgan, M.T., Charlesworth, D., 1993. The effect of deleterious
mutations on neutral molecular variation. Genetics 134, 1289–1303.
Charlesworth, D., Charlesworth, B., Morgan, M.T., 1995. The pattern of neutral
molecular variation under the background selection model. Genetics 141,
1619–1632.
Comeron, J.A., Kreitman, M., 2002. Population, evolutionary and genomic consequences of interference selection. Genetics 161, 389–410.
Comeron, J.M., Williford, A., Kliman, R.M., 2008. The hill-robertson effect:
evolutionary consequences of weak selection and linkage in finite populations.
Heredity 100, 19–31.
Etheridge, A.M., Griffiths, R.C., 2009. A coalescent dual process in a moran modelwith genic selection. Theoretical Population Biology 75, 320–330.
Etheridge, A.M., Griffiths, R.C., Taylor, J.E., 2010. A coalescent dual process in amoran model with genic selection, and the lambda coalescent limit. Theoretical
Population Biology 78, 77–92.
Ethier, S.N., Kurtz, T.G., 1987. The infinitely-many alleles model with selection
as a measure-valued diffusion. Stochastic Models in Biology, Lecture Notes in
Biomathematics 70, 72–86.
Ethier, S.N., Kurtz, T.G., 1994. Convergence to fleming-viot processes in the weakatomic topology. Stochastic Processes and their Applications 54, 1–27.
Ewens, W.J., 1972. The sampling theory of selectively neutral alleles. Theoretical
Population Biology 3, 87–112.
Ewens, W.J., 2004. Mathematical Population Genetics: I. Theoretical Introduction.
Springer, New York, NY.

157

Ewens, W.J., Li, W.-H., 1980. Frequency spectra of neutral and deleterious alleles ina finite population. Journal of Mathematical Biology 10, 155–166.
Gordo, I., Navarro, A., Charlesworth, B., 2002. Muller’s ratchet and the pattern ofvariation at a neutral locus. Genetics 161, 835–848.
Griffiths, R.C., 1983. Allele frequencies with genic selection. Journal of Mathematical
Biology 17, 1–10.
Grote, M.N., Speed, T.P., 2002. Approximate ewens formulae for symmetric
overdominance selection. Annals of Applied Probability 12, 637–663.
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’sratchet. Theoretical Population Biology 14, 251–267.
Hartl, D.L., 1988. A Primer of Population Genetics. Sinauer Associates, Sunderland,MA.
Hermisson, J., Redner, O., Wagner, H., Baake, E., 2002. Mutation-selection balance:
ancestry, load, and maximum principle. Theoretical Population Biology 62,
9–46.
Hudson, R., Kaplan, N., 1994. Gene trees with background selection. In: Golding, B.
(Ed.), Non-Neutral Evolution: Theories and Molecular Data. Chapman and Hall,New York, pp. 140–153.
Hudson, R., Kaplan, N., 1995. Deleterious background selection with recombination.
Genetics 141, 1605–1617.
Joyce, P., 1995. Robustness of the ewens sampling formula. Journal of Applied
Probability 32, 609–622.
Joyce, P., Tavare, S., 1995. The distribution of rare alleles. Journal of Mathematical
Biology 33, 602–618.
Kaiser, V.B., Charlesworth, B., 2009. The effects of deleterious mutations on
evolution in non-recombining genomes. Trends in Genetics 25, 9–12.
Kaplan, N.L., Darden, T., Hudson, R.R., 1988. The coalescent process in models withselection. Genetics 120, 819–829.
Karlin, S., McGregor, J., 1972. Addendum to a paper of W. Ewens. Theoretical
Population Biology 3, 113–116.
Kimura, M., Maruyama, T., 1966. The mutational load with epistatic geneinteractions in fitness. Genetics 54, 1337–1351.
Li, W.-H., 1977. Maintenance of genetic variability under mutation and selection
pressures in a finite population. PNAS 74, 2509–2513.
Li, W.-H., 1978. Maintenance of genetic variability under the joint effect of mutation,
selection and random drift. Genetics 90, 349–382.
Li, W.-H., 1979. Maintenance of genetic variability under the pressure of neutral anddeleterious mutations in a finite population. Genetics 92, 647–667.
Lohse, K., Harrison, R.J., Barton, N.H., 2011. A general method for calculating
likelihoods under the coalescent process. Genetics 189, 977–987.
McVean, G.A.T., Charlesworth, B., 2000. The effects of hill-robertson interference
between weakly selected mutations on patterns of molecular evolution andvariation. Genetics 155, 929–944.
Nordborg, M., Charlesworth, B., Charlesworth, D., 1996. The effect of recombination
on background selection. Genetical Research 67, 159–174.
O’Fallon, B.D., Seger, J., Adler, F.R., 2010. A continuous-state coalescent and theimpact of weak selection on the structure of gene genealogies. Molecular
Biology and Evolution 27, 1162–1172.
Sawyer, S.A., Hartl, D.L., 1992. Population genetics of polymorphism and divergence.
Genetics 132, 1161–1176.
Seger, J., Smith, W.A., Perry, J.J., Hunn, J., Kaliszewska, Z.A., Sala, L.L., Pozzi, L.,
Rowntree, V.J., Adler, F.R., 2010. Gene genealogies strongly distorted by weaklyinterfering mutations in constant environments. Genetics 184, 529–545.
Walczak, A.M., Nicolaisen, L.E., Plotkin, J.B., Desai, M.M., 2011. The structure ofgenealogies in the presence of purifying selection: an effective coalescent
theory. Genetics, doi:10.1534/genetics.111.134544.
Watterson, G.A., 1978. The homozygosity test of neutrality. Genetics 88, 405–417.
Williamson, S., Orive, M.E., 2002. The genealogy of a sequence subject to purifying
selection at multiple sites. Molecular Biology and Evolution 19, 1376–1384.

Supplementary Appendix
The Structure of Allelic Diversity in the Presence of Purifying Selection
Michael M. Desai,1 Lauren E. Nicolaisen,1 Aleksandra M. Walczak,2 and Joshua B. Plokin3
1

Departments of Organismic and Evolutionary Biology and of Physics and
FAS Center for Systems Biology, Harvard University, Cambridge, MA

2

CNRS-Laboratoire de Physique Théorique de l’École Normale Supérieure,3

Department of Biology, University of Pennsylvania

Correlations between the size of a lineage and the frequency of the class
In the main text, we asserted that our PRF calculations and our retrospective approach were equivalent, and thatthey reduce precisely to the neutral Ewens sampling formula when Ud = 0 and to the EPS approximation when
Ud /s  1. However, this is only approximately true: our PRF lineage structure calculations as described thus far areslightly different from the retrospective results (and reduce to the neutral ESF or the EPS approximation to leading
order in θ1k ). These discrepancies all stem from fluctuations in the size of the fitness classes hk . In this Appendix,
we explain how these discrepancies arise and describe an alternative analysis of lineage structure which allows us to
avoid them. We note that we have previously analyzed this issue in Walczak et al. [1]; here we largely summarize theresults as relevant for the present paper.
In our PRF calculations for sampling formulae in the main text, we wrote that the probability that two individuals
picked from class k come from the same lineage is
Z 1
0

x
hk

2
fk (x)dx.

(1)

The idea behind this equation is that fk (x)dx is the probability that there exists a lineage in class k at frequency
x, while hk is the total frequency of class k, so the probability that an individual in class k comes from this lineage(given the lineage exists) is x/hk . Similarly the probability two randomly chosen individuals come from this lineage is
R1
x2 /h2k , and summing up over all possible lineages (times the probability each exists) gives 0 x2 fk (x)/h2k dx. Analogous
expressions were constructed for the probabilities of other allelic configurations.
These expressions assume that the fact that there exists a lineage of frequency x in fitness class k does not impact
the expected size of the lineage hk . This is closely related to neglecting fluctuations in hk : we assume that all the
lineages in the class always add up to hk . We saw in the main text that this holds provided that N hk sk  1. As
we would expect, this is the same condition under which we can neglect fluctuations in the hk . We have therefore
restricted our analysis to this parameter regime.
Even within this parameter regime, this approximation introduces small corrections. These arise from the fact that,given that there exists a lineage that is at high frequency (a substantial fraction of hk ), it is likely that hk is larger
than average. In other words, there is a correlation between the size of a lineage and the frequency of the class, so theR 1  2
probability that two individuals picked from a class come from the same lineage is not precisely 0 hxk fk (x)dx.
The part of the integrand corresponding to larger x overestimates the probability two individuals are from the same
lineage, and the smaller x part of the integrand underestimates the probability two individuals are from the same
lineage.
These problems lead to only slight inaccuracies when our self-consistency equation holds. However, neglecting these
correlations obscures the precise relationship between our results and the neutral Ewens sampling formula (and to ourretrospective approach). Further, relaxing this approximation provides a useful comparison of the subtle differences

2
between the assumptions underlying the original neutral Ewens calculation and the PRF method we use here. Thus
we describe here an alternative approach to understanding the lineage structure in a fitness class which allows us to
avoid this approximation, along the lines we have previously described in Walczak et al. [1]. This eliminates the small
discrepancies between our PRF method, our retrospective approach, and (in the Ud → 0 limit) the neutral Ewens
formula.
Our alternative approach uses a continuous-time branching process to describe the dynamics of each lineage, rather
than a diffusion approximation. As before, we imagine that new lineages are created at a rate θk /2. In steady state
there will be some time-independent probability that there are σk total individuals across all the lineages in the class,
P (σk ). Note that on average we must have σk /N = hk , and that P (σk ) contains information on the fluctuations in
the hk . We first compute the generating function for P (σk ),Hk (z) ≡∞X

P (σk )z σk .

(2)

σk =0

To do so, we start by computing the generating function for the probability distribution of the number of individuals
from each lineage, as described by Eqs. (7-9) of Desai and Fisher [2]. We substitute this expression into Eq. (24) of
Desai and Fisher [2] and integrate. We find
Hk (z) ≡∞X

P (σk , t)z

σk

σk =0



sk
≡ hz i =1 − z(1 − sk )σkθk
 2(1−s
)
k

,

(3)

where angle brackets denote expectation values, and we have suppressed the k subscripts. Note that this equation
describes the fluctuations in the size of an individual fitness class: the mean, variance, and higher moments of σkcan be easily computed by taking derivatives of Hk (z). Note also that this calculation is based on a continuous-time
branching process, in which individuals have a different distribution of offspring number than in a 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,we have to take θ → 2θ (as we would in comparing Wright-Fisher to Moran models), as described by Desai and Fisher[2].
Assume there are B sites in the genome. New mutations found lineages in class k due to mutations at a specificθksite at rate 2B. We implement the infinite-alleles assumption by taking the infinite-B limit, in which each distinct
lineage in a given fitness class arises from a mutation at a different site. The generating function for the probability
that there are ni,k mutants at a particular site i in class k is then

skHi,k (z) =1 − z(1 − sk )θk
 B(1−s
)
k

.

(4)

Here again we have taken θ → 2θ to match to the Wright-Fisher model as described above.
We define ni,k as the total number of mutants at site i in class k, so
σk ≡

B
X

ni,k

(5)

i=1

is the total number of individuals in the class (note that on average we expect σk = N hk ). We now imagine that we
sample some number m individuals from class k. The probability that they are all the same allelic type is
* B
+
X nm
i,k
(k)
Qm =
.
(6)
σkm
i=1
nm

m
Note this has the same form as our PRF expression, except we are averaging over σi,km rather than averaging over ni,kkand then dividing by the average σkm . In other words, we are explicitly accounting for the correlations between x and3hk . Analogous expressions hold for all other possible allelic configurations; we simply take the PRF integrals relevant
for that sampling configuration and replace the hk by σk /N , replace the integrals over x with taking an expectation,and sum over all B sites.
(k)
We now wish to compute these sampling probabilities. For example, consider Qm . We can rewrite Eq. (6) using
the identity
Z ∞xm−1 −xσ
1
=
e
dx.
(7)
σm
(m − 1)!
0This identity can easily be verified by integrating the RHS by parts. Using this, and noting that lineages at each of
the B sites are independent, we find
* B
+
Z ∞
X
xm−1 −xσk
(k)
m
Qm =
e
ni,k
dx
(m − 1)!
0
i=1
Z ∞
xm−1
= B
hnm e−xσk idx(m − 1)! 1,k
0
Z ∞
xm−1
−xn1,khe−xni,k iB−1 hnm
idx.
(8)
= B
1 e
(m − 1)!
0
The first expectation value inside the integral can be computed by noting that
he−xni,k θk
1 − sk B(1−sk ).
i = Hi,k (z = 1 − x) = 1 + x
sk(9)−xn1,kDifferentiating this result m times with respect to x results in an expression for hnm
i. Plugging these results
1,k e
in and integrating, taking the limit B → ∞, and neglecting higher order terms in s, we find

Q(k)
m =θ

m−1
X

(−1)j

j=0



(m − 1)!
m−1
1= Qm−1=
θk + jjj=1 (θk + j)1.

θk +m−1θk(10)This result is exactly equivalent to the sampling formula described in the main text and to the result of theretrospective calculation. It is a correction to the PRF result,
Z 1
Q(k)
=
xm fk (x)dx,
(11)
m
0

which is identical only to leading order in 1/θk . Analogous manipulations for other possible sampling configurations
show that the method described in this section generally produces sampling formulae that are exactly equivalent
to our retrospective method, and which represent second-order corrections to the PRF results. All of the sampling
formulae quoted in the main text and shown in the figures incorporate this correction, which appropriately handles
the correlations between the frequency of an individual lineage and the size of the fitness class.

[1] A. M. Walczak, L. E. Nicolaisen, J. B. Plotkin, and M. M. Desai, Genetics In review (2011).
[2] M. M. Desai and D. S. Fisher, Genetics 176, 1759 (2007).

4

Algorithm to compute Q
In the main text we describe how to compute QESF
{Conf iguration},k for an arbitrary allelic configuration. The probabilityof any allelic configuration can be written as the sum over all possible class combinations that are consistent with agiven allelic configuration, where the probability of each configuration within a class is given by the neutral resultwith θ → 2N hk (sk + Un ). Below we provide an algorithm, written in mathematica, that can output a symbolic result
for an arbitrary allelic configuration and compute it numerically.
(* Symbolic Evaluation: *)Clear["Global‘*"];Needs["Combinatorica‘"];(* Insert configuration here: *)indices={2,2,1,1};
partitions=DeleteDuplicates[Sort/@SetPartitions[indices]];
multiplier=Times@@#&/@((Tally[#][[All,2]]&/@partitions)!);products=Table[Product[h[Subscript[k,m]]\^Total[partitions[[j,m]]]*Subscript[Q,{Subscript[k,m],partitions[[j,m]]}]/Total[partitions[[j,m]]]!,{m,Length[partitions[[j]]]}],{j,Length[partitions]}];result=Sum[Sum[products[[j]]/multiplier[[j]]*Total[indices]!,
Evaluate[Sequence@@Table[{Subscript[k,o],Infinity},{o,Length[partitions[[j]]]}]]],{j,Length[partitions]}];(* Numerical Evaluation: *)Clear["Global‘*"];Needs["Combinatorica‘"];(* Insert configuration here: *)indices={2,2,1,1};
Ud=0.004;s=0.001;Un=0.004;n=10\^5;
h[k\_]:=Exp[-Ud/s]*(Ud/s)\^k/k!
Q[k\_,conf\_]:=P[2*n*h[k]*(s*k+Un),Total[conf],BinCounts[conf,{1,Max[conf]+1,1}]]
P[th\_,nn\_,conf\_]:=1/Binomial[th+nn-1,th-1]*Product[(th/q)\^conf[[q]]/conf[[q]]!,{q,Length[conf]}]partitions=DeleteDuplicates[Sort/@SetPartitions[indices]];
multiplier=Times@@#&/@((Tally[#][[All,2]]&/@partitions)!);products=Table[Product[h[k[m]]\^Total[partitions[[j,m]]]*Q[k[m],partitions[[j,m]]]/Total[partitions[[j,m]]]!,{m,Length[partitions[[j]]]}],{j,Length[partitions]}];result=Sum[Sum[If[Length[Tally[Table[k[p],{p,Length[partitions[[j]]]}]]]!=
Length[partitions[[j]]],0,products[[j]]/multiplier[[j]]*Total[indices]!],
Evaluate[Sequence@@Table[{k[o],0,5*Ceiling[Ud/s]},{o,Length[partitions[[j]]]}]]],{j,Length[partitions]}];