Copyright Ó 2010 by the Genetics Society of America
DOI: 10.1534/genetics.110.115030

Some Consequences of Demographic Stochasticity in Population Genetics
Todd L. Parsons,* Christopher Quince† and Joshua B. Plotkin*,1*Department of Biology, University of Pennsylvania, Philadelphia, Pennsylvania 19104 and †Department
of Civil Engineering, University of Glasgow, Glasgow G12 8LT, United Kingdom
Manuscript received February 1, 2010
Accepted for publication May 5, 2010
ABSTRACT
Much of population genetics is based on the diffusion limit of the Wright–Fisher model, which assumes
a fixed population size. This assumption is violated in most natural populations, particularly for microbes.
Here we study a more realistic model that decouples birth and death events and allows for a stochasticallyvarying population size. Under this model, classical quantities such as the probability of and time before
fixation of a mutant allele can differ dramatically from their Wright–Fisher expectations. Moreover,inferences about natural selection based on Wright–Fisher assumptions can yield erroneous and evencontradictory conclusions: at small population densities one allele will appear superior, whereas at large
densities the other allele will dominate. Consequently, competition assays in laboratory conditions may
not reflect the outcome of long-term evolution in the field. These results highlight the importance ofincorporating demographic stochasticity into basic models of population genetics.

M

ATHEMATICAL descriptions of allele frequencies are typically built upon the Wright–Fisher
model (Wright 1931; Fisher 1958) or, more accurately, its diffusion limit (Kimura 1962; Ewens 2004).
This model forms the basis of Kimura’s work on fixation
probabilities (Kimura 1955), Ewens’ sampling formula(Ewens 1972; Lessard 2007), Kingman’s coalescent(Kingman 1982), tests of neutrality (Hudson et al.
1987; Tajima 1989; McDonald and Kreitman 1991;
Fu and Li 1992; Fay and Wu 2000), and techniques for
inferring mutation rates and selection pressures(Sawyer and Hartl 1992; Yang and Bielawski
2000; Bustamante et al. 2001).
The Wright–Fisher model has been generalized toaccount for a variety of complications, such as nonrandom mating, migration, and multiple loci, among
others (Ewens 2004; Durrett 2009). The standard
diffusion approximation of Kimura (Kimura 1962;
Ewens 2004) and its coalescent (Kingman 1982) areremarkably robust to violations of the model’s underlying assumptions. For example, when the population is stratified, or when the population size varies
rapidly and independently of the genetic compositionof the population, allele frequency dynamics are still
accurately approximated by the Wright–Fisher diffusion, under a suitable change of timescale or choice ofeffective population size (Ewens 1967; Otto andWhitlock 1997; Wakeley 2005, 2009). In fact, mostpopulation-genetic models—including the Moran
process (Moran 1958), Karlin’s conditional branching
1
Corresponding author: Department of Biology, University of Pennsylvania,
Philadelphia, PA 19104. Email: jplotkin@sas.upenn.edu

Genetics 185: 1345–1354 (August 2010)

process (Karlin and McGregor 1964), and someCannings processes (Cannings 1974; Ewens 2004)—
share the same diffusion limit as the Wright–Fisher
model (Möhle 2001). As a result, Kimura’s diffusion
approximation has had an enormous impact on thedevelopment of theoretical and applied population
genetics.
Despite its robustness, the Wright–Fisher diffusion is
not appropriate in all circumstances. Many natural
populations experience substantial stochastic variation.
However, with few exceptions (e.g., Kaj and Krone 2003;
Lambert 2005, 2006; Champagnat and Lambert2007), models typically assume fixed or deterministically varying population number (e.g., Ewens 1967;
Kimura and Ohta 1974; Donnelly 1986; Griffiths
and Tavaré 1994; Otto and Whitlock 1997).
Here, we consider an alternative approach, inspiredby Moran’s model (Moran 1958) and the Gause–Lotka–
Volterra model (Lotka 1925; Volterra 1926; Gause
1934): individuals give birth and die with rates that varywith the total population number. Populations are kept
finite by density-dependent factors (e.g., resource limitation) such as those that have been empirically verified
in microbial populations (Gause 1934; Vandermeer
1969; Pascual and Kareiva 1996). This approachproduces stochastic variation in the population size
because birth events are not immediately followed bydeath events. In this study, we formulate the simplest
population process that incorporates such demographic stochasticity. We show that our model admits adiffusion approximation that is qualitatively different
from the standard Kimura diffusion and that exhibits
novel behavior. In particular, we find that different life-

1346

T. L. Parsons, C. Quince and J. B. Plotkinhistory strategies employed by types of equal expected
lifetime reproductive output have markedly different
consequences for long-term survival. In this regard, ourresults complement Gillespie’s seminal work on fecundity variance and bet hedging (Gillespie 1974, 1975,1977) (see also Taylor 2009 for the corresponding
genealogical process.) Whereas Gillespie consideredvariance in individual reproductive output in a modelwith nonoverlapping generations, we consider allocation to maintenance vs. reproductive output as a formof temporal bet hedging. A similar trade-off is considered in Shpak (2007), via an explicit age-structured
model of fecundity variance where different genotypes
are expressed as different annual survival probabilities
and fecundities. In contrast with these works ofGillespie and Shpak, our model does not assume afixed population size, and it therefore allows us toconsider to consider novel forms of life-history tradeoffs and to contrast small, growing populations withpopulations near equilibrium. By comparing such
populations, we find that bet-hedging mutants can
exhibit an apparent contradiction in their short- vs.
long-term behavior: types that grow rapidly at low
population density are less likely to fix than types thatgrow slowly at low population density.
THE GENERALIZED MORAN MODEL

We consider a haploid population composed of K
allelic types. The number of each type changes as theresult of stochastic birth and death events in continuous
time. Each type i is characterized by its intrinsic birthrate, bi, and death rate, di. Instead of a deterministic
population size, the total population number is stochastic and governed by density dependence in the death
rates. We further allow offspring from type i to mutate totype j at rate mij; we assume mutations are rare; i.e.,
mij ¼ uij =N . We call this model the generalized Moran
process, and we represent it by a Markov chain in which
Xi, the number of individuals of type i, has transition
rates
0
Xi /Xi 1 1 @1 X

1
mij Abi Xi 1

j6¼i

di 1 1mji bj Xj

ð1Þ

j6¼i

PK
Xi /Xi  1

X

j¼1 Xj

N

!
Xi :

ð2Þ

The parameter N determines the typical population
size, as discussed below. This birth–death process isperhaps the simplest possible model of types competingin a population whose size is regulated internally.
Quasi-neutrality: We focus on the case when the ratioof the intrinsic death rate to the birth rate is the same foreach allele, so thatdi
[ a , 1:bi

ð3ÞAlthough one type may have an elevated birth rate, this iscompensated by a correspondingly elevated death rate.
We call this case ‘‘quasi-neutral’’: were the total population
size to be held fixed at Nfix by some exogenous mechanism, all individuals would have the same expected
lifetime reproductive output [i.e., reproductive value(Fisher 1958; Lande 1982; Charnov 1986)] and equalvariance in their lifetime reproductive successes, equal to1
að1 1 Nfix =N Þ

and

1
1
1
;
að1 1 Nfix =N Þ a2 ð1 1 Nfix =N Þ2

respectively [more generally, under this assumption, themoment-generating function for the number of offspring is M ðtÞ ¼ 1 1 ðe t  1Þ=ðað1 1 Nfix =N Þ 1 e t  1Þ,
which is type-independent]. We emphasize that in ourmodel, the total population size varies stochastically,
with the consequence that the different types are not ingeneral exchangeable (Cannings 1974).

ANALYTICAL APPROACH

We analyze the generalized Moran process by considering three phases, according to the size of Xi(t). In thefirst phase, describing the behavior at small initial
numbers, we approximate the process by a linear
birth–death process with Malthusian growth rate bi 
di. The second phase describes the dynamics at intermediate numbers, for which the abundance of each
type, rescaled by N, is asymptotically equivalent to a
deterministic system of Gause–Lotka–Volterra equations for logistic growth. The third phase describes thelong-term stochastic behavior of the process nearcarrying capacity; the relative abundance of each type,when time is measured in units proportional to N,
converges to a diffusion process as N approaches
infinity. The first two limits are valid for all values of bi
and di, while the third exists only when the quasi-neutral
condition (3) holds to Oð1=N Þ (for clarity of exposition,
we restrict our attention to the case when it holds
exactly). These three phases are illustrated in Figure 1
and discussed in greater detail below.
Birth–deathapproximation (phase one): For Xi ðtÞ ¼
p
ﬃﬃﬃﬃﬃ
oð N Þ and N large, density-dependent and mutationalterms in (1) and (2) are o(1) and may be neglected, so
Xi(t) may be approximated by a birth–death processwith birth rate biXi and death rate diXi, for which
E½Xi ðtÞ ¼ e ðbi di Þt :ð4Þ

A coupling argument (Andersson and Djehiche 1998;
Lindvall 2002) may be used to prove convergence
almost surely of the generalized Moran process to thisbirth–death process. Note that this approximation fails

Consequences of Demographic Stochasticity

1347

Figure 1.—Illustrative
sample paths of the generalized Moran process. The
three distinct phases described by our analytic approximations are indicated(seeanalyticalapproach).
Dashed lines show solutions to Equations 4 and 5.

pﬃﬃﬃﬃﬃ
when the total population size is $Oð N Þ, as thedensity-dependent componentof mortality, which isP
proportional to Xi ðtÞð Kj¼1 Xj ðtÞ=N Þ, may be $O(1) andwill not vanish with increasing N.
Deterministic approximation (phase two): For Xi(t) ¼
O(N ), we letdef1
Xi ðtÞ;N /‘ NYi ðtÞ ¼ lim

where the limit is in probability. Then, Yi(t) satisfies thedeterministic system of ordinary differential equations,!
KX_
Y i ¼ bi Y i  d i 1 1
Yj Yi
ð5Þ

P(t) is a diffusion process whose probability densityf(t, p) obeys a forward Kolmogorov equation
K
1
X
@f
@
¼
@t
@p
i
i¼1

"
ui bi pi 1KX

uji bj pj

j¼1

! #
K
X
1
bi pi
1 PK
ðb  bi Þbk pk f
2 ð k¼1 bk pk Þ2 k¼1 k


 
1 K
1
X
bj pj
1 KX
@2
P
bi pi dij  K1
f :2 i¼1 j¼1 @pi @pj
k¼1 bk pk
ð7Þ

j¼1(see Ethier and Kurtz 1986, Chap. 11, Theorem 2.1).
We observe that the rates of mutation are Oð1=N Þ andvanish in this limit; i.e., mutational effects can beignored in the transient dynamics.
When one type has a strictly greater expected reproductive output, i.e.,
djdi.
b i bj

and

dj
di

¼ Oð1ÞbibjPHere, uij ¼ N*mij and ui ¼ Kj¼1 uij . The generator for
P(t) is obtained from the generator of X(t) by applying
Itô’s formula (Protter 2004) to


1
N*
PK
t :
X
2
i¼1 Xi ððN *=2ÞtÞ

ð6Þ

for all j 6¼ i, that type rapidly sweeps to fixation for N large(Parsons and Quince 2007). However, in the quasineutral
P case Equation 5 has equilibria at all points suchthat Ki¼1 yi ¼ ð1=aÞ  1, so that, in the infinite-population
limit, we would expect coexistence of all types.
Diffusion approximation (phase three): On a longer
timescale, we consider the relative abundances of eachtype for large N,


1
N*
t ;
PðtÞ ¼ lim PK
X
*
N /‘
2
i¼1 Xi ððN =2ÞtÞ
where N * ¼ N ð1=a  1Þ is the shared carrying capacity.

Convergence to a diffusion follows by arguments similar
to those in Ethier and Kurtz (1986, Chap. 11), as
described in an associated preprint (Parsons 2010).
We also refer to Durrett and Popovic (2009) for a similar argument, in which the initial process is a diffusion
rather than a Markov chain; see also Katzenberger(1991). The diffusion equation above is used to obtain
the analytical results summarized in Table 1 below. These
analytical results are then illustrated in Figures 3 and4, below, along with comparisons to exact numericalsimulations.
It is important to note that when bi ¼ bj ¼ 1 for all i, j(i.e., when all types are strictly neutral), this diffusion
equation above simplifies to the standard diffusion for a
Wright–Fisher process:1348

T. L. Parsons, C. Quince and J. B. Plotkin
K
X
@f
@
¼
@t
@p
i
i¼1

1

"
ui pk 1

KX! #uji pj f

j¼1

K X
K
X

1
@2
½pi ðdij  pj Þf :
2 i¼1 j¼1 @pi @pj

The first and final terms of (7) correspond to
mutation and genetic drift, and they agree with thestandard diffusion equation of a Wright–Fisher process.
In the quasi-neutral regime, however, for each type ithere is an additional advective term in (7),KX
1bp
PK i iðb  bi Þbk pk ;2
2 ð k¼1 bk pk Þ k¼1 k

proportional to the variation of that type. This term
results in an effective selection coefficient favoring thetype with lowest birth and death rate, when the
population is near carrying capacity.
Results for two types: When K ¼ 2, P(t) is a diffusion
on [0, 1], and we need consider only the proportion oftype one, which we denote p.
In the absence of mutation (uij [ 0 for all i, j ) the
boundaries p ¼ 0 and p ¼ 1 are absorbing. Starting frominitial frequency p, we study the mean first absorption
time, T(p), the probability of fixation of type one, p(p),and the expected time to fixation of type one (i.e., theabsorption time conditioned on type one fixing), T1(p).
These quantities satisfy ordinary differential equations related to the forward Kolmogorov Equation (7)(Ewens 2004):
1 b1 b2 ðb2  b1 Þ
dT
pð1  pÞ2 b1 p 1 b2 ð1  pÞ
dp
1

1
b1 b2
d 2T
pð1  pÞ 2 ¼ 1;
2 b1 p 1 b2 ð1  pÞdp
T ð0Þ ¼ T ð1Þ ¼ 0;

and



 
1 b1 b2 ðb2  b1 Þ
pð1  pÞ f
u12 b1 p 1 u21 b2 ð1  pÞ 1
2 b1 p 1 b2 ð1  pÞ


1 d2
b1 b2
1
pð1

pÞf
¼ 0:
ð8Þ
2 dp 2 b1 p 1 b2 ð1  pÞ
d
dp



(Ewens 2004). For finite N, the generalized Moranmodel has only one absorbing state, corresponding to
extinction of all types. However, the time to global
extinction is exponential in N (and thus tends to infinity
as N / ‘). Prior to extinction, the generalized Moranmodel approaches a quasi-stationary distribution forwhich the relative frequencies are asymptotically distributed according to the solution of (8).
Expressions for p(p), T(p), and T1(p) obtained bysolving the differential equations above are summarizedin Table 1. See, e.g., Ewens (2004) for methods ofsolution.
NUMERICAL METHODS

The figures below show numerical results for the
fixation probabilities and times, along with the analytical approximations derived in the previous section. The
numerical results were obtained by solving, via anumerical library for sparse matrices (Davis 2004),
the matrix equations for the fixation probabilities andtimes as formulated from the transition rates (1) and (2)(Stewart 1994). The number of possible states isunbounded so it was necessary to apply a truncation T
on the maximum population number of both types. It
was found that as the truncation point was moved
beyond the carrying capacity N *, numerical results for
fixation times or probabilities starting from statesaround or below this point converged to within machine
precision by T ¼ 1.5N *; consequently this value was used.

RESULTS1 b1 b2 ðb2  b1 Þ
dp
pð1  pÞ2 b1 p 1 b2 ð1  pÞdp
1

1
b 1 b2d 2p
pð1  pÞ 2 ¼ 0;
2 b1 p 1 b2 ð1  pÞ
dp
pð0Þ ¼ 0;
pð1Þ ¼ 1;

and
1 b1 b2 ðb2  b1 Þdpð1  pÞ ½pðpÞT1 ðpÞ
2 b1 p 1 b2 ð1  pÞ
dp
1

In the presence of mutation (uij . 0), the boundary is
no longer absorbing and P(t) has a stationary distribution, which can be expressed as the solution to

1
b1 b2
d2
pð1  pÞ 2 ½pðpÞT1 ðpÞ ¼ pðpÞ;2 b1 p 1 b2 ð1  pÞ
dp
pð0ÞT1 ð0Þ ¼ pð1ÞT1 ð1Þ ¼ 0:

Figure 2 shows an illustrative sample path of thegeneralized Moran process, starting from a small
population, obtained by Monte Carlo simulation. Initially, the number of each type grows exponentially, andthe type with higher birth and death rates grows faster.
Over the long term, however, the total population size
stabilizes near its carrying capacity, and the type with thelower birth and death rates dominates, because it is lessprone to fluctuation. Thus, it appears that selection initially favors type one and subsequently favors type two—a
phenomenon that we explore in greater detail below.
We have developed a rigorous analytical frameworkfor studying the dynamics of the generalized Moranmodel (see analytical approach). As our analysis
shows, the behavior illustrated in Figure 2 arises because

Consequences of Demographic Stochasticity

1349

TABLE 1
A comparison of the neutral Wright–Fisher model and the quasi-neutral generalized Moran model
Wright–Fisher

Generalized Moran b1
p 1 bb2 1
b2 pð1  pÞ
1

Fixationprobability,p(p)pAbsorptiontimeN*½p ln p 1 ð1  p Þ ln ð1  p Þ"N*
ðb1 ð1 1 pÞ 1 b2 ð1  pÞÞð1  pÞlnð1  pÞ
b1 1 b2

Fixationtime

 Np* ð1  pÞlnð1  pÞðb  b2 Þ2
pð1  pÞ1 ðb1 p 1 b2 ð2  pÞÞp ln p 1 1
b1 1 b2"
1N*
ðb1 ð1 1 pÞ 1 b2 ð1  pÞÞð1  pÞlnð1  pÞpðpÞ b1 1 b2

Stationary
distribution

p u21 1 ð1  p Þu12 1#

1 b1  b2
pð1  pÞð2ðb21  b22 Þp 1 ðb1 1 2b2 Þð3b1 1 b2 ÞÞ1
3 ðb1 1 b2 Þ2#

p u21 ðb2 =b1 Þ1 ð1  pÞu12 ðb1 =b2 Þ1 e ðb1 b2 Þðu12 =b1 1u21 =b2 Þp

Our analytical, asymptotic results for the generalized Moran model are summarized, which are derived in the main text. Absorption and fixation times are given in generations. For the generalized Moran model the generation time is the average time toreproduction, 12 ð1=b1 1 1=b2 Þ, and the carrying capacity is N * ¼ N ð1=a  1Þ. p is the proportion of type one in a population starting from carrying capacity.

the population size is not deterministic. The total
population size will quickly approach and thereafter
fluctuate around the common carrying capacity
N * ¼ ð1  ð1=aÞÞN . Given, as we observed above, thetwo types have equal expected lifetime output at fixed
population size, one might therefore suspect that theclassical quantities, such as the probability and meantime before fixation of a mutant allele, would be givenby their neutral Wright–Fisher expectation with aneffective population size Ne ¼ N */2. This is false,
however. When two quasi-neutral types differ substantially in their vital rates (i.e., by order one), then there isno effective population size, including the harmonicmean (Ewens 1967; Otto and Whitlock 1997), forwhich the behavior of the stochastic population is
captured by a Wright–Fisher diffusion—even if the
population is started from carrying capacity (see analytical results in Table 1).
Under the Wright–Fisher model without mutation, aneutral allele will fix with probability given by its initialfrequency. In the generalized Moran model, however,the fixation probability differs qualitatively from the
Wright–Fisher expectation (Table 1). Starting from a
population at carrying capacity, the fixation probability of a mutant has a quadratic dependence on the
mutant’s initial frequency, and it favors fixation of thetype with low birth and death rates (Figure 3a).
Differences in generation time induce a form ofapparent selection favoring one type over another—
despite the fact that the two types would appear to bestrictly neutral variants under a Wright–Fisher modelwith Ne ¼ N */2.
A similar result holds for the mean time before amutant allele fixes or goes extinct. The behavior of thegeneralized Moran process is skewed relative to theneutral Wright–Fisher prediction: the mean absorption time is longer than the classical expectation when
the type with higher birth rate forms a larger initial
proportion of the population (Figure 3b). The conditional mean time before fixation of an allele is also
skewed: the type with low vital rates fixes more quickly
than a type with high vital rates. Again, the generalized
Moran model exhibits a form of selection that favors one
type over another, despite the fact that the two types
would have the same expected lifetime reproductive
output were the population size held fixed at N *.
When mutations are permitted between the two types,
then the population, conditioned on nonextinction,
will reach a stationary state determined by a balance ofmutation and drift. For the neutral Kimura diffusion
with effective population size Ne ¼ ðN *=2Þ, when mutation rates are symmetric (i.e., m12 ¼ m21 ¼ m), thestationary relative frequency, p, of type one follows thewell-known beta distributionf ðpÞ ¼ p u1 ð1  pÞu1(Ewens 2004), where u ¼ 2Nem. For a quasi-neutral
stochastic population, by contrast, the stationary distribution of type one is given by the solution to (8):1350

T. L. Parsons, C. Quince and J. B. Plotkin

Figure 2.—An illustrative stochastic realization of the generalized Moran process, starting from a population composedof 10 individuals of each type. Initially, the type with the highbirth and death rates grows faster (top), but over the longterm the type with low birth and death rates dominates (bottom). The two types have the same expected lifetime reproductive output, and they would behave neutrally in astandard Wright–Fisher model with fixed population size.
The two panels show output from the same simulation, withdifferent scales on the time axis. Parameters values: b1 ¼ 2,
b2 ¼ 1, a ¼ 0.5, N ¼ 1000, and m12 ¼ m21 ¼ 0.01.

p uðb2 =b1 Þ1 ð1  pÞuðb1 =b2 Þ1 e uðb1 b2 Þð1=b1 11=b2 Þp :

ð9Þ

As we demonstrate below, the stationary allele frequency
distribution of the generalized Moran process (Equation 9) is very similar to that of a Wright–Fisher modelwith selection favoring the type with low birth rate.
All of the analytical results summarized in Table 1 areasymptotic in N. These asymptotic approximations arehighly accurate for N ¼ 1000, as demonstrated in
Figures 3 and 4, and therefore also for N . 1000.

PRACTICAL IMPLICATIONS

The results above have important consequences forpractitioners who wish to make inferences about natural
selection on the basis of samples from wild or laboratory
populations. When life-history traits are under genetic
control, inference based on standard assumptions (i.e.,
Kimura’s diffusion) may yield erroneous and evencontradictory conclusions: at small population densities

Figure 3.—(A) The probability that type one fixes, as afunction of its initial proportion, for the Wright–Fisher model(blue line) and the generalized Moran model (exact numerical results, solid red line; analytical results from Table 1,dashed black line) starting at carrying capacity. (B) The absorption time of type one under the Wright–Fisher model(blue line) and the generalized Moran model (exact numerical results, solid red line; analytical results from Table 1,dashed black line). Parameter values: b1 ¼ 2, b2 ¼ 1, a ¼
0.5, and N ¼ 1000.

one allele will appear superior, whereas at large densities
the other allele will dominate. We illustrate this point
using two examples of the types of inferences typicallydrawn from population data.
In the first example, we analyzed polymorphism‘‘data’’ generated by simulating the generalized Moran
processes of competing, quasi-neutral types. We sampled 100 individuals from a large population nearcarrying capacity, and we recorded the polymorphism
frequency spectrum across 1000 independent sites(Figure 4). The observed spectrum, which describes
the number of sites at which the type-one allele
segregates at different frequencies, is typical of the polymorphism data observed in field studies of microbes
and higher eukaryotes (e.g., Hartl et al. 1994; Akashi
and Schaeffer 1997). If we analyze the simulated
polymorphism data using standard inference tech-Consequences of Demographic Stochasticity

Figure 4.—The stationary distribution of allele frequencies under reversible mutation. Bars show results from Monte
Carlo simulations of the quasi-neutral generalized Moran process (N ¼ 1000, a ¼ 0.5, b1 ¼ 1.6, b2 ¼ 2, and m ¼ 0.0005), in
frequency bins of size 10%. Dots show our analytic approximation to the quasi-neutral process, taken from the equation in
Table 1. The analytical approximation closely matches the numerical simulation. In addition, the quasi-neutral frequency
spectrum is very similar to the spectrum produced by a Wright–
Fisher model with selection (triangles). The parameters of abest-fit Wright–Fisher model (û ¼ 0:512; N̂e s ¼ 0:252) wereobtained by applying standard maximum-likelihood inference techniques (Sawyer and Hartl 1992; Bustamante
et al. 2001; Desai and Plotkin 2008) to the simulated data,
and they would indicate selection favoring the type withlow birth rate.

niques based on the Wright–Fisher diffusion (Sawyerand Hartl 1992; Bustamante et al. 2001; Desai and
Plotkin 2008), we would infer directional selection
favoring the type with low birth rate (Figure 4).
This result is important because alternative approaches to inferring selection pressures from data
would reach the opposite conclusion. Selection coefficients are often quantified by competing two types in anutrient-rich laboratory environment and measuring
the difference in their initial exponential growth rates(Lenski et al. 1991; Lenski and Elena 2003). Under thegeneralized Moran model starting from a small population, type one will grow at exponential rate b1  d1 ¼
b1(1  a) and type two at exponential rate b2  d2 ¼
b2(1  a), yielding an apparent selective advantage tothe type with high birth rate. Thus, the selection coefficient inferred from a laboratory competition experiment will have the opposite sign to the selection
coefficient inferred from the stationary polymorphism
frequency spectrum, as would be sampled in the field.
To take a textbook example, in a 12-hr competition
experiment from small initial numbers the frequency ofEscherichia coli strain gnd(RM77C) increased from 0.455
to 0.898, relative to strain gnd(RM43A) (Hartl and

1351

Figure 5.—Bottleneck frequency influences fixation probability under the quasi-neutral generalized Moran model (redcircles), but not under the neutral Wright–Fisher model (bluedashed line). Simulated populations were initiated at carrying
capacity with an equal number of individuals of type one andtype two. The populations were then periodically bottlenecked to 10 individuals, by binomial sampling. When bottlenecks are frequent, the type with larger birth rate is more
likely to fix and conversely when bottlenecks are infrequent.
Parameter values: N ¼ 1000, a ¼ 0.5, b1 ¼ 0.4, and b2 ¼ 2.

Dykhuizen 1981). The standard interpretation of thesedata suggests that strain RM77C is selectively advantageous (Nes ¼ 0.067) and that it will eventually fix with aprobability greater than its frequency at 12 hr (Hartl
and Clark 2006). But these data are equally compatible
with the quasi-neutral generalized Moran model, which
predicts that RM77C will be disadvantageous when the
population eventually reaches carrying capacity andthat it will fix with a probability less than its frequency
at 12 hr. Thus, the interpretation of such a competition
experiment depends on the trade-offs that occur between birth rate and other life-history traits and also on
assumptions about the underlying demographic process.
Our results also have implications for serial transfer
experiments used to study microbial evolution (e.g.,
Lenski et al. 1991; Bull et al. 1997; Holder and Bull
2001). In Figure 5 we consider the generalized Moran
process with the addition of periodic bottlenecks in the
population size, representative of periodic transfers.
Fixation probabilities for the process with bottlenecks
were obtained via Monte Carlo simulations. In this case,the probability that one quasi-neutral type will fix as
opposed to another depends on the bottleneck frequency: when bottlenecks are infrequent, the type withlow birth rate is more likely to fix, whereas the opposite
occurs when bottlenecks are frequent. This result
implies that the choice of serial transfer protocol induces
a form of selection for one or another life-history
phenotype, according to the transfer frequency—a phenomenon that should be considered when designing andinterpreting laboratory studies of microbial evolution.

1352

T. L. Parsons, C. Quince and J. B. Plotkin

Previous work on life history has reached similar
conclusions about fixation probabilities. For example,
Alexander and Wahl (2008) found that two different
beneficial mutations with the same expected lifetime
reproductive output may have different fixation probabilities relative to a disadvantageous wild type. However, such analyses were based on a branching process
approximation that applies only in the limit of strong
selection, when the fixation of the mutant is ensured
once it escapes genetic drift (Wahl and Gerrish 2001;
Wahl and Dehaan 2004; Alexander and Wahl 2008).
By contrast, we have demonstrated that an invading
mutant may be favored or disfavored compared to wild
type, depending on the population size when it isintroduced or the frequency of bottlenecks.
DISCUSSION

Here we have studied a simple stochastic model ofbirths and deaths in continuous time, which is more
realistic than the Wright–Fisher model in many circumstances. We have derived analytic expressions for fixation probabilities, absorption times, and steady-statepolymorphism distributions. We find that when the
population size is near carrying capacity, one type ispreferred, and when far from carrying capacity, theother type is preferred—despite the fact that both types
have the same expected lifetime reproductive output.
These results have implications for our understandingof molecular evolution and for the interpretation ofevolutionary studies in the field and in the laboratory.
The behavior of the generalized Moran process in thequasi-neutral regime can be understood as a form of
r-vs.-K selection (Pianka 1972; MacArthur and Wilson
2001). Quasi-neutrality enforces a life-history trade-off
in which an increased birth rate is compensated by anincreased death rate. The type with the lower mortality
rate is favored when the population size is large (i.e.,near carrying capacity), because such a type experiences
less dramatic fluctuations and is less likely to go extinct.
As a result, fixation probabilities, absorption times, andthe stationary polymorphism frequency spectrum favor
the type with low birth rate. In contrast, when the total
population size is small, so that density-dependent regulation is weak, rapid growth is preferable and the typewith higher birth rate is favored. Thus, in a model withstochastic population size there is no one-dimensional‘‘selection coefficient’’ that summarizes the dynamics
of competing types (see also Lambert 2005, 2006;
Champagnat and Lambert 2007).
A trade-off between mortality and fecundity may beinterpreted as a form of evolutionary bet hedging(Slatkin 1974), where high-fecundity strategies produce greater variability in lifetime reproductive success,but they are selectively disadvantageous in saturated
environments. Similar results have previously been
observed by Gillespie and others, who studied models

in which individuals produce variable numbers ofoffspring (Gillespie 1974, 1975, 1977; Proulx 2000;
Shpak 2007). Like the generalized Moran model, suchmodels are not equivalent to a Wright–Fisher diffusion,
for any constant value of the effective population size.
Nonetheless, our results differ from those of Gillespie,
who assumed that the population size is held constant bysome exogenous mechanism. We have shown that theseeffects arise in populations that are intrinsically regulated by density dependence. Indeed, differences inindividual fecundity variance appear only when the total
population size is allowed to vary. Moreover, when the
population size is stochastic, one type will be favored at
low density and another type at high density. Combinedwith exogenous disturbances, these effects can drastically change the probability of fixation (e.g., Figure 5).
Quasi-neutrality is a novel phenomenon that arises
when the total population number is allowed to varystochastically, in a hybrid regime between strong andweak selection. When one type has a significantly greater
lifetime reproductive success [i.e., bi =di . bi =di 1 Oð1Þ],then that type will sweep with probability approaching
one as the carrying capacity becomes large. By contrast,
when the birth and death rates differ by a small amount[b1  b2 ¼ Oð1=N Þ and d1  d2 ¼ Oð1=N Þ], we recover
the classical theory of weak selection. Quasi-neutrality
arises in an intermediate regime: when two or more types
have overall reproductive success that differs by at mostOð1=N Þ, achieved through birth and death rates thatdiffer by O(1).
Since quasi-neutrality occurs in a limited regime ofparameters, the question arises as to whether thisregime is found in natural populations. After all, any
mutations that confer a significant improvement inoverall reproductive success will quickly displace neutralor quasi-neutral segregating types. Indeed, such arguments have been raised against the biological relevanceof Gillespie’s model (e.g., Hopper et al. 2003 andreferences therein), which operates in a similar parameter regime. Nevertheless, there are both theoreticaland empirical reasons why quasi-neutral types might
regularly segregate in populations—arguments analogous to those that support widespread neutral variation.
In a fixed environment, mutations that confer strong
selective advantages will fix quickly and be depleted(Lenski and Travisano 1994); and so eventually mostsegregating variation will be neutral, weakly selected, or
quasi-neutral. Once strongly selected mutations havebeen exhausted, any mutations that allow an organism
to vary its allocation to reproduction vs. survivorship will
begin to segregate. Thus, quasi-neutral types will most
likely segregate in depleted or subdivided environments
where the carrying capacity of a given patch is small(Shpak 2005) or in populations where most mutations
that convey an absolute fitness advantage have already
fixed.

Consequences of Demographic Stochasticity

We have focused on the quasi-neutral case, bi =di [a tohighlight the contrast between the generalized Moran
process and the standard, neutral Wright–Fisher diffusion. Nevertheless, our analysis extends to competingtypes with different expected reproductive outputs—i.e.,values of bi =di that differ by Oð1=N Þ (Parsons 2010). In
other words, the r-vs.-K dynamic we have observed also
operates in the presence of traditional forms of weak
selection.
Fecundity/mortality trade-offs of the broad type wehave considered here may be ubiquitous in microorganisms, arising from multiple alternative mechanisms. Such trade-offs have been observed empiricallyin E. coli subject to prolonged starvation: lower intrinsic
mortality rates allowed adapted strains to persist in
minimal salts media, but such strains proved to beinferior competitors in fresh medium (Vasi and Lenski1999). Trade-offs between phage resistance and growthrate—arising because resistance involves the loss ofreceptors on the cell surface that are both exploitedby the phage to gain entry and involved in the uptake ofnutrients—are also well documented (Bohannan andLenski 2000; Bohannan et al. 2002). Similarly, antibiotic resistance often causes a reduction in doubling time,
because of the costs of actively transporting the antibiotic out of the cell through efflux pumps, e.g., in fluoroquinolone resistance in Salmonella enterica (Giraud et al.
2003). The exact shape of these trade-offs may differ
from the linear one considered here (Jessup and
Bohannan 2008), and our population dynamics—with
demographic rates that depend linearly on population
number—will apply only as an approximation. Nonetheless, the potential exists for conditions analogous toquasi-neutrality, such that the corresponding deterministic model of competition predicts stable coexistence.
Indeed, empirical studies of E. coli strains competing ina chemostat (Hansen and Hubbell 1980) under justthese conditions have produced time series in which
one type predominates at low density but is then lessabundant in saturated media, much as our model
predicts (Figure 2). Consequently the phenomena wehave investigated, i.e., stochastic contributions to thefitness of competing types, may be important in shaping
patterns of microbial genetic diversity.
We thank Peter Abrams and Warren Ewens for their insightful input;
Jeremy Draghi, Sergey Kryazhimsky, Jesse Taylor, and one anonymous
reviewer for comments that led to substantial improvements in the
manuscript; and Rick Durrett for bringing Durrett and Popovic(2009) to our attention. C.Q. is supported by an Engineering andPhysical Sciences Research Council Career Acceleration Fellowship(EP/H003851/1). J.B.P. acknowledges support from the Burroughs
Wellcome Fund, the David and Lucile Packard Foundation, the James
S. McDonnell Foundation, and the Alfred P. Sloan Foundation.

LITERATURE CITED
Akashi, H., and S. W. Schaeffer, 1997 Natural selection and thefrequency distributions of ‘‘silent’’ DNA polymorphism in Drosophila. Genetics 146: 295–307.

1353

Alexander, H. K., and L. M. Wahl, 2008 Fixation probabilities depend on life history: fecundity, generation time and survival in aburst-death model. Evolution 62: 1600–1609.
Andersson, H., and B. Djehiche, 1998 A threshold limit theorem
for the stochastic logistic epidemic. J. Appl. Probab. 35: 662–670.
Bohannan, B. J. M., and R. E. Lenski, 2000 Linking genetic change
to community evolution: insights from studies of bacteria andbacteriophage. Ecol. Lett. 3: 362–377.
Bohannan, B. J. M., B. Kerr, C. M. Jessup, J. B. Hughes and G. Sandvik,
2002 Trade-offs and coexistence in microbial microcosms. Antonie Van Leeuwenhoek 81: 107–115.
Bull, J. J., M. R. Badgett, H. A. Wichman, J. P. Huelsenbeck, D. M.
Hillis et al., 1997 Exceptional convergent evolution in a virus.
Genetics 147: 1497–1507.
Bustamante, C. D., J. Wakeley, S. Sawyer and D. L. Hartl,
2001 Directional selection and the site-frequency spectrum.
Genetics 159: 1779–1788.
Cannings, C., 1974 The latent roots of certain Markov chains arising in genetics: a new approach, i. haploid models. Adv. Appl.
Probab. 6: 260–290.
Champagnat, N., and A. Lambert, 2007 Evolution of discrete populations and the canonical diffusion of adaptive dynamics. Ann.
Appl. Probab. 17: 102–155.
Charnov, E. L., 1986 Dimensionless numbers and life history evolution: age of maturity versus the adult lifespan. Evol. Ecol. 4:
273–275.
Davis, T. A., 2004 Algorithm 832: UMFPACK, an unsymmetricpattern multifrontal method. ACM Trans. Math. Soft. 30: 196–199.
Desai, M. M., and J. B. Plotkin, 2008 The polymorphism frequency
spectrum of finitely many sites under selection. Genetics 180:
2175–2191.
Donnelly, P., 1986 A genealogical approach to variable-populationsize models in population genetics. J. Appl. Probab. 23: 283–296.
Durrett, R., 2009 Probability Models for DNA Sequence Evolution, Ed.
2. Springer, New York.
Durrett, R., and L. Popovic, 2009 Degenerate diffusions arisingfrom gene duplication models. Ann. Appl. Probab. 19: 15–48.
Ethier, S. N., and T. G. Kurtz, 1986 Markov Processes: Characterization and Convergence. Wiley Interscience, Hoboken, NJ.
Ewens, W. J., 1967 The probability of survival of a new mutant in afluctuating environment. Heredity 22: 438–443.
Ewens, W. J., 1972 The sampling theory of selectively neutral alleles.
Theor. Popul. Biol. 3: 87–112.
Ewens, W. J., 2004 Mathematical Population Genetics, Ed. 2. Springer,
New York.
Fay, J. C., and C.-I. Wu, 2000 Hitchhiking under positive Darwinian
selection. Genetics 155: 1405–1413.
Fisher, R. A., 1958 The Genetical Theory of Natural Selection. Dover,
New York.
Fu, Y. X., and W. H. Li, 1992 Statistical tests of neutrality of mutations. Genetics 133: 693–709.
Gause, G. F., 1934 The Struggle for Existence. Williams & Watkins,
Baltimore.
Gillespie, J. H., 1974 Nautural selection for within-generation variance in offspring number. Genetics 76: 601–606.
Gillespie, J. H., 1975 Natural selection for within-generation variance in offspring number ii. discrete haploid models. Genetics
81: 403–413.
Gillespie, J. H., 1977 Natural selection for variances in offspring
numbers: a new evolutionary principle. Am. Nat. 111: 1010–1014.
Giraud, E., A. Cloeckaert, S. Baucheron, C. Mouline andE. Chaslus-Dancla, 2003 Fitness cost of fluoroquinolone resistance in Salmonella enterica serovar Typhimurium. J. Med.
Microbiol. 52: 697–703.
Griffiths, R. C., and S. Tavaré, 1994 Sampling theory for neutral
alleles in a varying environment. Philos. Trans. R. Soc. Lond. B
344: 403–410.
Hansen, S. R., and S. P. Hubbell, 1980 Single-nutrient microbial
competition: qualitative agreement between experimental andtheoretically forecast outcomes. Science 207: 1491–1493.
Hartl, D. L., and A. G. Clark, 2006 Principles of Population Genetics,
Ed. 4. Sinauer Associates, Sunderland, MA.
Hartl, D. L., and D. E. Dykhuizen, 1981 Potential for selection
among nearly neutral allozymes of 6-phosphogluconate dehydrogenase in Escherichia coli. Proc. Natl. Acad. Sci. USA 78: 6344–6348.

1354

T. L. Parsons, C. Quince and J. B. Plotkin

Hartl, D. L., E. N. Moriyama and S. A. Sawyer, 1994 Selection intensity for codon bias. Genetics 138: 227–234.
Holder, K. K., and J. J. Bull, 2001 Profiles of adaptation in two similar viruses. Genetics 159: 1393–1404.
Hopper, K., J. Rosenheim, T. Prout and S. Oppenheim,
2003 Within-generation bet hedging: A seductive explanation?
Oikos 101: 219–222.
Hudson, R. R., M. Kreitman and M. Aguadé, 1987 A test of neutral
molecular evolution based on nucleotide data. Genetics 116:
153–159.
Jessup, C. M., and B. J. M. Bohannan, 2008 The shape of an ecological trade-off varies with environment. Ecol. Lett. 11: 947–959.
Kaj, I., and S. M. Krone, 2003 The coalescent process in a population with stochastically varying size. J. Appl. Probab. 40: 33–48.
Karlin, S., and J. McGregor, 1964 Direct product branching processes and related Markov chains. Proc. Natl. Acad. Sci. USA 51:
598–602.
Katzenberger, G. S., 1991 Solutions of a stochastic differential
equation forced onto a manifold by a large drift. Ann. Probab.
19: 1587–1628.
Kimura, M., 1955 Solution of a process of random genetic drift witha continuous model. Proc. Natl. Acad. Sci. USA 41: 144–150.
Kimura, M., 1962 On the probability of fixation of mutant genes ina population. Genetics 47: 713–719.
Kimura, M., and T. Ohta, 1974 Probability of gene fixation in an
expanding finite population. Proc. Natl. Acad. Sci. USA 71:
3377–3379.
Kingman, J. F. C., 1982 The coalescent. Stoch. Proc. Appl. 13: 235–248.
Lambert, A., 2005 The branching process with logistic growth. Ann.
Appl. Probab. 69: 419–441.
Lambert, A., 2006 Probability of fixation under weak selection: abranching process unifying approach. Theor. Popul. Biol. 15:
1506–1535.
Lande, R., 1982 A quantitative genetic theory of life history evolution. Ecology 63: 607–615.
Lenski, R., and S. F. Elena, 2003 Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation. Nat.
Rev. Genet. 4: 457–469.
Lenski, R., M. R. Rose, S. Simpson and S. C. Tadler, 1991 Longterm experimental evolution in Escherichia coli i. Am. Nat.
138: 1315–1341.
Lenski, R. E., and M. Travisano, 1994 Dynamics of adaptation anddiversification: a 10,000-generation experiment with bacterial
populations. Proc. Natl. Acad. Sci. USA 91: 6808–6814.
Lessard, S., 2007 An exact sampling formula for the Wright–Fisher
model and a solution to a conjecture about the finite-islandmodel. Genetics 177: 1249–1254.
Lindvall, T., 2002 Lectures on the Coupling Method. Dover, Mineola, NY.
Lotka, A. J., 1925 Elements of Physical Biology. Williams & Watkins,
Baltimore.
MacArthur, R. H., and E. O. Wilson, 2001 The Theory of Island Biogeography. Princeton University Press, Princeton, NJ.
McDonald, J., and M. Kreitman, 1991 Adaptive protein evolution
at the Adh locus in Drosophila. Nature 351: 652–654.
Möhle, M., 2001 Forward and backward diffusion approximations
for haploid exchangeable population models. Stoch. Proc. Appl.
95: 133–149.

Moran, P. A. P., 1958 Random processes in genetics. Proc. Camb.
Philos. Soc. 54: 60–71.
Otto, S. P., and M. C. Whitlock, 1997 The probability of fixationin populations of changing size. Genetics 146: 723–733.
Parsons, T. L., 2010 Limit theorems for competitive density dependent population processes. arXiv, http://arxiv.org/abs/1005.0010.
Parsons, T. L., and C. Quince, 2007 Fixation in haploid populations exhibiting density dependence I: the non-neutral case. Theor. Popul. Biol. 72: 121–135.
Pascual, M. A., and P. Kareiva, 1996 Predicting the outcome ofcompetition using experimental data: maximum likelihood
and Bayesian approaches. Ecology 77: 337–349.
Pianka, E. R., 1972 r and K selection or b and d selection? Am. Nat.
106: 581–588.
Protter, P. E., 2004 Stochastic Integration and Differential Equations,Ed. 2. Springer, Berlin.
Proulx, S. R., 2000 The ESS under spatial variation with applications to sex allocation. Theor. Popul. Biol. 58: 33–47.
Sawyer, S. A., and D. L. Hartl, 1992 Population genetics of polymorphism and divergence. Genetics 132: 1161–1176.
Shpak, M., 2005 Evolution of variance in offspring number: theeffects of population size and migration. Theor. Biosci. 124:
65–85.
Shpak, M., 2007 Selection against demographic stochasticity in agestructured populations. Genetics 177: 2181–2194.
Slatkin, M., 1974 Hedging one’s evolutionary bets. Nature 250:
704–705.
Stewart, W. J., 1994 Introduction to the Numerical Solution of Markov
Chains. Princeton University Press, Princeton, NJ.
Tajima, F., 1989 Statistical method for testing the neutral mutation
hypothesis by DNA polymorphism. Genetics 123: 585–595.
Taylor, J. E., 2009 The genealogical consequences of fecundity
variance polymorphism. Genetics 182: 813–837.
Vandermeer, J. H., 1969 The competitive structure of communities:
an experimental approach with protozoa. Ecology 50: 362–371.
Vasi, F. K., and R. E. Lenski, 1999 Ecological strategies and fitness
tradeoffs in Escherichia coli mutants adapted to prolonged starvation. J. Genet. 78: 43–49.
Volterra, V., 1926 Fluctuations in the abundance of a species considered mathematically. Nature 118: 558–560.
Wahl, L. M., and C. S. DeHaan, 2004 Fixation probability favors increased fecundity over reduced generation time. Genetics 168:
1009–1018.
Wahl, L. M., and P. J. Gerrish, 2001 Fixation probability in populations with periodic bottlenecks. Evolution 55: 2606–2610.
Wakeley, J., 2005 The limits of theoretical population genetics.
Genetics 169: 1–7.
Wakeley, J., 2009 Coalescent Theory: An Introduction. Roberts & Co.,New York.
Wright, S., 1931 Evolution in Mendelian populations. Genetics 16:
97–159.
Yang, Z., and J. P. Bielawski, 2000 Statistical methods for detecting
molecular adaptation. Trends Ecol. Evol. 15: 496–503.

Communicating editor: L. M. Wahl