Theoretical Population Biology 74 (2008) 302–310

Contents lists available at ScienceDirect

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

Absorption and fixation times for neutral and quasi-neutral populations withdensity dependence
Todd L. Parsons a,∗ , Christopher Quince b , Joshua B. Plotkin aaDepartment of Biology, University of Pennsylvania, Philadelphia, PA 19104, USAb

Department of Civil Engineering, University of Glasgow, Glasgow, G128LT, UK

articleinfo

Article history:
Received 21 April 2008
Available online 18 September 2008
Keywords:
Fixation
Absorption
Density-dependent
Logistic
Stochastic
NeutralAsymptotic
AlleleCompetition
Drift
Moran
Wright–Fishera b s t r a c t
We study a generalisation of Moran’s population-genetic model that incorporates density dependence.
Rather than assuming fixed population size, we allow the number of individuals to vary stochasticallywith the same events that change allele number, according to a logistic growth process with density
dependent mortality. We analyse the expected time to absorption and fixation in the ‘quasi-neutral’ case:both types have the same carrying capacity, achieved through a trade-off of birth and death rates. Suchtypes would be competitively neutral in a classical, fixed-population Wright–Fisher model. Nonetheless,we find that absorption times are skewed compared to the Wright–Fisher model. The absorption time islonger than the Wright–Fisher prediction when the initial proportion of the type with higher birth rate islarge, and shorter when it is small. By contrast, demographic stochasticity has no effect on the fixation orabsorption times of truly neutral alleles in a large population. Our calculations provide the first analytic
results on hitting times in a two-allele model, when the population size varies stochastically.
© 2008 Elsevier Inc. All rights reserved.

1. Introduction
The time to fixation or extinction of a mutant allele in a finite population is a classical problem of population genetics, andone in which there has been sustained interest (Kimura and Ohta,
1969a,b; Littler, 1975; van Herwaarden and van der Wal, 2002;
Ewens, 2004; Taylor et al., 2006). For reasons of mathematicaltractability, previous studies have typically assumed fixed population sizes. In such models with deterministic population size,
birth events are followed immediately by death events. In natural populations, however, birth and death events are independent,
and the total population size varies stochastically around a typicalvalue, often governed by density dependence. Here we study fixation times in a model that allows for stochastic population sizes,
which is therefore more realistic than the standard Wright–Fisher
formulation (Wright, 1931; Fisher, 1958; Ewens, 2004).
We have recently demonstrated (Parsons and Quince, 2007a,b)that when population size varies stochastically, the probability offixation differs qualitatively from that of fixed or deterministic
population size models (Kimura, 1957, 1962; Moran, 1958; Ewens,∗ Corresponding author.
E-mail address: tparsons@sas.upenn.edu (T.L. Parsons).
0040-5809/$ – see front matter © 2008 Elsevier Inc. All rights reserved.
doi:10.1016/j.tpb.2008.09.001

1967; Kimura and Ohta, 1974; Otto and Whitlock, 1997). Similar
results appear in Lambert (2006). In particular, in the case of
two types that have the same carrying capacity, but achieve thiscarrying capacity through different life-history strategies – what
we call the quasi-neutral case – the probability of fixation depends
not only on the initial density of the mutant type, as would be
expected from a fixed-population size model, but also on the initial
population size.
Here, we study the expected time to absorption or fixation
for two quasi-neutral types. We use the framework introduced
in Parsons and Quince (2007a): a two-dimensional Markov chainthat models competition between discrete individuals of two
types (mutant and resident) in continuous time. Competition
is modelled via density dependent mortality, whose strength
varies inversely with ‘system size’ N. The system size reflects theavailability of resources, and it determines the carrying capacityof all types. We use singular perturbation analysis (Grasman
and van Herwaarden, 1999) of a diffusion approximation to
obtain asymptotic expressions for the mean time to absorption
or fixation. These approximations increase in accuracy withincreasing N. Consistent with our previous findings, we findthat for types with different birth rates, the time to fixation
grows with the difference in birth rates. We also determine theexpected absorption and fixation times for an arbitrary number of

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

truly neutral types — all having the same birth and death rates.
In this case, we obtain the classical expressions with the total
population number replaced by the carrying capacity — thereby
demonstrating that under strict neutrality the classical results
remain valid even when the assumption of fixed population size isrelaxed. Lastly, we verify our approximations by comparison withexact numerical computations, and we investigate the accuracyof our approximations at small population sizes over a range ofparameters.
In what follows, we present our model and our asymptotic
results with discussion, with the details of the derivation in the
appendices.
2. A model of allele fixation with density dependent mortality
We consider a population consisting of K types, subject todensity dependent regulation. All types have the same carrying
capacity, which is proportional to the ‘system size’ N. The systemsize represents the pool of available resources; as the system size
is enlarged, carrying capacity of all types increases.
We denoteof individuals of each type by XN =
 the number
NNN
X1 , . . . , XK . Here, Xi is the number of individuals of type i.
We emphasise N, as we will be interested in the asymptotic
behaviour as the system size becomes large. We assume that the
number of the each type changes as a result of stochastic birthand death events in continuous time, for which the transition
rates — probabilities per unit time — depend only on the current

Kstate of the system. This gives a Markov chain on Z≥0
withtransitions only to neighbouring states, i.e., a K -dimensional birthand-death process (van Kampen, 1992), with transition rates asgiven in Table 1.
Here βi and δi represent the birth and intrinsic death rates oftype i, respectively, and the system size N determines the carrying
capacity. δi gives a lower bound for the per capita death rate, which
increases in proportion to the total number of individuals and isinversely proportional to the system size.
Our model naturally generalises the Moran model (Moran,1958) by making births and deaths independent events, thereby
allowing the total population number to vary stochastically withthe same events that cause changes in allele number.
We previously determined fixation probabilities in the caseK = 2 (Parsons and Quince, 2007a,b) for a closely related processwith density dependent fecundity. While the model here may
appear different, due to density dependence in the mortality term,
the arguments and results for fixation times presented here are
identical for both models (unpublished results).
In this study, we assume that all types are quasi-neutral –i.e. that the ratio of birth rates to intrinsic death rates is fixed,δi= α , for all types i. This corresponds to the case when all
βi
types have the same value for Kimura’s selection coefficient, s.
Thus, in a Wright–Fisher model assuming fixed population size,
all such types would be competitively neutral. However, when the
population size is allowed to vary stochastically, the population
dynamics correspond to a form of r vs. K selection (MacArthur
and Wilson, 1967; Pianka, 1972): types with higher birth rates
are favoured at small initial population sizes and types withlower mortality are favoured as the population size becomes large(Parsons and Quince, 2007b).
As N approaches infinity, N −1 XN approaches a deterministic
dynamical system,
Ẋi = βi Xi − δi Xi

1+

K
X

303

Table 1
Density dependent mortality
Transition

Rate

XiN → XiN + 1

βi XiN 

δi XiN 1 +XiN → XiN − 1

PKN
j=1 X jNsystem, corresponding to coexistence of all types. When N is finite,the total number of individuals will rapidly approach α1 − 1 N,
and will remain near this value for a long time, before demographic
stochasticity leads to the extinction of one of the types. We thus
refer to N ∗ = α1 − 1 N as the carrying capacity of the system. We
are primarily concerned with the case when α is much smaller than1, so that the population spends most of its time near the carrying
capacity.

3. Results
Strictly speaking, the only absorbing state of the model above
is extinction of all types. However, the total population number,
obtained by summing all types, obeys a logistic process, for which
the time before such extinction is exponentially large in the systemsize, N (Newman et al., 2004). Long before the entire population
goes extinct, the population will become monomorphic for one of
the K types.
We are concerned here with the the first time at which the
population is monomorphic, which we call the absorption time.
Focusing on one particular type, i, we define the fixation timeas the absorption time conditioned on monomorphism of type i.
We will obtain asymptotic expressions for the expected time toabsorption and the expected (conditional) time to fixation for two
quasi-neutral types, and for an arbitrary number of strictly neutral
types. These expressions become exact as N approaches infinity,
and, as we confirm with numerical analysis, they are already very
accurate for N of the order of hundreds. We present our results
below, with details of the derivation in the appendices.
Asymptotically, both hitting times depend only on the initialstate of the system. We express these hitting times as a functionof new co-ordinates σi , which parameterise the flows of thedeterministic system Eq. (A.6) (Appendix A.3). If the types arestrictly neutral or if the population size starts at carrying capacity,then σi is simply the initial proportion of type i in the population,
namelyσi =Xi (0)KP.

Xj (0)

j =1

In general, explicit expressions for the σi do not exist, although
they are well-defined and can be calculated numerically for any
parameters and initial population sizes (Appendix A.3).

3.1. Absorption and fixation times for two quasi-neutral types
For two quasi-neutral types, the expected time to firstabsorption of either type is:

!
Xj

.

(1)

j=1

Under the assumption of quasi-neutrality, every point on the linePK1
j=1 xj = α − 1 is a possible equilibrium point of the deterministic

T (σ1 ) = −N∗
2β1 β2(β1 (1 + σ1 ) + β2 (1 − σ1 )) (1 − σ1 ) ln(1 − σ1 )+ (β1 σ1 + β2 (2 − σ1 )) σ1 ln σ1 +
(β1 − β2 )2
σ1 (1 − σ1 ) ,β1 + β2(2)

304

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

while the conditional expected time to fixation for type 1 isT1 (σ1 ) = −1

N∗[(β1 (1 + σ1 )
π1 (σ1 ) 2β1 β2+ β2 (1 − σ1 )) (1 − σ1 ) ln(1 − σ1 )1 β1 − β2+
σ1 (1 − σ1 ) 2(β12 − β22 )σ1
3 (β1 + β2 )2+ (β1 + 2β2 )(3β1 + β2 ))] ,(3)

where π1 is the asymptotic probability that type 1 fixes (Parsons
and Quince, 2007b):

π1 (σ1 ) = σ1 +β2 − β1
σ1 (1 − σ1 ).
β1 + β2(4)

In subsequent sections we compare these hitting times to those
expected under a neutral model. For now, we note the time to
fixation of a single mutant of type 1 invading a resident populationof type 2 is:
1 N ∗ (5β12 + 14β1 β2 + 5β22 )12

β1 β2 (β1 + β2 ),(5)as N gets large. This is obtained by taking the limit σ → 0 in Eq. (3).
This expression is symmetric: a single mutant of type 2 invadinga large resident population of type 1 has the same expected time
to fixation. This result is interesting in light of Eq. (4), whichshows that the type with lower birth rate has higher probabilityof fixation. The symmetric behaviour of the fixation time arises
because we have conditioned on fixation. We discuss this issue
further below.
Our analysis requires that the population size spend a long
period of time near carrying capacity. This condition will be met
provided α is much less than unity. α1 is strictly larger than the
expected lifetime reproductive success of an individual, so as
α approaches unity, the carrying capacity approaches zero andthe entire population will rapidly go extinct. We are thus notinterested in the behaviour of the system for α near unity, as the
population behaviour degenerates in this regime. For intermediate
values of α , however, our results are in excellent agreement withnumerically determined times (neutral case in Fig. 2, quasi-neutral
absorption times in Fig. 3, and quasi-neutral fixation times inFig. 4). Our results agree with simulations for carrying capacities
as small as fifty individuals.
In Fig. 1A, we compare Eq. (3) with numerical calculationsof the fixation time for a single mutant for large and smallvalues of N. In agreement with our analytic expressions, fixation
time increases linearly in N. The slopes of both predicted andnumerically determined lines are equal to the limit of machine
precision. The lines differ only in their intercept – the result of
lower order terms – which we have not determined for Eq. (3),but which are negative for the numerical results. Consequently weare slightly overestimating the fixation time by a fixed amount.
This remains true as the population size is decreased, until aroundtwenty individuals (Fig. 1B), where the numerical results deviate
positively from the linear relationship to intercept the origin.
Much of this effect can be attributed to events where one or both
species vanish before the population arrives at the centre manifold.
Numerical investigations reveal that using the harmonic mean
population size in place of carrying capacity in our expressions
above provides a slight improvement in the intercept, with thedifference diminishing as N increases. Subsequent determinationof the quasi-stationary distribution shows that the harmonic mean
population size converges to the carrying capacity with increasing
values of N (Parsons and Plotkin, in preparation).

Fig. 1. A comparison of the analytic approximation (dashed lines) for fixation timeof a single mutant, Eq. (3) with exact numerical results (solid lines), for varying
values of N, (A) α = 0.5, β2 = 1.0, and β1 = 1.0, 2.0, and 5.0, as indicated on
the graph. (B) detailed view of the same figure for small N.

3.2. Absorption and fixation times for K neutral types
Under the assumption of true neutrality, βi ≡ β and δi ≡ δ
for all i, the results of the previous section can be used to obtainfixation probabilities and times for an arbitrary number of types:
when the types are truly neutral, we can calculate the probabilityof fixation and expected time to fixation for type i by considering
it as one type, and collapsing all other types into a single, second
type. Using Eqs. (3) and (4), we have

πi (σi ) = σiTi (σi ) = −N ∗ (1 − σi ) ln(1 − σi )

β

σi

.

(6)

Moreover, the absorption time is simply
T (σ1 , . . . , σK ) =KXπi (σi )Ti (σi )

(7)

i =1

=

K
X
i =1

−

N∗

β

(1 − σi ) ln(1 − σi )(8)where σi is exactly equal to the initial density of type i:σi =XiN (0)KP.

XjN (0)

j =1

Dividing by the expected age to first reproduction, β1 , we get the
expected number of generations until absorption or fixation. Our
expressions agree with the classical expressions (Ewens, 2004) for
absorption and fixation times in a population of fixed size equal
to the carrying capacity, N ∗ , and are in excellent agreement withnumerically exact results (Fig. 2A). In other words, in the case ofstrict neutrality, classical results on fixation times in constantsize populations still hold even when the population size varies
stochastically.
As we remarked above, we expect our estimated fixation time
to be less accurate for values of α near unity — in which casethe population spends significant periods away from carrying
capacity. In Fig. 2B, we compare our expression with numericallyT.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

Fig. 2. A comparison of the analytic approximation (dashed lines) for fixation time
Eq. (6) in the neutral case (β1 = β2 = 1) with exact numerical results (solid lines).
(A) N varied at fixed α = 0.1, . . . , 0.9. (B) α varied at fixed N = 500, 1000, 1500.

exact absorption times over a range of values of α , for threedifferent values of N. We see that for the intermediate systemsize, N = 1000, our results remain accurate for α < 0.8,
corresponding to a carrying capacity of only 250 individuals.
At smaller carrying capacities, our analytical expression overestimates the fixation time.
3.3. Comparison with fixed population size models
For strictly neutral types, (i.e. βi ≡ β ), we have shown that
our expressions for mean absorption and fixation time coincidewith the classical expression for a population with fixed total size.
In this section, we compare absorption and fixation times in thequasi-neutral case against equivalent fixed-size models. Defining acorresponding fixed-size model for the quasi-neutral case is more
complicated than for the strictly neutral case, as the two types havedifferent birth and death rates. We address this by using two fixedsize models. Both models are formulations of the continuous time
Moran model (Ewens, 2004) with types that are neutral and with
fixed population size equal to the carrying capacity. In the firstmodel, we assume all types are identical, with common birth rateequal to the harmonic mean of the birth rates in our quasi-neutralmodel:β̄ =2
1/β1 + 1/β2.

(9)For this choice of birth rate, each individual has mean time
to first reproduction equal to the average reproduction time inan equilibrium quasi-neutral population. To keep the population
size fixed, we assume that each birth is immediately followedby the death of a randomly chosen individual. The diffusion
approximation to this process, obtained in the large population size
limit, is equivalent to a haploid Wright–Fisher model with effective
population size one half of census size. In this limit the fixation andabsorption times are given by Eq. (6) and Eq. (7) respectively, withβ as defined in Eq. (9).
The model above is the fixed population size equivalent ofour density-dependent model when the two types are identical.
However, this confounds two different effects. To investigate theeffect of fixed population size alone, we use a variant of the Moranmodel that separates the fitness contributions of birth and death
rates (Parsons and Quince, 2007a). This is done by allowing the305

Fig. 3. Time to absorption for two quasi neutral types as a function of initialfrequency, σ , of type 1. Total population size starts at carrying capacity, N ∗ . The
graph shows exact numerical results (solid line) and the asymptotic approximation
Eq. (2) (dashed line). Parameters are β1 = 5, β2 = 1, α = 0.5 and N = 1000.
Results from the equivalent neutral fixed size models, with identical types and birthrate equal to the harmonic mean (dotted line) and different birth rates (dot-dash)
are also shown.

types to differ in their birth rates, and also in their probabilityof being replaced subsequent to a birth event. To parameterise
this model in the quasi-neutral case, we simply take the birthrates to be equal to those in our quasi-neutral model, β1 and β2 ,and then weight the probability of replacement by the birth rates.
This ensures that the two types are neutral. It differs from thestandard Moran model in that the normalisation for the deaths
results in an extra frequency dependence in the rate at which
events occur. As a result, in the limit of large population size
there is no corresponding Wright–Fisher model. Since it is a
one dimensional stochastic process, the absorption and fixation
times can be written down explicitly and efficiently calculated
numerically (van Kampen, 1992).
In Fig. 3, we plot the time to absorption for two quasi-neutral
types as a function of initial frequency σ , starting with an initial
total population at carrying capacity (solid line). Comparing thesewith the identical-type fixed-size model (dotted line), we see thatthe quasi-neutral absorption times are skewed, with a longer mean
absorption time when the type with higher birth rate forms a larger
initial proportion of the population. Examining the equivalent
fixed size model with different types (dot-dashed line) we seethat this skew does not depend on a fluctuating population size,
but is rather a result of having types with different birth rates.
Fluctuations in the total population size reduce absorption times
when the type with higher birth rate is initially at low frequency.
This reflects the lower fixation probability for this type in afluctuating population, and the consequent increased likelihood of
its relatively rapid elimination from the population.
In Fig. 4, we show the mean time to fixation of type 1 andtype 2 for the density-dependent model and the two fixed-size
models. When the types are identical (dotted line), the fixation
time is independent of type and a decreasing function of initial
density. The effect of having types with different birth rates at fixed
population size (dot-dashed lines) is to increase the fixation times
for the higher birth rate type 1 and decrease them for the lower
birth rate type 2. This counter-intuitive effect is most pronounced
at intermediate initial densities. This suggests that the dynamics offixation are controlled not by the birth rate of the fixating type, butby the death rate of the displaced type. Comparing these results toour density dependent model (solid and dashed lines) it is apparent
that demographic stochasticity reduces fixation times for bothtypes, an effect that is more pronounced when the initial frequency
is smaller.

306

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

Fig. 4. Time to fixation for two quasi neutral types as a function of initial frequency
σ , starting from carrying capacity: exact numerical results (solid line) and theasymptotic approximation Eq. (2) (dashed line). Parameters are β1 = 5, β2 = 1,
α = 0.5 and N = 1000. The times are shown for both type 1 and type 2 (typeindicated on graph). Results from the equivalent neutral fixed size models, withidentical types and birth rate equal to the harmonic mean (dotted line) and different
birth rates (dot-dash) are also shown.

Fig. 5. Generations to fixation for each of two quasi-neutral types (type indicatedon graph) as a function of initial population size, starting from equal numbers ofeach type: exact numerical results (solid line) and the asymptotic approximation
Eq. (3) (dashed line). Parameters are β1 = 5, β2 = 1, α = 0.5, and N = 1000.

3.4. Dependence on initial population sizes
In Fig. 5 we compare the time to fixation for the two
types as obtained from expression Eq. (3) with exact numericalcomputations. The two agree very well for all but extremely
small population sizes, a consequence of very small populations
tending to go extinct before reaching carrying capacity. More
importantly, the figure shows the consequence of the trade-off
on fixation times: types that allocate toward higher birth-rates atthe cost of higher death rates will on average fix very rapidly inunsaturated environments, but will take substantially longer thana type with fewer births and deaths when the population is nearcarrying capacity. This result likely has important consequences
when the population is subject to exogenous fluctuations, such as
bottlenecks. In particular, repeated extreme bottlenecks are likely
to favour the fixation of types with higher vital rates.
4. Discussion
The expressions in Section 3.1 are, to the best of ourknowledge, the first analytic results for the time to absorption

and fixation in a two-allele model with stochastically varying
population sizes. These are asymptotic results, and they can be
expected to be precise only for populations with large carrying
capacity. Nonetheless, numerical experiments confirm that they
are quite accurate for populations of one hundred or fewer haploid
individuals.
Our study focused on the quasi-neutral case: both types havethe same carrying capacity and would be considered neutral
in a fixed population size model. We have shown that in thequasi-neutral case, the mean absorption and fixation times forintermediate sized populations are skewed relative to a modelwith identical types. This effect is independent of the presence of afluctuating population size. The effect of demographic stochasticity
is to reduce fixation times and absorption times. For quasi-neutral
mutants, fixation times are substantially smaller than for theequivalent fixed size model, something we had noted previously
for deleterious mutants (Parsons and Quince, 2007a).
Interestingly, we observed that when a single individual
invades a resident population at equilibrium, the expected time to
fixation does not depend upon whether the invader has a higher orlower birth rate than the resident. In other words, the advantage ofa high birth rate at low population density is exactly compensatedby the disadvantage of high death rates at large population density.
By contrast, when the number of invaders is significantly larger
than one, invaders with a higher birth rate take longer to fix than
invaders with a lower birth rate. In Maruyama and Kimura (1974)
and van Herwaarden and van der Wal (2002) a similar result was
found for the Wright–Fisher model with selection: a single mutant
that is negatively selected will on average take the same time to fix
as a positively selected mutant. More recently, Antal and Scheuring(2006) have shown that fixation times are equal in a populationof fixed size consisting of two types competing according to a
2 × 2 matrix game, while Taylor et al. (2006) show that for all onedimensional discrete time Markov chains on {0, . . . , N } satisfying
a detailed balance condition, if 0 and N are absorbing states, thenthe distribution of first exit times conditioned on first exit through0 or through N are equal. Inspired by their elegant argument,in Appendix B, we show that for diffusions on an interval I =(r0 , r1 ), subject to a similar condition, the conditional first exit time
distributions across r0 and r1 are equal. The result in Taylor et al.
(2006) and our result are consequences of time-reversibility (Kent,
1978; Ewens, 2004): after a relaxation period, for all pairs of points
x, y ∈ {1, . . . , N − 1} (x, y ∈ (r0 , r1 ) respectively), the probability
that the process following a path of length t from x to y is equal tothe probability that the process goes from y to x in equal time. This
explains the equality in expected fixation times observed above,
and also ensures that all higher moments of the fixation time are
equal as well.
Lastly, we found that in the limit of true neutrality, when
all types have the same birth and death rates, our expressions
for times to absorption and fixation simplify to those derived
in the neutral case by Ewens, Littler, and Kimura and Ohta
(Kimura and Ohta, 1969a,b; Littler, 1975; Ewens, 2004), with thecarrying capacity in place of population size. In other words,
our results are compatible with classical theory. In addition,
our results demonstrate that the theory of neutral types with afixed population size remains asymptotically correct in predicting
fixation times for stochastically varying populations. This is duein large part to the fact that the total population will rapidly
grow to carrying capacity and will thereafter spend a long period
fluctuating about carrying capacity. This is also true for the quasineutral case, which indicates that one must use caution in making
separation of time-scale arguments. The differences between theneutral and quasi-neutral cases arise from these fluctuations:in the quasi-neutral case, the size of the fluctuations aboutcarrying capacity are not uniform, but rather vary according to the
proportion of the two types, which leads to fixation probabilities
and times that differ qualitatively from models with a fixed
population size.

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

307

Acknowledgments

A.2. Exit probabilities and times

The authors thank Sergey Kryazhimskiy and three anonymous reviewers for comments that greatly improved the final
manuscript; we are particularly grateful to the latter for bringing
the studies (Antal and Scheuring, 2006; Taylor et al., 2006) to ourattention. Christopher Quince was supported by the University ofGlasgow. Todd L. Parsons and Joshua B. Plotkin were supported bya grant from the Burroughs Wellcome Fund. Joshua B. Plotkin also
acknowledges support from the James S. McDonnell Foundation.

motion of a particle in R≥0 . Extinction of a type i corresponds
to the particle exiting this region through one of the faces Γi ={x|xi = 0, xj 6= 0, j 6= i}.
The problem of exit from a domain has been long studied, and
it is a well known result (Gardiner, 2004) that the probability ofexiting through Γi starting from an initial point x, πi (x), satisfies
the Komolgorov equation

The diffusion process XN may be visualised as describing the

K

K
X

Appendix A. Deriving the asymptotic expressions

K1X

bi (x)∂xi πi (x) +2 i,j=1i =1aij (x)∂x2i πi (x) = 0(A.3)

In this appendix, we give a sketch of the mathematical details of
the derivation of equations Eqs. (2) and (3). Our argument proceedsby three steps. We

with boundary conditions(1) approximate the Markov chain XN by a diffusion process XN ,(2) express the expected absorption and fixation times as solutions to a Dirichlet problem,(3) and lastly, reduce the Dirichlet problem to an analyticallysolvable one-dimensional boundary value problem.

while the expected time to first exit, across any face, starting fromx, T (x) satisfies

πi |Γj = δij ,KXKbi (x)∂xi T (x) +i =1

T |Γj = 0

aij (x)∂x2i T (x) = −1

(A.4)

for all j.

Lastly, the expected time to exit, conditioned on exit through Γi , Ti ,satisfies

A.1. Diffusion approximation
For large N, we can approximate the process N −1 XN by adiffusion process XN (Ethier and Kurtz, 1986; van Kampen, 1992;
Gardiner, 2004). If XN has probability density p(y, t |x) – for all
A ⊆ RK ,Px {XN (t ) ∈ A} =2 i,j=1

with boundary condition

We elaborate on each of these steps below.

Z

1X

K
X

K

bi (x)∂xi (πi (x)Ti (x)) +

i =1

1X
2 i,j=1

aij (x)∂x2i (πi (x)Ti (x)) = −πi (x)(A.5)

with boundary condition

p(y, t |x) dy,

for all j.

(πi Ti ) |Γj = 0

Awhere Px denotes the probability conditioned on initial value x –then p must satisfy the Fokker–Planck equation
A.3. Projection on the centre manifold∂t p(y, t |x) = −

K
X

K
 1X


∂y2i aij (y)p(y, t |x) ,∂yi b (y)p(y, t |x) + i

2 i,j=1

i=1



1


−
1
σi e−βi τα!
xi = 
K −1
X
1


1−
σi e−βK τ α −1i=1

where
bi (y) = βi yi − δi yi

1+

N
X

!
yj

j =1


!
N
X


βi yi + δi yi 1 +yjNaij (y) =
j =1


0

(A.1)
if i = j,σ(x) =p(y, 0|x) = δ(x − y).
Under the assumption of quasi-neutrality, in the limit as N →∞, XN approaches a deterministic dynamical system
1−α

1+

N
X

!!
Xj

j =1

(A.2)



x ∈ RK |

K
X
j =1

xj =

1

α

)
−1 .

(A.6)
if i = K .

G(t , x) = −1
−1lim X(t , x),

1

α
Z ∞

t →∞

1−α

1+

t

K
X

!
Xj (s, x)

ds,

j =1

and

τ (x) = G(0, x).
Then if Ẏi = βi Yi has solution Y(t , x), a simple calculation showsthatwith an attracting centre manifold along the hyperplane(

if i < K

We obtain σ and τ via a change of time coordinates in the
dynamical system Eq. (A.2): let X(t , x) denote the solution startingfrom x at time t = 0, and define

otherwisewith initial conditionẊi = βi Xi

To proceed, we introduce co-ordinates σi and τ such thatX(t , x) = Y(G(t , x),



1

α



− 1 σ (x)).

(A.7)

Thus, the two systems follow the same solution curves, but atdifferent rates, with trajectories of Y beginning at the endpoints

308

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

of the trajectories of X at t = ∞. Now Yi (t , x) is simply an
exponential function; combining this with Eq. (A.7) and evaluating
at t = 0 yields Eq. (A.6). Intuitively, τ is a measure of the distancePKof the point x from the centre manifold i=1 xi = α1 − 1, while σ(x)maps each trajectory of Eq. (A.2) to a unique point in the standard
simplex

(
∆ =
K

K

x∈R |

K
X

Appendix B. Equality of fixation times in symmetric diffusions
In this appendix, we sketch a proof of the equality of the
distribution of fixation times for a class of diffusions X (t ) definedon an open interval I = (r0 , r1 ) that includes our process Σ (t ).
Let X have distribution p(t , y|x) satisfying the Fokker–Planck
equation)xj = 1 ,∂t p(t , y|x) =j =1

giving a convenient parameterisation of the flows. In the neutral
case, we have

σi =

xi
K
P

,

2

∂y2 [a(y)p(t , y|x)] − ∂y [b(y)p(t , y|x)]with a(y) > 0 on I, with absorbing boundary condition
lim W (y)a(y)p(t , y|x) = 0

xj

W (y) = ewhile in general, the dependence of σi on x is more complicated,
and the value of σ must be determined numerically.
In (Parsons and Quince, 2007b), we showed that by making achange of variable τ = √T and looking for an asymptotic seriesN1

expansion in powers of N − 2 , the solution, π1 , of Eq. (A.3) could be1

determined to order O(N − 2 ), from a lower dimensional problem
on the centre manifold. An identical argument applies here; werefer the reader to our previous paper for details.
We now focus on the case K = 2, for which we are able
to find analytical asymptotic expressions for the exit times andprobabilities. In this case, the fixation of type 1 corresponds to thefirst exit across Γ2 , while the centre manifold is simply the set ofpoints {(σ1 , τ )|0 ≤ σ1 ≤ 1, τ = 0}. After projection to the centre
manifold, our process reduces to a one dimensional diffusionΣ (t ) on (0, 1), while Eqs. (A.3)–(A.5) become ordinary differential
equations:

β1 β2 (1 − σ1 )σ1− 1 N (β1 σ1 + β2 (1 − σ1 ))2
α


dπ1
d2 π1
=0
× (β2 − β1 )+ (β1 σ1 + β2 (1 − σ1 ))dσ1dσ12
π1 (0) = 0
π1 (1) = 1
for the fixation probability (see (Parsons and Quince, 2007b),)−2

R b(y)
a(y)

dy

.

(B.2)

Thus, p(t , y|x) also satisfies Komologorov’s equation∂t p(t , y|x) =12

a(x)∂x2 p(t , y|x) − b(x)∂x p(t , y|x)with boundary condition
lim p(t , y|x) = 0

(B.3)

x→ri

(Feller, 1954a; McKean, 1956).
We assume additionally that a(y) and b(y) are in C 2 (I ), p(t , y|x)
is in C 2 ([0, ∞) × I × I ), and that for i = 0, 1b(y)0 < lim

a(y)

y→ri

< ∞.

(B.4)

Note that we do not require limy→ri a(y) or limy→ri b(y) to be nonzero or finite, so our assumptions include the process Σ (t ). We
refer the reader to Feller (1954a,b) and McKean (1956) for detailed
discussions of conditions under which these postulates apply.
For the exit time distributions to be sensible, we also require
both boundaries of I to be accessible according to Feller’sclassification scheme (Feller, 1954a,b; Ethier and Kurtz, 1986;
Ewens, 2004): letm(y) =

Z

1
a(y)W (y)dyandp(y) =

Z

W (y) dybe the speed measure and scale function, respectively. We thenrequire thatZβ1 β2 (1 − σ1 )σ1−
1
N (β1 σ1 + β2 (1 − σ1 ))2
α


dT
d2 T
× (β2 − β1 )+ (β1 σ1 + β2 (1 − σ1 )) 2 = −1
dσ1
dσ1
T (0) = T (1) = 0

u(y) =

for the expected time to absorption, andthen for all x ∈ I,

β1 β2 (1 − σ1 )σ11−1
N (β1 σ1 + β2 (1 − σ1 ))2
α


d
d2
× (β2 − β1 )(π1 T1 ) + (β1 σ1 + β2 (1 − σ1 )) 2 (π1 T1 )
dσ1
dσ1
= −π1
π1 (0)T (0) = π1 (1)T (1) = 0

Px τri < ∞ > 0.

for the expected time to fixation of type 1.
These linear equations may be readily solved to give theexpressions presented in 3.1.

where

1

(B.1)

y→ri

where

j =1

1

1

m(y) dp(y)

be finite in the limit as y → ri (so both boundaries are exit or
natural.) This implies that with non-zero probability, the process
will exit across one or the other boundary in finite time: if

τy = inf{t ≥ 0|X (t ) = y},It may be readily verified that both 0 and 1 are accessible
boundaries for Σ .
Under these assumptions, the diffusion is symmetric withrespect to w(y) (Kent, 1978):p(s, y|x)w(y)w(y) ==p(s, x|y)w(x)1a(y)W (y).

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

Moreover, the probability current,1

J (t , z |x) = lim

y→z 2

and, provided the limits in the numerator and denominator exist,∂y [a(y)p(t , y|x)] − b(y)p(t , y|x)lim ∂y a(y)w(y)is continuous and finite everywhere on Ī = [r0 , r1 ]. Note that we
require the use of a limit in our definition of the probability current
so that it can be applied at the endpoints r0 and r1 of I. We also notethat Eqs. (B.1) and (B.4) together imply that1

J (t , ri |x) = lim

y→ri 2∂y [a(y)p(t , y|x)] ,

=

y→r1



lim ∂y a(y)w(y)R∞lim ∂y a(y)w(y)

R∞

y→r1



y→r1

lim J (t , ri |x) = 0

for i 6= j.
Under our smoothness assumptions on p(t , y|x), we may
express the conditional exit-time distributions in terms of theprobability current J (Gardiner, 2004):0

J (s, ri |x) ds
J (s, ri |x) ds

lim Px {τ > t |X (τ ) = r1 } = lim Py {τ > t |X (τ ) = r0 } .
y→r1

(B.6)

The right hand side of this equality isR∞
lim Rt∞J (s, r1 |x) ds

J (s, r1 |x) ds
R∞a(x)w(x)J (s, r1 |x) ds= lim Rt∞.
x→r0a(x)w(x)J (s, r1 |x) ds
0
x→r0

0

Now, from Eq. (B.3), both numerator and denominator vanish as
x → r0 , so applying l’Hôpital’s rule,

= lim

x→r0

∂x

R ∞

∂x

R ∞

a(x)w(x)J (s, r1 |x) ds



a(x)w(x)J (s, r1 |x) ds

0Using Eq. (B.5),

∂x
= lim

x→r0

∂x
∂x

= lim

x→r0

∂x


R∞
t



lim 21 ∂y [a(y)a(x)w(x)p(s, y|x)] ds

y→r1


R∞
0

t

0

lim

∂



( ) ( )w(x)p(s, y|x)] ds

1[a y a x
2 y

y→r1


R∞

∂

lim

1
[a y a x
2 y

y→r1


R∞

( ) ( )w(y)p(s, x|y)] ds

lim 21 ∂y [a(y)a(x)w(y)p(s, x|y)] ds




y→r1

p(s,y|x)p(s,x|y)where we have used the symmetry of X (t ): w(y) = w(x) . We
have assumed p(t , x|y), a(x), and a(y)w(y) are twice continuously
differentiable, and their product has a continuous and bounded(and thus uniformly continuous) extension to Ī, so we caninterchange limits, derivatives and integrals with impunity toobtainlim ∂y a(y)w(y)R∞ 1lim ∂y a(y)w(y)

R∞ 1



= lim

x→r0

y→r1



y→r1

t
0

0

,

lim 12 ∂x [a(x)p(s, x|y)] dsx→r0J (s, r0 |y) ds

J (s, r0 |y), ds

.

∂( ) ( , x|y)] ds∂( ) ( , x|y)] ds[a x p s
2 x[a x p s
2 x

= lim

a(y)w(y)R∞
Rt∞

J (s, r0 |y) ds( )w(y) 0 J (s, r0 |y) ds= lim Py {τ > t |X (τ ) = r0 } ,y→r1 a yy→r1

12∂y2 [a(y)w(y)] − ∂y [b(y)w(y)] = 0,

so w(y), when integrable, is a stationary distribution for X (t ),which is the case provided limy→ri a(y) > 0. When a(y) degenerates at the boundaries, w(y) is often a good approximation to thequasi-stationary distribution for X (t ), except near to the endpointsof I (Gardiner, 2004) — in this case, the only true stationary
distribution is 0, as X will exit through one of the endpoints withprobability 1.
We may use this fact to give an intuitive interpretation of Eq.
(B.6): after a transient period, the distribution of X converges tow(x), except perhaps very near the boundaries. Thus, X is timereversible with respect to w : the probability of moving from x
to y in time t is w(x)p(t , y|x), which, by symmetry, is equal tothe probability w(y)p(t , x|y) of X moving from y to x in the sametime (Kent, 1978). In particular, the probabilities of the process
traversing a path from r0 to r1 in time t or from r1 to r0 in the sametime are equal, yielding Eq. (B.6).
References

.

t

t

x→r0

yielding our desired result.
Note that.

J (s, ri |x) is the flow of probability across ri at time s, i.e. theprobability that the diffusion process will exit in [s, s + ∆s). Thus,
the numerator is the cumulative probability of exit at ri for all timess > t, while the denominator is the flow across ri for all time,
and thus the probability of exit across ri ; the fraction gives theconditional probability of exit across ri at some time after t. We
will use this expression to show that for a symmetric diffusion,
x→r0

0



lim 12 ∂x [a(x)p(s, x|y)] ds

Since the flux is well defined at r0 , this last line allows us toconclude that the limits do indeed exist. We conclude by applying
l’Hôpital’s rule again:x→rj

Px {τ > t |X (τ ) = ri } = Rt∞tlim ∂y a(y)w(y)



=

R∞
R∞

y→r1

(B.5)

while from Eq. (B.3) we have

R∞

309


,

Antal, T., Scheuring, I., 2006. Fixation of strategies for an evolutionary game in finite
populations. Bull. Math. Biol. 68, 1923–1944.
Ethier, S.N., Kurtz, T.G., 1986. Markov Processes: Characterization and Convergence.
Wiley, New York.
Ewens, W.J., 1967. The probability of survival of a new mutant in a fluctuating
environment. Heredity 22, 438–443.
Ewens, W.J., 2004. Mathematical Population Genetics. Springer-Verlag, New York.
Feller, W., 1954a. Diffusion processes in one dimension. Trans. Amer. Math. Soc. 77(1), 1–31.
Feller, W., 1954b. The general diffusion operator and positivity preserving semigroups in one dimension. Ann. of Math. 60 (3), 417–436.
Fisher, R.A., 1958. The Genetical Theory of Natural Selection, Second ed.. Dover, NewYork.
Gardiner, C.W., 2004. Handbook of Stochastic Methods for Physics, Chemistry and
the Natural Sciences. Springer, Berlin Heidelberg.
Grasman, J., van Herwaarden, O.A., 1999. Asymptotic Methods for the
Fokker–Planck Equation and the Exit Problem in Applications. Springer,
Berlin Heidelberg.
Kent, J., 1978. Time-reversible diffusions. Adv. Appl. Prob. 10, 819–835.
Kimura, M., 1957. Some problems of stochastic processes in genetics. Ann. Math.
Statist. 28, 882–901.
Kimura, M., 1962. On the probability of fixation of mutant genes in a population.
Genetics 47, 713–719.
Kimura, M., Ohta, T., 1969a. The average number of generations until the extinctionof an individual mutant gene in a finite population. Genetics 63 (3), 701–709.
Kimura, M., Ohta, T., 1969b. The average number of generations until the fixation ofa mutant gene in a finite population. Genetics 61 (3), 763–771.
Kimura, M., Ohta, T., 1974. Probability of gene fixation in an expanding finite
population. Proc. Natl. Acad. Sci. 71, 3377–3379.

310

T.L. Parsons et al. / Theoretical Population Biology 74 (2008) 302–310

Lambert, A., 2006. Probability of fixation under weak selection: A branching process
unifying approach. Theor. Pop. Biol. 69, 419–441.
Littler, R.A., 1975. Loss of variability at one locus in a finite population. Math. Biosci.
25 (1-2), 151–163.
MacArthur, R.H., Wilson, E.O., 1967. The Theory of Island Biogeography. Princeton
University Press, Princeton, NJ.
Maruyama, T., Kimura, M., 1974. A note on the speed of gene freqency changes inreverse directions in a finite population. Evolution 28, 161–163.
McKean, H.P., 1956. Elementary solutions for certain parabolic partial differential
equations. Trans. Amer. Math. Soc. 82 (2), 519–548.
Moran, P.A.P., 1958. Random processes in genetics. Proc. Camb. Phil. Soc. 54, 60–71.
Newman, T.J., Ferdy, J.-B., Quince, C., 2004. Extinction times and moment closure in
the stochastic logistic process. Theor. Pop. Biol. 65 (2), 115–126.
Otto, S.P., Whitlock, M.C., 1997. The probability of fixation in populations ofchanging size. Genetics 146, 723–733.

Parsons, T.L., Plotkin, J.B., Quasi-stationary distributions for quasi-neutral andneutral populations exhibiting density dependence and symmetric mutation.
(in preparation).
Parsons, T.L., Quince, C., 2007a. Fixation in haploid populations exhibiting density
dependence I: The non-neutral case. Theor. Pop. Biol. 72 (1), 121–135.
Parsons, T.L., Quince, C., 2007b. Fixation in haploid populations exhibiting density
dependence II: The quasi-neutral case. Theor. Pop. Biol. 72 (4), 468–479.
Pianka, E.R., 1972. r and K selection or b and d selection?. Am. Nat. 106, 581–588.
Taylor, C., Iwasa, Y., Nowak, M.A., 2006. A symmetry of fixation times in evolutionary
dynamics. J. Theor. Biol. 243, 245–251.
van Herwaarden, O.A., van der Wal, N.J., 2002. Extinction time and age of an allele
in a large finite population. Theor. Pop. Biol. 61, 311–318.
van Kampen, N.G., 1992. Stochastic Processes in Physics and Chemistry. Elsevier,
Amsterdam, San Diego, Oxford, London.
Wright, S., 1931. Evolution in Mendelian populations. Genetics 16, 97–159.