Evolutionary branching in a stochastic population model with discrete mutational steps

# Evolutionary branching in a stochastic population model with discrete mutational steps

## Abstract

Evolutionary branching is analysed in a stochastic, individual-based population model under mutation and selection. In such models, the common assumption is that individual reproduction and life career are characterised by values of a trait, and also by population sizes, and that mutations lead to small changes in trait value. Then, traditionally, the evolutionary dynamics is studied in the limit . In the present approach, small but non-negligible mutational steps are considered. By means of theoretical analysis in the limit of infinitely large populations, as well as computer simulations, we demonstrate how discrete mutational steps affect the patterns of evolutionary branching. We also argue that the average time to the first branching depends in a sensitive way on both mutational step size and population size.

###### keywords:
Evolutionary branching, genetic drift, selection, adaptation
1\cortext

[cor1]Corresponding author

## 1 Introduction

The development of populations subject to frequency- or density-dependent selection is an important and central but not easily analysed topic in theoretical biology. Many different models and approximation schemes have been developed in order to understand how trait values change and how the intriguing possibility of ‘evolutionary branching’ may arise, where a population of individuals with a well-defined trait value or genotype (monomorphic, in other words) experiences a split leading to two or more coexisting trait values or genotypes in the population. This process is important because it gives rise to increased biodiversity.

A variety of population-genetic and game-theoretic methods, as well as the body of ideas and approximation schemes known as adaptive dynamics, cf. Metz et al. (1996); Geritz et al. (1998); Diekman (2004); Waxman and Gavrilets (2005); Metz et al. (2012), have been employed in order to study evolutionary branching. Several simplifications are commonly resorted to. First, complications due to mating and recombination in diploid populations are disregarded by ascribing to each individual a single trait value that is transferred by clonal reproduction unless mutation occurs (see however the recent interesting approach to adaptive dynamics with Mendelian inheritance by Collet et al., 2011). Selection is taken to act through a trait- and density-dependent fitness function, identified with mean reproduction. Second, it is assumed that fitness is a smooth function of the trait value and the population density. The third postulate is that the mutation rate is so small that only few trait values are represented in the population at any one time. Mutations thus occur so rarely that we can talk of a separation of time scales between long-term evolutionary dynamics (the sequence of mutations) and short-term population dynamics, defined by birth and death events, an evolutionary versus an ecological time scale. Fourth, it is assumed that mutations lead only to small changes in the trait value (that mutational step sizes are small). Fifth, the population size is taken to be large, and effects of random genetic drift are neglected (apart from incorporating the fixation probability of an advantageous mutation).

In the adaptive-dynamics literature a further step is taken by letting , , and . Then evolution of a monomorphic population to the first branching of traits is analytically described by the so-called ‘canonical equation’ (Dieckmann and Law, 1996), or more generally by deterministic evolution on the fitness landscape (Waxman and Gavrilets, 2005a; Lande, 1976; Iwasa et al., 1991). How these deterministic adaptive dynamics may arise from a strict stochastic population model has been described in detail in the elegant treatment of Champagnat (2006) and later papers.

As a monomorphic (single-trait) population approaches a local fitness maximum where evolutionary branching may occur in the sense that populations with two different trait values may coexist from there onwards, the canonical equation fails to describe the dynamics. It becomes necessary to ask how stochasticity in a finite population affects the possibility of evolutionary branching. Further questions are: what is the shape of the bifurcation diagram of trait values in the wake of a branching? How do the values of and influence the evolutionary dynamics in a finite population? How long does it take (in the mean) for evolutionary branching to occur?

We address these matters analytically and by computer simulations of a stochastic, individual-based model for the evolution of trait values in a finite population subject to density-dependent fitness function symmetric with respect to its maximal value at the trait . In the limit of , , and small, we show that the evolution of a monomorphic population towards the fitness optimum at is governed by a canonical equation. For our model, standard stability analysis in the limit of would predict evolutionary branching of a monomorphic population at . We demonstrate that at small but fixed mutational step sizes , evolutionary branching typically occurs at non-zero values of . In the limit of infinite population size and vanishing mutation rate we compute the critical value of where the first branching is most likely to occur, as a function of . In addition we show that the evolutionary dynamics in this limit gives rise to a non-symmetric bifurcation diagram after the first branching (for a symmetric model such as ours, the standard theory predicts a symmetric bifurcation diagram in the limit ). Moreover, we state geometrical conditions determining the trait values where the second evolutionary branching event is most likely to occur.

Finding the average time until evolutionary branching occurs in our model requires a more refined analysis. It is necessary to specify how and . This determines the stochastic population dynamics in the critical stage preceding the first branching. That stage must last long enough to render branching possible. We discuss results of computer simulations relating to this question briefly, but a precise mathematical description of the problem is left for future work.

The remainder of this paper is organised as follows. In Section 2 we introduce the stochastic, individual-based model used in simulations. Section 3 summarises standard results obtained in the limit , , and . These are contrasted with numerical experiments at relevant values of , , and in Section 4. Discrepancies between the predictions and the results of numerical experiments are discussed. Our results (and the calculations supporting them) addressing the issues raised in Section 4 are summarised in Section 5. Section 6 contains our conclusions. Three appendices summarise details of our calculations.

## 2 Model

The stochastic, individual-based model used in the numerical experiments discussed below describes a finite population of individuals competing for resources. Each individual is assigned a trait value . Individuals with the same trait value (belonging to the same ecological ‘niche’) compete for a limited amount of resources. The scarcity of resources is parametrised by a trait-dependent carrying capacity . There is an ‘optimal’ trait value corresponding to the value of where assumes its global maximum. Individuals belonging to different niches (corresponding to different trait values, and , say) may also interact. The strength of interaction depends upon the difference in trait values.

Mutations give rise to changes in the trait value, . Here is a fixed mutational step size and the set of possible traits occurring in the population is . Individual fitness is measured by the mean offspring number, which is determined by the trait values and the intensity of competition. The latter, in its turn, is controlled by the scarcity of resources and the number of individuals. In other words, the model describes the evolution of a trait value, or of trait values, in a population subject to density-dependent selection. The population would be expected to evolve towards the optimal trait value. Intense competition between individuals with trait values close to the optimal one may however give rise to sub-populations with well-defined but distinct trait values. This phenomenon is referred to as evolutionary branching, and is the subject of this paper.

The number of individuals with a given trait value is denoted by . It is assumed that each individual has either none or two offspring (Klebaner et al., 2011). The offspring of individuals in one generation constitute the next generation. We denote the fitness function by which is the mean offspring number for an individual with trait value , and where is the set of trait values represented in the population. This fitness is a function of the co-existing sub-population sizes (. We take it to be of the form:

 Mx(Zy,y∈D)=21+(NKx)−1∑y∈DγxyZy. (1)

Here the carrying capacity of the subpopulation with trait value is , where is a population-size parameter to be thought of as large. Competition between subpopulations corresponding to different trait values and is regulated by the function .

The offspring usually inherit the parental trait , but with a small probability they mutate (independently) either to or to , with equal probabilities. Since time is measured in generations, the model is thus discrete both in time and trait space. It can be viewed as the simplest possible structure (‘the bare bones’, Klebaner et al., 2011) behind adaptive dynamics, and as such is an interesting mathematical object in its own right. Since all individuals have the same life span, one, it displays an unyielding age dependence which can be seen as a diametrical counterpart and natural complement to the equally artificial non-aging individuals of birth-and-death processes, underlying the sequence of papers by Champagnat and Méléard (2011) and others, and also classical deterministic approaches. For more realistically age-structured models cf. Méléard and Tran (2009) as well as Jagers and Klebaner (2011).

The form (1) is certainly ad hoc, and should be viewed as an explorative illustration. But it is not without historical precedence (Christiansen and Loeschke, 1980). As there, we take and to be Gaussian functions. Indeed, this is a choice often encountered in the adaptive-dynamics literature,

 Kx=e−x2andγxy=e−(x−y)2(1+α). (2)

The parameter regulates the competition between subpopulations with different traits. A larger value of implies weaker competition. We take . This ensures that competition between subpopulations is weak enough to allow for evolutionary branchings. It is convenient to replace population sizes by densities , so that the fitness function (1) takes the form

 Mx(fy,y∈D)=21+∑y∈DRxyfy. (3)

Here

 Rxy=γxyKy/Kx=e−(x−y)(αx−(2+α)y) (4)

is the fitness kernel which combines the effects of the interaction and the carrying capacity .

Eqs. (3) and (4) allow for explicit formulae for the equilibrium densities at which subpopulations may coexist in the absence of mutations. Note however that even this coexistence is temporary, albeit potentially of long duration (Jagers and Klebaner, 2011). We shall still refer to them as equilibrium densities. They are defined by the condition , where

 MDz≡21+∑y∈DRzyfDy. (5)

The equation describes the critical regime of reproduction. It corresponds to the condition

 ∑y∈DRxyfDy=1, x∈D. (6)

In the monomorphic case we find that and in the absence of mutations, the density stabilises at the value . Substituting this result into Eq. (5) we obtain

 M{x}z=21+e−(z−x)(αz−(2+α)x). (7)

Importantly, the function given by Eq. (5) plays a crucial role in the invasion analysis. Suppose an individual of type appears in a stable population with the set of traits as a mutant of . Its fitness has the form

 Mz(fDy,y∈D;f∗z) =21+∑y∈DRzyfDy+f∗z≈MDz

because and the small initial frequency of the mutant can be neglected. If , there is a chance for the new trait to establish itself in the set of coexisting traits. If , the mutation can not establish itself and is wiped out by genetic drift.

In summary, we have chosen an extremely simple model which however retains so much of the structure of evolutionary dynamics that the questions raised in the introduction can still be formulated. It is described by only four parameters, the population-size parameter , the mutation rate , the mutational step size , and the interaction parameter .

## 3 Standard predictions in the limit ϵ→0

In spite of its simplicity, the model introduced in the previous Section is not easily analysed. As pointed out, the standard adaptive-dynamics approach (Geritz et al., 1998) studies continuous-time population models in the limit and , and then . In this Section we state the corresponding results for our discrete-time model.

If can be thought of as infinite, law-of-large-number effects replace random processes by their expectations, and if both and , then only a small number of trait values are represented in the population at any specific time: the basic cases studied below are those of monomorphic populations, , and dimorphic ones, .

In this limit, a time-scale separation appears between the evolutionary dynamics (referring to long-term changes of the trait value) and the ecological time scale of population dynamics. Letting makes it possible to describe the dynamics of the trait value by local stability analysis, central to adaptive dynamics.

Fig. 1 illustrates the standard predictions obtained in these limits. First, starting from an initially positive trait value, is expected to evolve deterministically towards the local optimum (at in our model). The path is described by a differential equation referred to as a ‘canonical equation’ (Dieckmann and Law, 1996; Champagnat et al., 2001) and related to similar equations studied in population genetics (Lande, 1976; Iwasa et al., 1991). Second, the first evolutionary branching occurs at . Third, after the first branching, the pair of traits evolves deterministically in a symmetric fashion, and the second and third branchings occur simultaneously.

### 3.1 Deterministic evolution towards the first branching

Assume that the population starts from a positive initial trait value which is large compared to the small value of the mutational step size : . Classical adaptive dynamics then predicts that the initial trait value is consecutively replaced by smaller values of , driving towards . The resulting function of time is described by a differential equation for , known as the canonical equation (of trait evolution), the standard reference being Dieckmann and Law (1996). In the notation of Champagnat et al. (2001, see Eq. (7a) ) the classical canonical equation of trait evolution has the form

 dxdt=μ(x)σ202n(x)∂1f(x,x). (8)

In this expression, is the mutation rate for the trait value , is the mutational variance, is the equilibrium population size, and is the fitness function. Recall that within a model formulated in terms of stochastic birth and death processes, the fitness function is the difference between the birth and death rates for rare mutants with trait in an -population at its quasi-stable size.

In our case, the mutation rate is independent of the trait value , is simply , and . The fitness function corresponding to our discrete-time model should be calculated as

 f(z,x)=logM{x}z, (9)

where is given by Eq. (7). This formula relies on a well-known property of a linear birth-death process stemming from a single individual: if is the difference between the birth and death rates per individual, the population size at time has mean . Eq. (9) stipulates that the expected population size after one unit of continuous time is equal to the expected size of the first generation in our discrete-time model.

From Eqs. (9) and (7) we conclude that the fitness gradient is given by

 ∂1f(x,x)=−x,

since . In our case we would thus expect the canonical equation to take the form

 dxdt=−μϵ22Ne−x2x. (10)

However, the correct canonical equation for our model Eq. (19) derived in Section 5.1 reveals that the trait substitution process goes twice as fast in the discrete-time model compared to the continuous-time model. This is reminiscent of the well-known phenomenon of genetic drift running twice as fast in the Moran model than in the Wright-Fisher model of the same size, see for example Wakeley (2008), though the latter phenomenon has a different explanation.

### 3.2 Location of the first evolutionary branching

The canonical equation suggests that the substitution process of the trait value slows down as the point is approached. In the limit of , , and , standard adaptive-dynamics analysis predicts the first evolutionary branching to occur at , giving rise to two branches, and .

### 3.3 Second and third evolutionary branchings

Moreover, in the very same limit, the deterministic evolution after the first branching is predicted to be symmetric. In other words, the two trait values and after the first branching are expected to evolve as . The second and third evolutionary branchings are expected to occur simultaneously when reaches a critical value, which in our case can be computed explicitly as

 xα=12√log(2α+1)1+α. (11)

As a function of the value of reaches its maximum at .

## 4 Computer simulations

Figs. 2, 3 and 4 summarise results of direct numerical simulations of the model described in Section 2 for large values of , small values of , and for a small value of (equal to ). In the algorithm, binomial distributions of numbers of mutants have been approximated by the appropriate Poisson distributions.

We now compare simulation results with the predictions from Section 3. Fig. 2a shows ten realisations of the evolution of for the mutational step size from an initially large value . (The additional parameters are given in the figure caption). We see that decreases towards the first branching (which in all cases occurs for positive values of and not at ). Fig. 2b shows a plot of the trait value versus the average time of staying at , conditional on no branching having occurred. Also shown is an approximate solution of the correct canonical equation derived in Section 5.1 as Eq. (19). This is expressed in terms of the expected time for to reach level :

 t(x)=x0∑y=x−ϵλ−1y, (12)

where (see below) is the rate of the asymptotically exponential holding time at level . For large values of we observe good agreement between (12) and the average of ten independent realisations of the substitution process. But as is approached, the numerical average falls below the predicted average. The neighbourhood of the first branching event is not described by the canonical equation.

Fig. 3 shows results of two computer simulations starting in the vicinity of . In Fig. 3a the first branching does not occur at but at . Further, evolution after the first branching is not symmetrical. This can also be seen from Fig. 4 which displays the data of Fig. 3a in the - plane (the standard adaptive-dynamics approach predicts, as pointed out above, that the pair () moves on the diagonal shown in Fig. 4 as a dashed line). Note also that the second and third branchings do not occur simultaneously, and not symmetrically in the simulations.

Fig. 3b shows results of a direct numerical simulation with the same overall substitution rate as Fig. 3a, but for a smaller value of and a larger value of . In this example, branching did not occur during the simulation time . By contrast, the trait value is seen to diffuse around .

In summary we find that the standard adaptive-dynamics approach describes the initial time evolution of the trait value well, provided is large enough. In the vicinity of and after the first branching, however, the observed patterns of evolutionary dynamics (Figs. 3 and 4) differ from the standard predictions. In the following we show that the asymmetric patterns observed in Figs. 3a and 4 are due to the positive mutational step size . In the limit and we were able to characterise the branching patterns at small but fixed values of . To determine how long it takes on average until evolutionary branching occurs in a finite population requires a more refined analysis. We return to this in Section 5.3.

## 5 Results

### 5.1 Deterministic evolution towards the first branching

In the monomorphic case () Eq. (7) implies

 M{x}z−1≈x(x−z)+(α/2−x2)(z−x)2, (13)

for . For it follows that

 M{x}x−ϵ>1,M{x}x+ϵ<1,M{x−ϵ}x<1,andM{x+ϵ}x>1. (14)

These inequalities show that the expected trend of adaptive evolution in the monomorphic phase is simple: the initial positive trait value is consecutively replaced by smaller values driving towards zero. The dynamics of this trait substitution process is described by the canonical equation whose correct form for our model is derived next.

Observe that the expected number of mutations in the -population, assumed to be at its quasi-stable size, , is per generation, where is divided by 2 since we discard the -mutations. Thus the expected waiting time between two consecutive replacements is

 λx≈μ2Ne−x2Px−ϵ(x). (15)

Here stands for the probability that a population from one -individual survives and takes over from the monomorphic -population. To find the probability of the mutant invading and replacing a resident population with trait we can approximate the population dynamics of the invader by a binary branching process starting from a single individual. The survival probability of this branching process is

 Pz(x)=2M{x}z−1M{x}z, (16)

see for example Eq. (5.66) in Haccou et al. (2005). Eq. (13) implies

 M{x}x−ϵ≈1+ϵx, (17)

and we obtain , at least for far from zero. (The corresponding classical result for the Wright-Fisher model is the following. Suppose that an advantageous allele is introduced in the population with its mean offspring number being larger compared to the wild type by a factor . Then the fixation probability of this allele is approximately given by for small values of .)

Substituting into Eq. (15) results in . Now, since the evolution of the trait value as a function of time is approximately given by

 −dtdx≈t(x−ϵ)−t(x)ϵ≈1ϵλx, (18)

we conclude that the canonical equation takes the form

 dxdt=−μϵ2Ne−x2x (19)

which is similar (but not identical) to Eq. (10). Indeed, Eqs. (19) and (10) differ by a factor of .

This factor can be explained as follows. The standard derivation of Eq. (8) (Dieckmann and Law, 1996; Champagnat et al., 2001) refers to birth-and-death processes in continuous time (or corresponding deterministic formulations), whereas our model considers binary splitting in discrete, non-overlapping generations. Our form of the equation appears explicitly for Poisson reproduction in discrete time in Metz (2011). The fact that different effective population sizes can appear in the canonical equation was first noted by Durinx, Metz, and Meszéna (2008).

Note that the speed of the trait-substitution process given by Eq. (8) would not change if we were to take different values of the birth and death rates, as long as their difference remains the same. We argue that the proper choice of the birth rate in the monomorphic stable population is unity, matching the birth rate in our discrete-time model: one birth per generation per individual, . Then, in accordance with Eq. (9), the corresponding survival probability of a single mutant is given by . Comparing this with its discrete-time counterpart , we conclude that the survival probability in the discrete model is two times larger.

An important assumption behind Eq. (19) is that successful mutations should not come too soon after each other in order to make sure that selective sweeps do not overlap (‘clonal interference’), see Su-Chan Park et al. (2010). As a very rough approximation, we estimate the average fixation time from the equation

 (M{x}x−ϵ)Tfix≈N. (20)

This equation simply suggests to neglect the competition among individuals and equate the target population size with the mutant reproduction factor for consecutive generations. Combining Eqs. (17) and (20) we find in the absence of competing mutations. The average time between successful mutations is . So the condition ensuring non-overlapping sweeps is that the product must be very small so that

 μNlogN≪1. (21)

In other words, condition Eq. (21) guarantees that after a mutation the new set of types settles in an equilibrium before the next successful mutation event.

### 5.2 Location of the first evolutionary branching

The trait substitution process considerably slows down as approaches zero, so that it must be considered as a sequence of discrete jumps of size towards zero, rather than a continuous process. At levels close to zero a more refined analysis of fitness asymptotics is required.

In the following we derive an approximation for the location of the first evolutionary branching. We show that it happens close to zero, but typically not at . More precisely, we demonstrate that evolutionary branching becomes possible in an interval of width of order around . It is thus necessary to distinguish between trait values separated by a few mutation steps. For simplicity we assume in the following that initially , so that the trait substitution process approaches from above (as in the simulation results shown in Fig. 3). The problem of determining the location of the first branching event is symmetric and corresponding expressions for are easily found.

The aim is to find trait values where mutations to and from are mutually invasive. The corresponding condition is: and . Using Eq. (13) and dropping the terms much smaller than we arrive at two conditions:

 M{x}x−ϵ−1 ≈ xϵ+αϵ2/2>0,and M{x−ϵ}x−1 ≈ −(x−ϵ)ϵ+αϵ2/2>0. (22)

We conclude that the first evolutionary branching is possible for , so that populations with trait values and are mutually invasive provided:

 ϵ≤x≤j∗ϵwithj∗=⌈α/2⌉. (23)

Here denotes the smallest integer greater than . (To avoid complications arising in the situation when and for , we will assume that is not an even integer.) The fact that the critical trait value is larger for larger values of is very intuitive: as competition becomes weaker, initially impossible invasion of into becomes favorable once the benefit of reduced competition outweighs the cost of having a ”worse” trait value.

Fig. 5 shows where the first evolutionary branchings occurred in an ensemble of 300 simulations of the stochastic, individual-based model. The parameter was chosen to be which implies that . The parameters were chosen so that the overall substitution rate was the same for all runs, . All simulations were started at , and were run for the same total time , given by . In a small number of runs evolutionary branching did not occur during this time. Fig. 5 shows where the first branching occurred for the remaining runs. First, we see that apart from a small number of cases, branching occurs in the region predicted. Second, the larger the population size is, the more likely it is that evolutionary branching occurs at the boundary of the branching region ( in this case). To explain this behaviour requires a more refined analysis. We return to this question below in Section 5.3. Third, in the limit the condition (23) is consistent with the standard prediction in this limit (that the first branching occurs at ). Fourth, some branchings occurred one mutation step earlier than expected. This can be explained by the fact that occasionally it can take a long time for a slightly advantageous mutant to stabilise, long enough for the next mutation to initiate the first evolutionary branching.

### 5.3 Critical state of coexistence prior to first branching

In the preceding section we have derived conditions determining where the first evolutionary branching may occur, Eq. (23). However, results of our computer simulations of the stochastic model described in Section 2 also demonstrate that while evolutionary branching may occur (and often does occur) when these conditions for evolutionary branching are met, branching need not necessarily take place immediately. This raises the question how long it takes on average in a finite population for evolutionary branching to occur (a question also raised by Boettiger (2010)). What is the probability that evolutionary branching happens when the substitution process reaches the boundary of the region where the first evolutionary branching becomes possible (Fig. 5)? Fig. 6 shows a realisation of a computer simulation of this situation for (corresponding to ). In Fig. 6a we see how the trait value reaches the critical value from above, and how subpopulations with trait values and coexist until evolutionary branching occurs. The corresponding population sizes are shown in Fig. 6b. For evolutionary branching to occur, the critical state of coexistence must last sufficiently long for a favourable mutation to occur. It is possible, however, that the critical state of coexistence may spontaneously disappear by a fluctuation. We see in Fig. 6b that the population-size fluctuations in the coexistence state are substantially larger than the population-size fluctuations in the monomorphic states.

A preliminary analysis shows that while the fluctuations of the monomorphic population sizes are of order , the fluctuations of the subpopulation sizes in the critical state of coexistence are much larger: of order in the limit of large values of and small values of . The smaller value takes, the shorter is the average life time of the critical state of coexistence, lowering the probability that evolutionary branching occurs in this state. Without going into details, we just briefly sketch our argument (which remains to be made rigorous). Let be the densities of two populations coexisting at time and having the neighbouring trait values and with . The evolution of the vector can be approximated by a stochastic dynamical system

 ⎧⎪ ⎪⎨⎪ ⎪⎩fk+1=2fk1+fk+agk+1√Nuk,a=e−ϵ2(α+2−2j),gk+1=2gk1+bfk+gk+1√Nwk,b=e−ϵ2(α+2j), (24)

where and are independent normal random variables with zero means and standard deviations and respectively. Provided the deterministic part has a nontrivial stable point we consider a linearised version of this system around this stable point. Our analysis indicates that while the sum behaves as a stationary autoregression process with fluctuations of the order , the difference behaves as a stationary autoregression process with fluctuations of order . Given that the fluctuations in the sizes of coexisting populations featuring the trait values are of order a rough estimate of the probability of the sudden loss of coexistence is obtained using the approximate tail probability for the normal distribution , where and are positive constants. Thus, the life time of the critical state of coexistence is expected to be inversely proportional to the probability of sudden loss of coexistence so that

 logT∼SαNϵ2 (25)

for some constant depending on through (23). Fig. 7 shows simulation results for the life time as a function of and . The results of the simulations are consistent with the expectation Eq. (25). We find . The precise mathematical derivation of Eq. (25) and the calculation of the constant are interesting questions for further work.

### 5.4 Evolution after the first branching

Condition Eq. (23) demonstrates that if the first evolutionary branching occurs, then it happens within a region of size around . We now turn to the evolution of the pair of trait values after the first evolutionary branching. Standard theory suggests that the dimorphic population evolves symmetrically (along the dashed line in Fig. 4). Our simulations of the stochastic, individual-based population model with discrete mutational steps show that the pair of coexisting trait values followed through the chain of consecutive replacements does not develop in a fully symmetrical fashion. The random path on the --plane (see Fig. 4) follows the trajectory of a random walk with steps from to either or depending on which of the two branches the next replacement has been successful. The repulsion between the two branches is due to the pressure of competition. It is advantageous to stay further apart reducing interspecies competition.

Observe that the consecutive steps in this random walk are weakly negatively dependent: the further the walk deviates from the diagonal , the larger is the probability for the next step to decrease this deviation. To see this, consider a pair of coexisting populations having trait values and with and for definiteness . For a given branch, the corresponding replacement rate is a product of two factors (see the discussion leading to Eq. (15) in the monomorphic case): the stable population size and the probability of fixation for a single mutant. We show that given the lower branch is closer to zero, , both factors of the replacement rate for the lower branch are larger making the move of the lower branch (always off zero, see Appendix B) more probable. Indeed, as the negative trait value is closer to zero, its stable population size must be larger than the size of the population with the positive trait value. On the other hand, as computed in Appendix B, the probability of fixation for a single mutant in the -population is also larger, Eq. (49).

The last observation implies that with high probability the locations of the pair of branches satisfy

 x1=−x2+O(√ϵ) (26)

prior to the second evolutionary branching. This means that the the expected location of the second evolutionary branching, if any, should not deviate from the diagonal more than a constant times . Indeed, the predicted deviation from the diagonal can only be made larger if we neglect the negative dependence and consider a symmetric random walk with independent steps. Referring to the classical central limit theorem for the simple symmetric random walk we observe that the quantity (which is always zero in the symmetric case) is of order , where is the number of jumps in the random walk (representing the number of replacements since the first evolutionary branching). Since the second evolutionary branching takes place at values of and of order unity, one must wait for a large number of jumps, namely of order for the second branching to occur. We conclude that the path taken by the two branches in the - plane will deviate from the diagonal by a distance of at most of order .

### 5.5 ’Triggering cross ’for the second branching

We now discuss the location of the second evolutionary branching. As in the previous section, this location is determined by mutual invasibility analysis, but now for the dimorphic case . The second evolutionary branching may occur on either of the two branches stemming from the first evolutionary branching. We first discuss the condition for the case where the second evolutionary branching occurs on the upper branch of the dimorphic population (an example is shown in Fig. 3a). With the corresponding conditions for mutual invasibility are and . Analysis of the signs of and , reported in Appendix B, shows that in view of Eq. (26) the second evolutionary branching of the upper branch becomes possible in the region

 x1=−xα+O(√ϵ),x2=xα+O(√ϵ), (27)

where is given by Eq. (11). More precisely, we show (see Appendix C) that the second evolutionary branching of the upper branch with high probability occurs in the region

 (x1+xα)cα−x2+xα=O(ϵ) (28)

which is a neighbourhood of the straight line

 (x1+xα)cα=x2−xα. (29)

The coefficient is given by

 cα=(1+2α)2log(1+2α)−2α(1+α)(1+2α)(2α−1)log(1+2α)+2α(1+α). (30)

It is shown in Appendix A that satisfies . The minimum value of equals and is achieved at .

Turning to the possibility that the second branching occurs on the lower branch we find a similar condition

 (x2−xα)cα−x1−xα=O(ϵ) (31)

corresponding to an -neighbourhood of another straight line

 x1+xα=(x2−xα)cα. (32)

Taking these results together, we have shown that the second evolutionary branching is expected to occur in the vicinity of the lines Eqs. (29) and (32). If the pair of diverging branches comes close to the line Eq. (29), then a second evolutionary branching of the upper branch becomes possible. Conversely, when the pair comes close to the line Eq. (32), a second evolutionary branching may occur on the lower branch. Note that since is assumed to be independent of the trait value, our answers Eqs. (29) and (32) depend only on . The cross formed by the pair of lines (29), and (32) is centered at the point in the --plane.

Eqs. (29), (32) are consistent with results of direct numerical simulations of the model. Fig. 8 shows the lines (29) and (32) for and , as well as the locations of second evolutionary branchings for an ensemble of simulations of the stochastic, individual-based model. Six sets of dimorphic initial conditions were used, and . We remark that the last three initial conditions deviate substantially from the diagonal (dashed line in Fig. 4). As shown above, paths in the --plane typically exhibit deviations of order from the diagonal. We have nevertheless included dimorphic initial conditions far from the diagonal in order to clearly exhibit the locations of the second evolutionary branching.

In Fig. 8, blue dots correspond to cases where the second evolutionary branching occurred on the upper branch of the dimorphic state, while red dots correspond to cases where the second evolutionary branching occurred on the lower branch. The coordinates of the dots show the location of the second branching event in the --plane (we determined this location by taking the average of the trait values of the critical state of coexistence prior to the second branching, as well as the coordinate of the second branch when coexistence first occurred). Fig. 8 is consistent with the theoretical expectation that second branchings of the upper branch are triggered by the line Eq. (29), while the second branchings of the lower branch are triggered by the line Eq. (32).

## 6 Conclusions

We have analysed evolutionary branching in a stochastic, individual-based model for the evolution of a trait value subject to small but non-negligible mutational step sizes . We found that the branching patterns are in general asymmetric at distinct, non-null mutational steps, complementing the standard adaptive dynamics predictions Geritz et al. (1998) valid in the limit . In particular we found conditions describing the locations where the first branching may occur, Eq. (23), and where the second branching may take place, Eqs. (27), (28), and (31). Results of simulations of a stochastic, individual-based model at small mutational step sizes were seen to be consistent with these conditions. Our results are derived in the limit of large population sizes and small mutation rates , but allow for fixed mutational step sizes. In this respect our results Eqs. (23), (27), (28), and (31) complement the established approach.

We have also demonstrated that population-size fluctuations in the critical state of coexistence prior to a branching may erase this state. Indeed, its sojourn time depends sensitively on mutational step size. Evolutionary branching occurs typically only when this time is much longer than the time between successful mutations.

## Acknowledgments

Financial support from the Swedish Research Council (SS, BM), the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine (BM), by the Bank of Sweden Tercentenary Foundation (SS) is gratefully acknowledged. VV was partially supported by the grant RFBR 11-01-00139 and the program ”Dynamical systems and control theory” of the Russian Academy of Sciences.

## Appendix A

This appendix is divided into three sections summarising material needed in appendices B and C.

1. We discuss the properties of the three functions:

 C(α)=α2−((1+α)2−α2)x2α, (33) c1(α)=(1+2α)22α(1+α)log(1+2α)−1, (34) c2(α)=(1+2α)(2α−1)2α(1+α)log(1+2α)+1. (35)

Here is defined by Eq. (11). We show that , and for all positive values of . The latter inequalities in view of imply that the constant , defined by Eq. (30) and determining the slopes of the triggering lines Eqs. (29), (32), satisfies .

To see that we observe that

 C(α) =α2−1+2α4(1+α)log(1+2α)=h(1+2α)8(1+α), (36)

where . The function takes positive values for all , since and for . Thus for all positive .

Similarly, the fact that follows from

 c1(α)=h1(1+2α)2α(1+α),c2(α)=c1(α)+4αC(α), (37)

where is positive for (since and ).

2. We discuss properties of the equilibrium densities in the dimorphic case . In this case, the equilibrium densities and satisfying the system Eq. (6):

 Rx1x2f{x1,x2}x2+f{x1,x2}x1=1, f{x1,x2}x2+Rx2x1f{x1,x2}x1=1,

can be represented as and with

 ϕ(x1,x2)=1−Rx2x11−Rx1x2Rx2x1=1−e−(x2−x1)(αx2−(2+α)x1)1−e−2(1+α)(x1−x2)2. (38)

In particular, in the symmetric case and , we obtain

 ϕ(−x,x)=ϕ(x,−x)=11+e−4(1+α)x2. (39)

3. We quote the Taylor expansion for the function defined by Eq. (5) (recall that for )

valid for small values of . Here the terms

assume the form of first and second moments of the coexisting trait values using the weights satisfying condition Eq. (6).

## Appendix B

In this appendix we discuss in which region of the --plane the second evolutionary branching becomes possible. It may occur on either of the two branches stemming from the first evolutionary branching. We start by discussing the evolution of the upper branch of the dimorphic population on its way towards the second branching (an example is shown in Fig. 3a). With the corresponding conditions for mutual invasibility are and . The signs of and are determined by Eq. (40) with

 A{x1,x2}x2 = x2ϕ(x1,x2)+x1Rx2x1ϕ(x2,x1) (41) = x1+(x2−x1)ϕ(x1,x2), B{x1,x2}x2 = x22ϕ(x1,x2)+x21Rx2x1ϕ(x2,x1) (42) = x21+(x22−x21)ϕ(x1,x2),

where is given by Eq. (38). Using Eqs. (40), (41), and (42) we find

 M{x1,x2}z−1≈C1(x1,x2)(z−x2)+C2(x1,x2)ϵ2, (43)

where

 C1(x1,x2) = αx2−(1+α)x1 (44) − (1+α)(x2−x1)ϕ(x1,x2), C2(x1,x2) = α2−{αx2−(1+α)x1}2 (45) + (1+α)(x2−x1){(α−1)x2 − (1+α)x1}ϕ(x1,x2).

Using a counterpart of Eq. (16) for the survival probability of a single -mutant in a stable -dimorphic population we find that

 Px2+ϵ(x1,x2)≈2ϵ(x2−x1)(α−(x2−x1)−1x1−(1+α)ϕ(x1,x2)). (46)

Similarly, consider the possibility that the second branching occurs on the lower branch. In this case we have for

 M{x1,x2}z−1≈C1(x2,x1)(z−x1)+C2(x2,x1)ϵ2 (47)

which yields

 Px1−ϵ(x1,x2) ≈ 2ϵ(x2−x1)(α+(x2−x1)−1x2 (48) −(1+α)ϕ(x2,x1)).

Close inspection of Eqs. (46) and (48) for using Eq. (38) reveals that

 Px1−ϵ(x1,x2)>Px2+ϵ(x1,x2) given |x1|

As explained in Section 5.4, this is one of the factors ensuring the negative dependence among the steps of the random walk in the --plane leading to the condition Eq. (26).

Applying Eq. (26) we can further refine our previous analysis and deduce from Eq. (43)

 M{x1,x2}z−1≈C1(x1,x2)(z−x2)+C2(−x2,x2)ϵ2. (50)

Replacing by in this expression incurs error of order . This implies:

 M{x1,x2}z−1=C1(−x2,x2)(z−x2)+O(ϵ3/2). (51)

It follows from Eq. (39) that

 C1(−x,x) =(1+2α)x−2(1