INVESTIGATION

Evolutionary Branching in a Finite Population:
Deterministic Branching vs. Stochastic Branching
Joe Yuichiro Wakano*,†,1 and Yoh Iwasa‡*Meiji Institute for Advanced Study of Mathematical Sciences, Meiji University, Kawasaki 214-8571, Japan, †Japan Science andTechnology, Precursory Research for Embryonic Science and Technology, Tokyo 102-8666, Japan, and ‡Department of Biology,
Kyushu University, Fukuoka 812-8581, Japan

ABSTRACT Adaptive dynamics formalism demonstrates that, in a constant environment, a continuous trait may ﬁrst converge toa singular point followed by spontaneous transition from a unimodal trait distribution into a bimodal one, which is called “evolutionary
branching.” Most previous analyses of evolutionary branching have been conducted in an inﬁnitely large population. Here, we study
the effect of stochasticity caused by the ﬁniteness of the population size on evolutionary branching. By analyzing the dynamics of trait
variance, we obtain the condition for evolutionary branching as the one under which trait variance explodes. Genetic drift reduces thetrait variance and causes stochastic ﬂuctuation. In a very small population, evolutionary branching does not occur. In larger populations, evolutionary branching may occur, but it occurs in two different manners: in deterministic branching, branching occurs quickly
when the population reaches the singular point, while in stochastic branching, the population stays at singularity for a period before
branching out. The conditions for these cases and the mean branching-out times are calculated in terms of population size, mutational
effects, and selection intensity and are conﬁrmed by direct computer simulations of the individual-based model.

I

N addition to models using quantitative genetics (Lande 1981;
Barton and Turelli 1991; Iwasa et al. 1991; Pomiankowski
et al. 1991), adaptive dynamics formalism has revealed diverse behaviors of evolutionary dynamics associated with
continuous traits in various organisms (Metz et al. 1992;
Dieckmann and Law 1996; Geritz et al. 1997, 1998). Thisapproach allows us to analyze the evolution of many ecological traits for which frequency-dependent selection isprevalent without specifying the derlying genetic systems
controlling the trait. Even in a perfectly constant environment, a population may show intriguing temporal behavior.
For example, if the trait evolves by accumulation of mutations of small magnitude and if the ﬁtness is frequency
dependent, a continuous trait may show convergence to a
singular point followed by spontaneous splitting of a unimodal
trait distribution into a bimodal (or multimodal) one, referred
to as “evolutionary branching” (Metz et al. 1992, 1996; Geritz
et al. 1997). Evolutionary branching is predicted to occur

Copyright © 2013 by the Genetics Society of America
doi: 10.1534/genetics.112.144980
Manuscript received August 17, 2012; accepted for publication October 22, 2012
Available freely online through the author-supported open access option.
1
Corresponding author: Meiji Institute for Advanced Study of Mathematical Sciences,
Meiji University, 1-1-1 Higashi Mita, Tama-ku, Kawasaki 214-8571, Kanagawa, Japan.
E-mail: joe@meiji.ac.jp

at an evolutionarily singular point that is approachingstable (or convergence stable, CS) but not evolutionarily
stable (ES), and it is actually observed in individual-based
simulations in many models for the evolution of ecological traits.
In evolutionary branching, an evolving population spontaneously changes its trait distribution from unimodal tobimodal (or multimodal). Although evolutionary branchinghas been regarded as a model of sympatric speciation bysome authors (Dieckmann and Doebeli 1999; Doebeli andDieckmann 2000), the relationship between evolutionary
branching in adaptive dynamics formalism and speciation
events has been the focus of much debate (for review, see
Waxman and Gavrilets 2005).
Most previous analyses of adaptive dynamics haveassumed an inﬁnite population size and that the ﬁtness
landscape acts as the deterministic driving force of evolution. However, in a ﬁnite population, stochastic ﬂuctuations
are unavoidable. For example, an individual with a lower
ﬁtness may be able to produce more offspring just by
chance. Thus, a ﬁnite population can evolve against the
selection gradient. As a result, stochasticity allows a population to jump from a low local peak to a higher peak. Hence,stochastic ﬂuctuation can be an important factor in theoutcome of evolution. The importance of stochasticity, orGenetics, Vol. 193, 229–241

January 2013229

genetic drift, has been studied in detail in the context ofquantitative population genetics (e.g., Lande 1976; Tazzyman
and Iwasa 2000). On the other hand, in population genetics, random drift is known to decrease the genetic diversityof a population (e.g., Ewens 2004), which tends to slow
down the action of natural selection. Hence, the effect ofgenetic drift is important to understand evolutionary dynamics, including evolutionary branching, in ﬁnite populations. In fact, Claessen et al. (2007) analyzed evolutionary
branching in a model including three species and reported
that genetic drift tends to cause a delay in evolutionary
branching. However, a more detailed study of the relationship between evolutionary branching and genetic drift isneeded.
In this article, we study the effect of stochastic ﬂuctuation
caused by the ﬁniteness of population size (random genetic
drift) on evolutionary branching. As an illustrative example,
we adopt a model of a public goods game studied by Doebeli
et al. (2004), who investigated evolutionary branching byusing a corresponding individual-based simulation. The gamemodel deals with a situation that is both biologically andsocially interesting: a self-organized differentiation of players into a highly cooperative group and a noncooperative
group. Concerning the genetics, we in effect assume thatthe population is haploid and that the focal trait is controlled by a locus with many alleles. By focusing on thissimplest case, we can analyze the effect of random genetic
drift on the evolutionary outcome most clearly. We ﬁrst
show the results of individual-based simulations, using different population sizes. We show that genetic drift is especially important when the population size is small and thatevolutionary branching does not necessarily occur evenwhen the standard theory of adaptive dynamics predicts thatit should. We also derive deterministic and stochastic versions of the dynamics of the trait variance. We reveal that
there are two different manners in which evolutionary branching occurs: in deterministic branching, branching may occur
quickly after the population reaches the singular point, whilein stochastic branching, the population stays at a singularity
for a period before branching out suddenly. We derive theconditions and mean branching-out times for these cases andconﬁrm our predictions by direct computer simulations of theindividual-based model.

Evolutionary Branching in a Game in a Finite Population
Individual-based model

As a full stochastic model of a population under natural
selection, mutation, and random genetic drift, we consider
the following stochastic process. The population consists of
N individuals reproducing asexually. Each individual has
a genetically determined trait value z 2 ½0; 1. At each discrete time step t, an individual i obtains a fecundity Fi frompairwise interactions with the other N 2 1 individuals. Eachinteraction is a two-player game where a payoff given to an

230

J. Y. Wakano and Y. Iwasa

individual i matched with an individual j is given by f ðzi ; zj Þ.
Its fecundity Fi is denoted byX
Fi ¼
f ðzi ; zj Þ:
(1)
j6¼i

Our framework is valid for any payoff function f ðzi ; zj Þ. As an
example, we adopt a model by Doebeli et al. (2004) inwhich trait z represents the amount of investment into public goods shared by both players. In this game, the payoff
function is represented byf ðzi ; zj Þ ¼ 1 þ Bðzi þ zj Þ 2 Cðzi Þ;(2a)where the beneﬁt function isBðz1 þ z2 Þ ¼ b1 ðz1 þ z2 Þ þ b2 ðz1 þ z2 Þ2 ;(2b)and the cost function isCðzÞ ¼ c1 z þ c2 z2 :(2c)

To keep the population size constant, the ﬁtness of anindividual with trait zi is denoted aswðzi Þ ¼

FiP :ð1=NÞ j Fj

(3)When the population size is inﬁnitely large, Doebeli
et al. (2004) have shown that the singular strategy
z* ¼ ðc1 2 b1 Þ=ð4b2 2 2c2 Þ is convergence stable when
QCS ¼ 4b2 2 2c2 , 0 and evolutionarily stable when
QES ¼ 2b2 2 2c2 , 0. The singular strategy is referred to as
an evolutionary branching point when QCS , 0 and QES . 0.
They have observed evolutionary branching in an individualbased simulation of a large population (N = 10,000).
In a ﬁnite population, the model includes random genetic
drift, and the details of the updating rule and the order ofdifferent stages in population dynamics become important. Here, we assume a Wright–Fisher process. At eachtime step, fecundity is calculated for all N individuals.
Then all individuals die and are replaced by their offspring, which are sampled multinomially, with the probability of an individual being the parent of each offspring
individual proportional to the fecundity of the focal
individual.
After reproduction, there is a small chance of mutation m,
and the mutant has a trait value slightly different from thatof the parent. To be speciﬁc, it follows a normal distribution
with a mean equal to the parent’s and a variance s2 . Mutants
are constrained to be between 0 and 1, which conﬁnes trait z
within an interval [0, 1].
Evolutionary branching depends on the population size

All simulations were run with parameter values for whicha unique singular point exists that is convergence stable andevolutionarily nonstable (i.e., the evolutionary branching

Figure 1 Evolutionary dynamics. (A) A simulation run with N = 8000. (B and C) Two different simulation runs with N = 1000. (D) A simulation run with
N ¼ 200. Parameters: m 5 0:01; s 5 0:02; ðb1 ; b2 ; c1 ; c2 Þ 5 ð6:0; 21:4; 4:56; 21:6Þ.

point). Thus, evolutionary branching is expected when the
population size is inﬁnitely large.
When the population size is large (N = 8000), evolutionary branching occurs as soon as the trait evolves to reach the
CS value z* (Figure 1A). For an intermediate population size
(N = 1000), the population lingers around the CS value, z*,before it ﬁnally shows evolutionary branching. The waitingtime until the evolutionary branching varies considerably
among different simulation runs with the same parametervalues (Figure 1, B and C). For a small population size
(N = 200), we do not observe evolutionary branching
and the population seems to stay around z*, even thoughz* is not evolutionarily stable (Figure 1D) (see also Wakano
and Lehmann 2012).
The evolutionary branching dynamics are characterized
by an abrupt increase or explosion of the variance of the
focal phenotype (the trait variance). If all the parameters forthe interactions are the same, a small population does notbranch, an intermediately large population shows delayedbranching, and a large population branches immediately.
The intuitive reason for the smallness of the population size
to suppress evolutionary branching is that genetic diversity
is reduced by genetic drift every generation. In the following
sections, we provide approximate dynamics for the trait
variance to understand these simulation results.

Dynamics of Trait Variance
The full (individual-based) model is analytically intractable.
Therefore, we derive a reduced simpliﬁed dynamical systemhere, which is tractable and gives predictions on the behaviorof the full model. One generation in the reduced model
consists of two substeps.
Substep 1: Natural selection and mutation

Let uðzÞ denotea distribution of trait z, which is normalizedR
to satisfy uðzÞdz ¼ 1. When a ﬁtness function wðzÞ is given
and selection and mutation are sufﬁciently weak, the trait
distribution in the offspring pool after mutation, denoted byxðzÞ, can be approximated by the following replicator equation plus a mutation term,
ms2 @ 2 uðzÞ þuðzÞ;xðzÞ ¼ uðzÞ þ wðzÞ 2 w
2 @z2

(4)

 is the population average of ﬁtness values. The
where wsecond term on the right-hand side represents the effect of
selection and the third term represents the effect of mutation. The Taylor expansion of the ﬁtness function around theaverage trait value z is denoted aswðzÞ ¼ wðzÞ þ w1 ðz 2 zÞ þw2
ðz2zÞ2 ;2(5)where the coefﬁcients w1 and w2 are determined by thecurrent state z. The coefﬁcient w1 represents the selection
gradient. The coefﬁcient w2 indicates stabilizing selection if
it is negative and indicates disruptive selection if it is positive. The derivation of w2 from the payoff function f ðzi ; zj Þ is
explained in Appendix A. For parameter values used in our
simulations, we have w2  0:1, implying that disruptive selection is operating. The population average and the variance of trait values in the offspring trait distribution xðzÞ arerespectively denoted as

Evolutionary Branching in a Finite Population

231

Avx ½z ¼ z þ w1 m2 þ
Varx ½z ¼ m2 2 ðw1

2

w2
m3 ;
2

(6a)

w2
w2
m4 þ ms2 ;þ Þm2 2 þ w1 m3 þ
2
2
(6b)

where mk denotes the kth centralized moment of uðzÞ. The
above approximation is derived by neglecting m5 and higherorder moments (see Appendix B).
Substep 2: Random sampling (or random genetic drift)

From the offspring pool xðzÞ, N individuals are independently sampled to become adults. Hence, even for a ﬁxed
xðzÞ, different distributions of adult phenotypes are realizedaccording to a probability distribution (probability measure)over all possible trait distributions. After sampling N individuals independently from xðzÞ, the expectation of the average
trait value is identical to the average trait value in xðzÞ. In
contrast, the expectation of the trait variance is decreased
by a factor (N 2 1)/N (Ewens 2004). Hence, we have thefollowing:
E½Av½Z ¼ Avx ½z
E½Var½Z ¼

N 21
Varx ½z:N(7)

Combining the two substeps, the expected changes in the
population average and the variance in one generation are
given by
E½DAv½Z ¼ w1 m2 þ

E½DVar½Z ¼

w2
m3 ;
2

(8a)

i
N 21h
w2
m2 þ w1 m3 2 w1 2 m2 2 þ ðm4 2 m2 2 Þ þ ms2 2 m2 :
2
N
(8b)

We are particularly interested in the dynamics of the varianceof the trait distribution. To analyze Equation 8 further, wehere assume that the trait distribution is close to a normal
distribution whose mean and variance obey the dynamics
given by Equation 8. For normal distributions, m3 ¼ 0 andm4 ¼ 3m2 2 . Using these relationships, we have the following recursive equations for the mean and variance of thetrait.:

Dm2 ¼

Dz ¼ w1 m2

(9a)

N 21
1
N 21 2
ðw2 2 w1 2 Þm2 2 2 m2 þ
ms :
N
N
N

(9b)

The ﬁrst, second, and third terms on the right-hand side ofEquation 9b represent the effect of natural selection, the lossof diversity due to genetic drift, and the increase of diversity
due to mutation, respectively. Without selection, the variance m2 converges to ðN 2 1Þms2 .

232

J. Y. Wakano and Y. Iwasa

Condition for the Explosion of Trait Variance
As long as the selection gradient does not vanish (w1 6¼ 0),
the mean trait value changes according to Equation 9a, which
is a standard result of population genetics (Ewens 2004) andcorresponds to the canonical equation in adaptive dynamics
theory (Dieckmann and Law 1996; Champagnat et al. 2006).
If a singular point z*, where the selection gradient vanishes(w1 ¼ 0), is convergence stable, then the system approaches
the singular point with time (Geritz et al. 1997). Once themean trait (i.e., the ﬁrst moment of trait distribution) reaches the singularity, it will not change any longer. However,
the trait variance (i.e., the second centralized moment of the
distribution) may change with time. See Figure 2 for a schematic illustration. Here we focus on the recursion for thetrait variance,
Dm2 ¼

N21
1
N21 2
w2 m2 2 2 m2 þ
ms ;
N
N
N

(10)

which indicates that the one-generational change in the trait
variance is equal to the sum of three terms: the natural selection term, the genetic drift term, and the mutation term.
The right-hand side of Equation 10 represents the “force”
acting on the trait variance (Figure 3). Note that the ordersof the three terms in Equation 10 are m2 2 , m2 , and 1, respectively. When the variance is small, the mutation (the thirdterm) is dominant, and it increases the trait variance. As thevariance increases, the genetic drift (the second term) increases its magnitude and it reduces the trait variance. Thebalance between these two processes would result in theequilibrium m2 *  Nms2 , which represents the variance determined by mutation–drift balance. When the trait varianceis large, selection (the ﬁrst term) becomes dominant. Selection decreases the trait variance when w2 , 0 (stabilizing
selection), but it increases the trait variance when w2 . 0(disruptive selection).
Combining these, we have three different cases, depending on the parameters.
Case i: When w2 , 0, stabilizing selection is at work. Then the
dynamics of the trait variance (Equation10) have a single
positive equilibrium
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1 þ 1 2 4ðN21Þ2 w2 ms2;
m2 *52ðN 2 1Þw2which is globally stable, and the trait variance converges tothis value for any initial value. At this equilibrium, the trait
variance is maintained by mutation, which increases variance,and by genetic drift and stabilizing selection, which decrease
variance. Especially when the stabilizing selection is weak,
trait variance converges to the level maintained by mutation–
drift balance (m2 *  ðN 2 1Þms2 ) when jw2 j is small.
Case ii: When 0 , w2 , 1=4ms2 ðN21Þ2 , disruptive selection
is at work, but it is not very strong. Then the dynamics ofthe trait variance have two positive equilibria: the smaller

Figure 2 A schematic illustration of phenotype distribution uðzÞ andﬁtness function wðzÞ when the mean trait value z has reached the singular strategy z*. The trait variance is Var½uðzÞ 5 m2 so the standard
pﬃﬃﬃﬃﬃﬃﬃ
deviation is m2 .

one is stable, but the larger one is unstable. If the trait
variance starts from a small value, it converges to thesmaller equilibrium, where the trait variance is maintained
mainly by the mutation–drift balance. If the trait variance starts from a larger value, then disruptive selection
dominates, leading to the unlimited growth of the trait
variance.
Case iii: When w2 . 1=4ms2 ðN21Þ2 , disruptive selection isstrong. The dynamics of the trait variance no longer havea stable equilibrium. Starting from any initial value, thetrait variance keeps increasing without limit.
Figure 4 shows a bifurcation diagram of the variance dynamics given by Equation 10, with the stabilizing/disruptive
selection coefﬁcient (w2 ) as a bifurcation parameter. Figure
5 shows a bifurcation diagram in which the population size(N) is a bifurcation parameter. We have a similar bifurcation
diagram where mutational effects (ms2 ) serve as a bifurcation
parameter (not shown). We regard the explosion of the trait
variance in the dynamics of Equation 10 as evolutionary
branching. Then, these bifurcation diagrams show that evolutionary branching occurs as a result of a saddle-node bifurcation. Evolutionary branching is predicted to occur at
a non-ES singular point when disruptive selection intensity,
population size, mutation rate, and mutation step size are
sufﬁciently large. Speciﬁcally, the deterministic dynamics
of Equation 10 predict that evolutionary branching occurs
when 4ms2 ðN21Þ2 w2 . 1 holds.

Stochastic Calculation of the Explosion of Variance
The deterministic dynamics of the trait variance represented by Equation 10 indicate only the expectation of onegenerational change. There must be stochasticity caused by
the genetic drift. Results of an individual-based model, suchas those illustrated in Figure 1, clearly show stochasticity.
To account for stochasticity, the deterministic differenceof Equation 10 needs to be modiﬁed by an additional term.
In this section, we derive the stochastic differential equa-

Figure 3 “Force” plotted against the trait variance based on Equation
10. Dotted, solid, and thick solid curves represent stabilizing selection(w2 ¼ 20.2), weak disruptive selection (w2 = 0.2), and strong disruptive
selection (w2 ¼ 0.6), respectively. They also correspond to cases i, ii, andiii in the main text. Parameters: N = 400, m 5 0:01; s 5 0:02.

tion model considering the stochasticity of trait variancedynamics.
The deterministic dynamics of Equation 10 can be approximated by the ordinary differential equation (ODE)dm2 ¼ f ðm2 Þ;
dt

(11a)

ðN 2 1Þw2 x2 2 x þ ðN 2 1Þms2
:
N

(11b)

where
f ðxÞ ¼

If stochasticity has no correlation over generations, and if
the paths can be approximated as continuous, then we canincorporate the Brownian motion term, and the differential
Equation 11 becomes the stochastic differential equation
(SDE)
rﬃﬃﬃﬃ
2
M2 dB;
dM2 ¼ f ðM2 Þdt þN(12)

where dB denotes Brownian motion and M2 denotes a random variable representing the trait variance. f ðM2 Þ is thesystematic change per generation in variable M2 and is given
by Equation 11b.
The second term on the right-hand side of Equation 12
represents stochasticity in trait variance caused by random
genetic drift (see Appendix C for derivation). The stochastic
ﬂuctuation term depends on the trait variance itself (geometric Brownian motion), because the amplitude of ﬂuctuation in the trait variance in each generation of offspring islarge when the trait variance in each parental generation islarge. The stochastic term also depends on population size,
because ﬂuctuation is averaged out in large populations.
Hence, random genetic drift becomes less important in larger
populations. The system is similar to but different from the

Evolutionary Branching in a Finite Population

233

Figure 4 Stable (solid) and unstable (dashed) branches of the trait variance m2 with disruptive selection intensity w2 as a bifurcation parameter.
Arrows indicate the direction of systematic change in m2 . Parameters: N =1000, m 5 0:01; s 5 0:02.

motion of a particle moved by force shown in Figure 3 underheat noise.
Stationary distribution

The corresponding Fokker–Planck equation (or forward Kolmogorov equation) to Equation 12 is given by@@
1 @2 2
pðm2 ; tÞ ¼ 2
m2 2 pðm2 ; tÞ
½ f ðm2 Þpðm2 ; tÞþ2@t@m2
2 @m2 N

(13)

(Feller 1971), where pðm2 ; tÞdm2 ¼ Pr½m2 # M2 ðtÞ , m2 þ
dm2  is the probability density for the trait variance m2 .
The stationary distribution of the trait variance is~ðxÞ ¼ C x23 eðN21ÞaðxÞ ;
p

(14a)

where C is a normalizing constant andaðxÞ ¼

w2 x2 2 ms2
:
x

(14b)

Figure 6 illustrates an example of stationary distribution of
the trait variance, showing that the system tends to show
slightly smaller trait variance than the locally stable value inEquation 11, because the ﬂuctuation is weaker when thevariance is smaller (geometric Brownian motion).
Comparison with the results of simulations: We compared
the stationary distribution of the trait variance between
the theoretical prediction given by Equation 14 and theobserved distribution obtained from the individual-based
simulation (Evolutionary Branching in a Game in a Finite
Population section). Figure 7 shows good agreement between the diffusion approximation (i.e., the SDE model)
and the corresponding simulation result with a small population size. In particular, the approximation captures theshift of the mode—the peak of the distribution lies at thevalue smaller than Nms2 ¼ 0:0008. The average trait variance in an interval shown in Figure 7 is close to Nms2 forboth analytic prediction (Equation 14) and simulation
results.

234

J. Y. Wakano and Y. Iwasa

Figure 5 Stable (solid) and unstable (dashed) branches of the trait variance
m2 with population size N as a bifurcation parameter. A Log plot is shown.
Arrows indicate the direction of systematic change in m2 . Parameters:
m 5 0:01; s 5 0:02; ðb1 ; b2 ; c1 ; c2 Þ 5 ð6:0; 2 1:4; 4:56; 2 1:6Þ.

Waiting time for evolutionary branching

In case ii, there are stable and unstable equilibria for thetrait variance, the latter being larger than the former. If theinitial value of trait variance is smaller than the unstableequilibrium, the deterministic model Equation 11 predicts
that variance would be maintained indeﬁnitely. In contrast,
if the initial value is larger than the unstable equilibrium, it
quickly increases and grows without limit. Hence, whether
trait variance eventually explodes depends on the initialcondition. In the presence of stochasticity, however, the trait
variance would not be maintained at smaller values than theunstable equilibrium; instead, it would ﬂuctuate around thestable equilibrium for a period, but then eventually increasebeyond the unstable equilibrium and diverge to inﬁnity.
Hence, under the stochastic differential equation, the trait
variance is predicted to explode eventually even if the initialvalue is small. The timing of this explosion is controlled by
the stochasticity.
The mean waiting time prior to explosion is a typical“ﬁrst passage time” problem in diffusion theory. This has
also been used in calculating the mean time of ﬁxation in
population genetics (Kimura and Ohta 1968) and the meantime to extinction in conservation biology (e.g., Lande 1993;
Hakoyama and Iwasa 2000) and physiology (e.g., Ricciardi
and Smith 1977).
To study the expectation of waiting time, we assume thata system has a reﬂecting boundary at m2 ¼ l  0 and anabsorbing boundary at m2 ¼ r (note that the left boundary
can never be achieved because the deterministic force is increasing the variance while stochastic ﬂuctuation vanishes
at m2 ¼ 0; i.e., Equation 12 is singular at m2 ¼ 0). Let TðhÞbe the expected waiting time until the system reaches theright absorbing boundary m2 ¼ r when the initial value ofthe trait variance is m2 ð0Þ ¼ h. To obtain TðhÞ, we use thebackward Kolmogorov equation
d
h2 d2
TðhÞ;
21 ¼ f ðhÞ TðhÞ þ
dh
N dh2

(15)

~ðm2 Þ of the trait variance m2 (EquaFigure 6 The stationary distribution ption 14). A Log-Log plot is shown, and thus the probability density ismuch higher at the peak than it appears. The locally stable value of thevariance according to Equation 11 is m2 *  0.00081, which is close toNms2 = 0.0008 (dashed line). Note that a local peak lies at a value,0.0008 because of geometric Brownian motion. Parameters: N = 200,m 5 0:01; s 5 0:02; ðb1 ; b2 ; c1 ; c2 Þ 5 ð6:0; 2 1:4; 4:56; 2 1:6Þ.

with the boundary conditions TðrÞ ¼ 0 and T9ðlÞ ¼ 0 (seeAppendix E). The solution is
Z rZ yx23 yeðN21ÞðaðxÞ2aðyÞÞ dxdy;(16)TðhÞ ¼ Nhlwhere aðxÞ is given by Equation 14b. The expected waitingtime is smaller when the population size is larger, disruptive
selection is stronger, or mutational effects are larger (seeAppendix D).
Comparison with the results of simulations: We compared
the mean waiting times observed in the computer simulationsof the individual-based model and the prediction by the approximate stochastic differential equation model given by
Equation 16. Once the system reaches a convergence stable
point, the trait variance m2 is expected to be close to the
locally stable value Nms2 (see also our results on stationary
distribution). Thus, we numerically calculate TðNms2 Þ andcompare the corresponding result in individual-based simulation (Figure 8). We also derive the probability distributionof the ﬁrst passage time by using Kolmogorov’s backward
equation in Appendix E, but the ﬁt of the shape of the probability distribution of waiting times to the distribution observed in individual-based simulations is less than perfect(Figure 9).
As we decrease N, the waiting time increases exponentially, and it eventually becomes impossible to escape fromthe locally stable unimodal state in a reasonable time (Figure8). This implies that evolutionary branching is unlikely tooccur in a small population.
When N is very large, genetic drift can be neglected
and the system is approximated by the ODE: dm2 =dt ¼
w2 m2 2 þ ms2 . Thus, the waiting time can be approximatedby the time required to move from m2 ¼ Nms2 to the rightboundary m2 ¼ r. This time is 1500 for N = 1000 and960 for N = 2000. This calculation, in conjunction with

Figure 7 The stationary distribution of the trait variance m2 . Curves
represent an analytic result (Equation 14), while bars are a histogram ofrealized trait variances in a simulation run for 106 time steps. The dashedline represents Nms2 = 0.0008. In the simulation, the initial trait distribution was set as monomorphic at z = 0.2. The transition from this initialstate to the stationary state was negligible as it took only 5000 time
steps. Parameters: N = 200, m 5 0:01; s 5 0:02; ðb1 ; b2 ; c1 ; c2 Þ 5
ð6:0; 21:4; 4:56; 21:6Þ.

Figure 8, shows that the stochastic ﬂuctuation plays a dominant role even in a population of 1000 individuals.
Our reduced model measures genetic diversity in trait
distribution. An alternative measure of genetic diversity is
the number of segregating alleles, which depends only on m
and is predicted by Ewens’ sampling formula when selection
is absent (Ewens 2004). Figure 10 shows the effect of mutation rate m and mutation step size s2 . Simulations conﬁrmed our prediction that only the value of ms2 matters,even though the numbers of segregating alleles in the simulation with N = 1800 were 50 and 150 for m ¼0.0025
and 0.01, respectively. It is also conﬁrmed that smaller mutational effects require a longer waiting time until evolutionary branching.

Discussion
In this article, we have studied the effect of stochastic
ﬂuctuation due to ﬁnite population size on evolutionary
branching. We have derived an ordinary differential equation and then a stochastic differential equation for the dynamics of trait variance. The analysis has revealed that there
are two different manners in which evolutionary branchingoccurs: deterministic branching and stochastic branching. In
the case of deterministic branching, evolutionary branching occurs relatively soon after the mean trait reaches thebranching point. This occurs when disruptive selection in2
2
tensity w2 exceeds a threshold value, w2 .
p1=ð4ms
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ðN21Þ Þ,which can be approximated by 1=N , 2 w2 ms2. This condition for deterministic branching implies that the intensityof random genetic drift 1=N is smaller than disruptive selection intensity w2 and the mutational effects ms2 (see later inarticle). If the opposite inequality holds, i.e., when disruptive
selection is weaker than the threshold, then a unimodal distribution is locally stable in the deterministic model and

Evolutionary Branching in a Finite Population

235

Figure 9 The distribution of waiting times. Curves represent an analytic
result (see Appendix E), while bars represent a histogram of the waitingtime obtained in 1000 different runs of individual-based simulation.
Parameters: N = 1000, m 5 0:01; s 5 0:02; ðb1 ; b2 ; c1 ; c2 Þ 5
ð6:0; 21:4; 4:56; 21:6Þ.
Figure 8 A Log plot of waiting times for evolutionary branching fordifferent population sizes (N). The thick and thin curves represent analytic
results based on a SDE model (Equation 16) and those based on an ODE
model (Equation 11), respectively. Each open circle represents the resultof a single run in individual-based simulation. The ﬁrst time the trait
variance exceeded 0.1 is shown. In the simulation, the initial trait distribution was set as monomorphic at z = 0.2 and each run consisted of 106
time steps. Simulation runs in which the variance did not explode areplotted at 106 . Time required to reach a singular point in the simulation
was 5000 when N = 200 and 2000 when N = 1800, which is included
to show open circles. Thus, open circles are expected to lie slightly
above the curve. Parameters: m 5 0:01; s 5 0:02; ðb1 ; b2 ; c1 ; c2 Þ 5
ð6:0; 21:4; 4:56; 21:6Þ.

stochastic ﬂuctuation is necessary for the population toundergo evolutionary branching. We call this latter typeof evolutionary branching “stochastic branching.” We observeda relatively long waiting time until evolutionary branchingoccurs in simulations. The distribution and expectation of
the waiting time were mathematically calculated and thesepredictions were consistent with computer simulations of theindividual-based model. Considering the several assumptions
we have adopted, the consistency is better than expected,
although not perfect. As population size or mutational effectdecreases, the waiting time exponentially increases. Mathematically, this is because these factors appear in the exponent
in the expression of the expected waiting time (Equation
16). Although the increase of the waiting time is a continuous
phenomenon resulting from the underlying stochasticity
in individual-based simulations, we can understand evolutionary branching by classifying it as either deterministic branching in which stochastic ﬂuctuation is unnecessary
or stochastic branching in which stochastic ﬂuctuation isnecessary.
Our analysis has revealed that evolutionary branching isalmost impossible in some biological situations in which theeffective population size or mutational effect is sufﬁciently
small. As a population size increases, our condition for deterministic evolutionary branching approaches and converges tothe standard branching-point condition. The condition can be
written as 4w2 N 2 ms2 . 1, which can be used as an approximate indicator; for example, a population with size N = 108

236

J. Y. Wakano and Y. Iwasa

with mutational effects ms2 ¼ 10214 is predicted to show
deterministic branching when disruptive selection intensity
is w2 ¼ Oð1Þ, although such individual-based simulation isalmost impossible.
The condition separating deterministic and stochastic
branching is understood intuitively by considering the force
suppressing variance (genetic drift) and the two forces
promoting the variance (disruptive selection and mutation).
The balance between them determines which type ofbranchingoccurs. The condition for deterministic branching,pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1=N , 2 w2 ms2 can be rewritten as ð1=NÞ2 , 4w2  ms2 ,which meansðgenetic driftÞ2 , 4 · ðselection effectÞ · ðmutation effectÞ:This implies that deterministic branching occurs when theintensity of genetic drift is smaller than double the geometric average of the selection and mutation intensities. It is
important to note that the condition is determined by N 2
instead of N.
To mimic the full individual-based model by using thedeterministic and stochastic equations of the trait variance,
we have adopted several simplifying assumptions in deriving
the dynamics: i.e., trait distribution is assumed to be close tonormal, and the trait variance is small enough to justify
series expansion. From simulation results, it is clear thatthe trait distribution ﬁnally becomes multimodal when evolutionary branching occurs. Thus, our analytic results arevalid only for the initial growth of the trait variance where
unimodal distribution (approximated by a normal distribution) becomes broader at the singular point. This initial explosion is followed by highly nonlinear and complex dynamics
that realize a multimodal distribution. Nevertheless, the agreement of our analytic prediction and simulation results suggeststhat our reduced model captures the onset of evolutionary
branching.
Based on individual-based simulations, Claessen et al.
(2007) reported delayed evolutionary branching and proposed two reasons: random genetic drift and extinction ofadaptive dynamics (e.g., Dieckmann and Doebeli 1999;
Sasaki and Dieckmann 2011; Mirrahimi et al. 2012).
We have focused on the dynamic processes in which the
population escapes from the locally stable value of the trait
variance, while keeping its trait distribution as a unimodal
distribution. Conceptually, if we allow all possible trait distributions, we might have two locally stable states, one ofwhich corresponds to a unimodal state and a second that
corresponds to a fully branched state. Stochastic branchingin the present model considers only the case in which asystem located at the unimodal state escapes from this local
optimum by stochastic ﬂuctuation (random genetic drift).
However, a branched state can also evolve back to theunimodal state by genetic drift; thus, a stochastic analysis
including both escaping from and evolving back to theunimodal distribution would give us an even deeper understanding of evolutionary branching.

Acknowledgments

Figure 10 A Log plot of waiting times for evolutionary branching fordifferent population sizes and mutation parameters. The top panel(m 5 0:0025; s 5 0:02) depicts a result for a smaller mutation rate
than in Figure 8, while the bottom panel (m 5 0:01; s 5 0:01) depicts
a result for a smaller mutation step size. As ms2 5 1026 is the same,the analytic prediction curves are the same. The two simulation results
are also undistinguishable. Parameters: ðb1 ; b2 ; c1 ; c2 Þ 5 ð6:0; 21:4;
4:56; 21:6Þ.

incipient branches. In contrast to Claessen et al. (2007), ourmodel has no ﬂuctuation in the total population size. Thus,
the delayed or absent evolutionary branching is caused bydemographic stochasticity—the ﬂuctuation in the number ofoffspring each individual produces. Our simple model has
enabled us to derive the condition for evolutionary branching mathematically, with random genetic drift taken into
consideration, which has been only verbally proposed before. In addition, our mathematical analysis and simulation
results have revealed that the mutation rate and step size
are as important as population size.
In this article, to illustrate the effect of genetic drift on
the evolutionary branching, we used as an example a pairwise interaction game that incorporated a Wright–Fisher
update. We conjecture that the extension of our model toinclude a Moran update or to a multiplayer game is readily
possible. In the case of including a Moran update, the loss ofgenetic diversity (Equation 7) should be appropriately modiﬁed and also the unit of time should be rescaled since N
Moran steps correspond to a single Wright–Fisher step. The
result of our model should be easily extended to resourcecompetition models, which have been intensively studied in

We thank the following people for their very helpful comments: C. Hauert, Y. Kobayashi, C. Miura, W. Nakahashi,
H. Ohtsuki, and A. Sasaki. This work was supported by Japan
Science and Technology Precursory Research for Embryonic
Science and Technology (J.Y.W.), by a Grant-in-Aid for
Scientiﬁc Research (B) from Japan Society for Promotionof Science (to Y.I.), and also by the support of the Environment
Research and Technology Development Fund (S9) of the
Ministry of the Environment, Japan.

Literature Cited
Barton, N. H., and M. Turelli, 1991 Natural and sexual selection
on many loci. Genetics 127: 229–255.
Champagnat, N., R. Ferriere, and S. Meleard, 2006 Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models. Theor. Popul. Biol. 69: 297–321.
Claessen, D., J. Andersson, L. Persson, and A. M. de Roos,
2007 Delayed evolutionary branching in small populations.
Evol. Ecol. Res. 9: 51–69.
Dieckmann, U., and R. Law, 1996 The dynamical theory of coevolution: a derivation from stochastic ecological processes. J.
Math. Biol. 34: 579–612.
Dieckmann, U., and M. Doebeli, 1999 On the origin of species bysympatric speciation. Nature 400: 354–357.
Doebeli, M., and U. Dieckmann, 2000 Evolutionary branching andsympatric speciation caused by different types of ecological interactions. Am. Nat. 156: S77–S101.
Doebeli, M., C. Hauert, and T. Killingback, 2004 The evolutionary
origin of cooperators and defectors. Science 306: 859–862.
Ewens, W. J., 2004 Mathematical Population Genetics. SpringerVerlag, New York.
Feller, W., 1971 Introduction to Probability Theory and Its Application, Vol. 2, Ed. 3. John Wiley & Sons, New York.
Geritz, S. A. H., J. A. J. Metz, E. Kisdi, and G. Meszena, 1997 Dynamics of adaptation and evolutionary branching. Phys. Rev. Lett.
78: 2024–2027.
Geritz, S. A. H., E. Kisdi, G. Meszena, and J. A. J. Metz,
1998 Evolutionarily singular strategies and the adaptive growthand branching of the evolutionary tree. Evol. Ecol. 12: 35–57.

Evolutionary Branching in a Finite Population

237

Hakoyama, H., and Y. Iwasa, 2000 Extinction risk of a densitydependent population estimated from a time series of population size. J. Theor. Biol. 204: 337–359.
Iwasa, Y., A. Pomiankowski, and S. Nee, 1991 The evolution of
costly mate preferences. 2. The handicap principle. Evolution
45: 1431–1442.
Karlin, S., and H. M. Taylor, 1981 A Second Course in Stochastic
Processes. Academic Press, London.
Kimura, M., 1957 Some problems of stochastic processes in genetics. Ann. Math. Stat. 28: 882–901.
Kimura, M., and T. Ohta, 1969 The average number of generations until ﬁxation of a mutant gene in a ﬁnite population. Genetics 61: 763–771.
Lande, R., 1976 Natural selection and random genetic drift inphenotypic evolution. Evolution 30: 314–334.
Lande, R., 1981 Models of speciation by sexual selection on polygenic traits. Proc. Natl. Acad. Sci. USA 78: 3721–3725.
Lande, R., 1993 Risks of population extinction from demographic
and environmental stochasticity and random catastrophes. Am.
Nat. 142: 911–927.
Metz, J. A. J., R. M. Nisbet, and S. A. H. Geritz, 1992 How should
we deﬁne “ﬁtness” for general ecological scenarios? Trends
Ecol. Evol. 7: 198–202.
Metz, J. A. J., S. A. H. Geritz, G. Meszena, F. J. A. Jacobs, and J. S.
van Heerwaarden, 1996 Adaptive dynamics, a geometrical
study of the consequences of nearly faithful reproduction,

pp. 183–231 in Stochastic and Spatial Structures of Dynamical
Systems, edited by S. J. van Strien and S. M. Verduyn Lunel.
North-Holland Publishing, Amsterdam.
Mirrahimi, S., B. Perthame, and J. Y. Wakano, 2012 Evolution of
species trait through resource competition. J. Math. Biol. 64:
1189–1223.
Pomiankowski, A., Y. Iwasa, and S. Nee, 1991 The evolution of
costly mate preferences. 1. Fisher and biased mutation. Evolution 45: 1422–1430.
Ricciardi, L. M., and C. E. Smith, 1977 Diffusion processes andrelated topics in biology, in Lecture Notes in Biomathematics, Vol.
14. Springer-Verlag, Berlin.
Sasaki, A., and U. Dieckmann, 2011 Oligomorphic dynamics foranalyzing the quantitative genetics of adaptive speciation. J.
Math. Biol. 63: 601–635.
Tazzyman, S. J., and Y. Iwasa, 2000 Sexual selection can increase
the effect of random genetic drift: a quantitative genetic modelof polymorphism in Oophaga pumilio, the strawberry poisondart frog. Evolution 64: 1719–1728.
Wakano, J. Y., and L. Lehmann, 2012 Evolutionary and convergence stability for continuous phenotypes in ﬁnite populations
derived from two-allele models. J. Theor. Biol. 310: 206–215.
Waxman, D., and S. Gavrilets, 2005 Twenty questions on adaptive dynamics. J. Evol. Biol. 18: 1139–1154.
Communicating editor: W. Stephan

Appendices
Appendix A: Fitness Function and the Coefﬁcients
Because individuals are playing a game, their payoff and thus their ﬁtness depend not only on an individual’s own trait but also on
their opponents’ traits. In principle, it is impossible to determine a ﬁtness function for an individual with trait z by an average trait
z because more information is needed [i.e., the distribution of the trait uðzÞ]. However, when the trait variance Var½uðzÞ ¼ m2 issmall, ﬁtness is mainly determined by the deviation of the trait value from the average, z 2 z and can be approximated bywðzÞ ¼ wðzÞ þ w1 ðz 2 zÞ þw2
ðz2zÞ2 ;2(A1)

where the derivatives w1 and w2 are evaluated at a situation for which all individuals (except self) carry the average trait z.
In other words, Equation A1 represents invasion ﬁtness. The coefﬁcients are functions of z, so that w1 ¼ w1 ðzÞ andw2 ¼ w2 ðzÞ. According to adaptive dynamics, a singular point z* satisﬁes w1 ðz*Þ ¼ 0 and it is convergence stable ifðd=dzÞw1 ðzÞjz¼z* , 0. It is evolutionarily stable if w2 ðz*Þ , 0. When a singular point is convergence stable but not evolutionarily stable, evolutionary branching is known to occur.
Here we mainly focus on a case in which the average trait remains at z ¼ z*, and thus we treat w1 and w2 as if they wereconstant coefﬁcients. They can be explicitly calculated for ﬁnite N, but when N .100, it is well approximated by the value
obtained by an inﬁnite population limit. In this limit, ﬁtness function is denoted aswðzÞ ¼f ðz; z*Þ;f ðz*; z*Þ

(A2)and thus
w1 ¼




1
@
f ðz; z*Þf ðz*; z*Þ @z
z¼z*

(A3a)

w2 ¼


 2

1
@

f
ðz;
z*Þ

2
f ðz*; z*Þ @z
z¼z*

(A3b)

and

238

J. Y. Wakano and Y. Iwasa

hold. To calculate the intensity of disruptive selection, not only the sign but also the magnitude is important. For the modelof Doebeli et al. (2004), this becomesw2 ¼

2ðb2 2 c2 Þ;f ðz*; z*Þ

(A4)where z* ¼ ðc1 2 b1 Þ=ð4b2 2 2c2 Þ. We numerically conﬁrmed that the value of w2 can be well approximated by w2  0:1 forthe parameter values ðb1 ; b2 ; c1 ; c2 Þ ¼ ð6:0; 21:4; 4:56; 21:6Þ and for the population sizes N .100. The magnitude of
disruptive selection depends on b1 and c1 , as well as on b2 and c2 . In the vicinity of our choice of parameter values usedin all ﬁgures, disruptive selection is stronger for smaller values of b1 and larger values of c1 .

Appendix B: Derivation of Deterministic Approximation
When both selection and mutation are weak, their effects can be approximated as additive and we obtain the replicator
Equation 4. By integrating Equation 5, which is the Taylor expansion of w(z) aroundZz [ zuðzÞdz;(B1)
we haveZw
 [ uðzÞwðzÞdz ¼ wðzÞ þ 2 m2 ;
w
2

(B2)

where
Z
mk [ ðz2zÞk uðzÞdz

(B3)

is the kth centralized moment around z. Putting this into Equation 4, we haveih w
w2
ms2 @ 22
m2 þ w1 ðz 2 zÞ þ ðz2zÞ2 uðz; tÞ þxðzÞ ¼ uðzÞþ 2
uðzÞ:
2
2
2 @z2
Imposing zero-ﬂux boundary conditions @u=@z ¼ 0, we haveRxðzÞdz ¼ 1
Rw2
m3
xðzÞðz 2 zÞdz ¼ w1 m2 þ
2
R
w2
w2
m2 m2 þ w1 m3 þ
m4 þ ms2 :xðzÞðz2zÞ2 dz ¼ 1 2
2
2
We denote statistical values of the offspring distribution byR
Avx ½z [ xðzÞzdz
R
Varx ½z [ xðzÞðz2Avx ½zÞ2 dz:

(B4)

(B5)

(B6)

Then we have Avx ½z ¼ z þ w1 m2 þ ðw2 =2Þm3 , which is Equation 6a.
To calculate the variance in the offspring pool, we use
R
2
Varx ½z ¼ xðzÞ z2zÞ2ðAvx ½z2zÞ dz
R
w2
2¼ xðzÞðz2zÞ dz 2 ðw1 m2 þ
m3 Þ2
2

(B7)

w2 2
w2
w2
m2 þ w1 m3 þ
m4 þ ms2 2 ðw1 m2 þ
m3 Þ2 :
2
2
2

(B8)

and we have
Varx ½z ¼ m2 2When we neglect the ﬁfth and higher-order terms, we have Varx ½z ¼ m2 2 ðw1 2 þ w2 =2Þm2 2 þ w1 m3 þ ðw2 =2Þm4 þ ms2 ,which is Equation 6b.

Evolutionary Branching in a Finite Population

239

Appendix C: Derivation of SDE Approximation
We want to approximate the dynamics by an SDE with the form
dM2 ¼ f ðM2 Þdt þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
gðM2 ÞdB;(C1)

where dB denotes Brownian motion and M2 denotes a random variable representing the variance in trait distribution. To
obtain gðxÞ, we use an established result in statistics that a random variable

X¼

N 
X
Zi 2z

2

(C2)

v

i¼1

obeys the chi-square distribution where Zi ’s are N independent random variables each of which is sampled from a normal
distribution with mean z and variance v2 . The variance of the chi-square distribution is 2N. When N is large, we can The trait variance M2 can be written byapproximate z in Equation C2 by Z.
M2 [

N
2
1 X
 2 ¼ v X;
Zi 2Z
N i¼1
N

(C3)

where M2 is a random variable, due to random genetic drift. Its variance over probability distribution is given by
Var½M2  ¼ 2N

 2
v
N

2

¼

2v4
:
N

(C4)

According to our model assumptions, we approximate the variance v2 by the trait variance in offspring distribution xðzÞ.
Then, we have the approximation
2m2 2
;
N

(C5)

2ðN 2 1Þm2 2
;
N2

(C6)

Var½M2  ¼
which yields gðxÞ ¼ 2x2 =N and Equation 12.
By a more detailed calculation, we can show
Var½M2  ¼~ðxÞ ¼ x 2ð3N22Þ=ðN21Þ eNaðxÞ , but for our purpose Equation C5 is precise enough.
which results in p

Appendix D: Properties of Waiting Times
Here, we show that the expectation of waiting time
T¼N

Z rZ yhx23 yeðN21ÞðaðxÞ2aðyÞÞ dxdy

(D1)

l

is a decreasing function of w2 , N, and ms2 . First, T is an increasing function of aðxÞ 2 aðyÞ, which satisﬁesaðxÞ 2 aðyÞ ¼ 2 w2 ðy 2 xÞ 2 ms

2

1 1
2 :
x y

(D2)

Since y . x holds in the region of integration in Equation D1, the waiting time T is a decreasing function of w2 and ms2 . To
show the dependence on N, we use
d
T¼
dN

240

J. Y. Wakano and Y. Iwasa

Z rZ yhlx 23 eðN21ÞðaðxÞ2aðyÞÞ Bðx; yÞdxdy;(D3a)whereBðx; yÞ ¼ y þ Nðx 2 yÞðwxy þ ms2 Þ
:
x

(D3b)

When N is large, B(x, y) has the same sign as x 2 y, which is negative. Thus, T is a decreasing function of N for large N. In
addition, since aðxÞ 2 aðyÞ is a negative constant independent of N, we havelimN/N T ¼ 0;

(D4)with all the other parameters ﬁxed at positive values. This is true even for the vanishing mutation rate (i.e., Nms2 ¼ const:
when N/N). However, it does not necessarily go to zero if NðaðxÞ 2 aðyÞÞ does not go to 2N. For example, if Nw2 andNms2 remain constant as N/N, the waiting time T can be an increasing function of N.

Appendix E: Distribution of Waiting Times
Let cðxÞ be an indicator function such that cðrÞ ¼ 1 and cðxÞ ¼ 0 for l # x , r. Then a random variable cðM2 ðtÞÞ is unity if
absorption at m2 ¼ r has occurred until time t, and it is zero otherwise. Let uðh; tÞ ¼ Eh ½cðM2 ðtÞÞ, where Eh denotes the
expectation in the stochastic process with the initial value M2 ð0Þ ¼ h. Then the following Kolmogorov backward equation
holds (Karlin and Taylor 1981):
@u
@u 1
@2u
¼ f ðhÞ þ gðhÞ 2 :
@t
@h 2
@h

(E1)

Let Q be a random variable representing the ﬁrst passage time; then it is intuitively clear that
uðt; hÞ ¼ Pr½Q # tjm2 ð0Þ ¼ h(E2)is the probability that absorption at m2 ¼ r has occurred on or before time t. The probability density that absorption occurs
exactly at time t is given by @u=@t. A similar logic and calculation were adopted in a classical article studying the ﬁxation
probability in population genetics (Kimura 1957). The expected ﬁrst passage time is
Z N
@u
t dt:
TðhÞ ¼(E3)@t
0
After some calculations, we obtain
21 ¼ f ðhÞdh2 d2TðhÞ þ
TðhÞ;
dh
N dh2

which is Equation 15. To obtain the distribution of waiting time, we numerically solve the partial differential Equation E1
with domain h 2 ½l; r. We impose a reﬂecting boundary at left and Dirichlet boundary at right, uðt; rÞ ¼ 1. In the numericalalgorithm, the initial value is set as uð0; hÞ ¼ 0 except uð0; rÞ ¼ 1. The Forward-Time Central-Space method with small dt~ðt; hÞ. To check the validity of this result, we numerically put it into Equation E3
and dh values gives the numerical solution u
and compared the result with that of the numerically obtained value of Equation 16. The expected waiting times obtained by
the two different approaches were almost the same. The distribution of the waiting time t from the initial trait variance Nms2
~ðt; Nms2 Þ, which is shown in Figure 9.
is given by u

Evolutionary Branching in a Finite Population

241