A mapping of the stochastic Lotka-Volterra model to models of population genetics and game theory

TL;DR

This paper establishes a theoretical connection between the stochastic Lotka-Volterra model and various models in population genetics and game theory, enabling new analytical insights into species extinction dynamics.

Contribution

It provides precise conditions under which the stochastic Lotka-Volterra model maps onto Moran models with different selection regimes, bridging ecological and genetic modeling frameworks.

Findings

Derived conditions for mapping SLVC to Moran models

Calculated extinction probabilities and times

Linked ecological competition to genetic and game-theoretic models

Abstract

The relationship between the M-species stochastic Lotka-Volterra competition (SLVC) model and the M-allele Moran model of population genetics is explored via timescale separation arguments. When selection for species is weak and the population size is large but finite, precise conditions are determined for the stochastic dynamics of the SLVC model to be mappable to the neutral Moran model, the Moran model with frequency-independent selection and the Moran model with frequency-dependent selection (equivalently, a game-theoretic formulation of the Moran model). We demonstrate how these mappings can be used to calculate extinction probabilities and the times until a species' extinction in the SLVC model.

Tables

Table 1. Table 1: Summary of the mappings between the Moran model in various forms and the SLVC model. Mappings are valid at long times in the limit where N 𝑁 N is large but finite, and s 𝑠 s and ϵ italic-ϵ \epsilon are small.

	All cases	Neutral	Constant selection	Frequency-dependent selection
Moran type	Eqs. (14), (17)	Eq. (16), ($\underset{¯}{A}\left(\underset{¯}{x}\right)=0$)	Eqs. (15), (16)	Eqs. (20), (16)
equations
SLVC type	Eqs. (3), (5)	Eq. (4)	Eq. (4)	Eqs. (4)
equations
Mapping	$x_a^{\left(M\right)}=c_0{x}_a^{\left(\mathrm{LV}\right)}/\left(b_0-d_0\right)$	$ϵ=s\equiv 0$	${\rho }_{\alpha }=\left(b_0{\beta }_{\alpha }-d_0{\delta }_{\alpha }\right)/\left(b_0-d_0\right)$	$g_{\alpha M}=\left(b_0{\beta }_{\alpha }-d_0{\delta }_{\alpha }\right)/\left(b_0-d_0\right)$
	$N^{\left(M\right)}=\left(b_0-d_0\right)V/c_0$		$-{\gamma }_{aM}$	$-{\gamma }_{aM}$
	$\equiv ⟨N^{\left(\mathrm{LV}\right)}⟩$		$s=ϵ\left(b_0-d_0\right)/b_0$	$s=ϵ\left(b_0-d_0\right)/b_0$
	$t^{\left(M\right)}=\left[b_0\left(b_0-d_0\right)/c_0\right]t^{\left(\mathrm{LV}\right)}$		${\gamma }_{ab}-{\gamma }_{aM}-{\gamma }_{Mb}+{\gamma }_{MM}=0$	${\gamma }_{ab}-{\gamma }_{aM}-{\gamma }_{Mb}+{\gamma }_{MM}=$
	${\tau }^{\left(M\right)}=b_0{\tau }^{\left(\mathrm{LV}\right)}$			$-\left[g_{ab}-g_{aM}-g_{Mb}+g_{MM}\right]$

Equations

$$T_{\alpha +}(n_{\alpha }+1|n_{\alpha })$$

$$T_{\alpha -}(n_{\alpha }-1|n_{\alpha })$$

$$\frac{\mathrm{d}{P}_{\mathbf{n}}(t)}{\mathrm{d}t}$$

$$-T_{\alpha +}(\mathbf{n}$$

$$\frac{\partial P(\mathbf{x},t)}{\partial t}$$

$$A_{\alpha }(\mathbf{x})$$

$$B_{\alpha \alpha }(\mathbf{x})$$

$$\frac{\mathrm{d}{x}_{\alpha }}{\mathrm{d}\tau }=A_{\alpha }(\mathbf{x})+\frac{1}{\sqrt{V}}{\eta }_{\alpha }(\tau ),\alpha =1,\dots ,M,$$

$$⟨{\eta }_{\alpha }(\tau ){\eta }_{\beta }({\tau }^{\prime })⟩=B_{\alpha \beta }(\mathbf{x})\delta (\tau -{\tau }^{\prime }),$$

$$T(n_1,\dots ,n_a+1,\dots ,n_b-1,\dots ,n_{M-1}|\underline{n})=\frac{n_a}{N}\frac{n_b}{N},$$

$$T(n_1,\dots ,n_a\pm 1,\dots ,n_{M-1}|\underline{n})=\frac{n_a}{N}\frac{N-{\sum }_{b=1}^{M-1}{n}_b}{N},$$

$$T(n_1,\dots ,n_a+1,\dots ,n_b-1,\dots ,n_{M-1}|\underline{n})=\frac{W_a{n}_a}{\mathcal{W}}\frac{n_b}{N},$$

$$T(n_1,\dots ,n_a+1,\dots ,n_{M-1}|\underline{n})=\frac{W_a{n}_a}{\mathcal{W}}\frac{N-{\sum }_{b=1}^{M-1}{n}_b}{N},$$

$$T(n_1,\dots ,n_a-1,\dots ,n_{M-1}|\underline{n})=\frac{n_a}{N}\frac{W_M\left[N-{\sum }_{b=1}^{M-1}{n}_b\right]}{\mathcal{W}},$$

$$\mathcal{W}=\sum _{a=1}^{M-1}{W}_a{n}_a+W_M\left[N-\sum _{b=1}^{M-1}{n}_b\right].$$

$$W_{\alpha }=1+s{\rho }_{\alpha },$$

$$\frac{\partial P}{\partial t}$$

$$A_a(\underline{x})=s{x}_a\left[{\rho }_a-\sum _{b=1}^{M-1}{\rho }_b{x}_b-{\rho }_M\left(1-\sum _{b=1}^{M-1}{x}_b\right)\right],$$

$$B_{ab}(\underline{x})=2\left(x_a{\delta }_{ab}-x_a{x}_b\right)+\mathcal{O}\left(s\right).$$

$$\frac{\mathrm{d}{x}_a}{\mathrm{d}\tau }=A_a(\underline{x})+\frac{1}{\sqrt{N}}{\eta }_a(\tau ),a=1,\dots ,M-1,$$

$$⟨{\eta }_a(\tau ){\eta }_b({\tau }^{\prime })⟩=B_{ab}(\underline{x})\delta (\tau -{\tau }^{\prime }).$$

$$W_{\alpha }(\underline{n})=1+s\left[\sum _{b=1}^{M-1}{g}_{\alpha b}\frac{n_b}{N}+g_{\alpha M}\left(1-\sum _{b=1}^{M-1}\frac{n_b}{N}\right)\right],$$

$$A_a(\underline{x})$$

$${\mathcal{G}}_{a\beta }\equiv g_{a\beta }-g_{M\beta };G_{ab}\equiv {\mathcal{G}}_{ab}-{\mathcal{G}}_{aM}.$$

$$A_{\alpha }(\mathbf{x})$$

$$B_{\alpha \alpha }(\mathbf{x})$$

$$\sum _{\beta =1}^M{x}_{\beta }=(b_0-d_0)c_0^{-1}.$$

$$\frac{\mathrm{d}{y}_{\alpha }}{\mathrm{d}\tilde{\tau }}={\tilde{A}}_{\alpha }(\mathbf{y})\equiv \frac{c_0}{(b_0-d_0{)}^2}{A}_{\alpha }(\mathbf{x})=y_{\alpha }\left[1-\sum _{\beta =1}^M{y}_{\beta }\right],$$

$$J_{\alpha \beta }=\frac{\partial {\tilde{A}}_{\alpha }}{\partial y_{\beta }}={\delta }_{\alpha \beta }\left[1-\sum _{\gamma =1}^M{y}_{\gamma }\right]+y_{\alpha }\left[-1\right],$$

$$\displaystyle\bm{u}^{(M)}=\left(\begin{array}[]{c}1\\ \vdots\\ 1\end{array}\right),\ \ \bm{v}^{(M)}=\left(\begin{array}[]{c}y_{1}\\ \vdots\\ y_{M}\end{array}\right),$$

$$\frac{x_a^{\mathrm{CM}}}{{\sum }_{\beta =1}^M{x}_{\beta }^{\mathrm{CM}}}=\frac{x_a^{(0)}}{{\sum }_{\beta =1}^M{x}_{\beta }^{(0)}}.$$

$$y_a^{\mathrm{CM}}=\frac{y_a^{(0)}}{{\sum }_{\beta =1}^M{y}_{\beta }^{(0)}}.$$

Full text

I Introduction

Perhaps the most important models in ecology, population genetics and game theory are respectively the generalized Lotka-Volterra (LV) model Roughgarden (1979), the Moran model Moran (1957) (and its discrete-generation variant the Wright-Fisher model Fisher (1930)), and the replicator equations Hofbauer and Sigmund (1998). The generalized LV model describes the dynamics of an arbitrary number of species interacting in a pairwise fashion according to an interaction matrix (which can be used to describe competitive, mutualistic and predatory interactions), and is almost always treated deterministically Pielou (1977). The Moran model describes the evolution of a population of individuals carrying different alleles in a way that accounts for genetic drift, and is therefore inherently stochastic Halliburton (2004). The replicator equations describe the time-evolution of the frequency of players playing a given strategy in a pairwise game, where strategies increase according to the average payoff players receive from that strategy when playing against the population Hofbauer and Sigmund (1998). Historically these models were viewed deterministically, however the last decade has seen a surge of interest in incorporating and analyzing stochasticity in these models Nowak et al. (2004); Traulsen et al. (2005, 2008), typically through casting this as a birth-death process analogous to the Moran model.

In this paper, we begin by developing a stochastic analog of the generalized LV model for $M$ species. This stochastic Lotka-Volterra competition (SLVC) model does not have a fixed population size (the number of individuals in the system is free to vary). We then set about determining the conditions under which its stochastic behavior can be seen at long time to be equivalent to the neutral haploid multiallelic Moran model Kimura (1956), the haploid multiallelic Moran model with constant selection Zeng et al. (1989); Mallet (2012) (variously termed the Moran model with directional selection Otto et al. (2008) or frequency-independent selection Mallet (2012)) and a Moran version of a game-theoretic model of pairwise games with multiple strategies Antal et al. (2009). The Moran model in each instance features a population of fixed size $N$ (the number of individuals in the system is fixed). We will show that the SLVC model can be mapped onto each of the above mentioned processes under conditions which we summarize in tabular form. Our analysis relies on timescale separation arguments, and is dependent on the process of population regulation occurring on a much faster timescale than that of the change in population composition. Our results are thus valid when selection between species/alleles/strategies is weak and $N$ is large but finite.

Although the links between the models which we identify have not been discussed previously, other ways in which they are related to each other have been explored, and we will now review these. It is perhaps not surprising that the deterministic versions of these models have been studied far more thoroughly than their stochastic analogs. It is well known that the deterministic LV model in $M$ variables can be mapped to the replicator equations in $M+1$ variables by the inclusion of an additional variable in the LV model that keeps track of population size Hofbauer and Sigmund (1998). While this first result by Hofbauer Hofbauer (1981) allowed for a plethora of results to be obtained, these were all entirely in the deterministic limit (see Bomze (1995) for a comprehensive review).

For our purposes, the most relevant work was not until the publication of Ref. Nowak et al. (2004), when interest was ignited in incorporating demographic stochasticity into the replicator equations. Here it was first demonstrated that the Moran model could be formulated with reproduction rates that vary with population composition to give analogous dynamics to the replicator equations in the infinite population size and weak selection limit. Since then various choices have been explored for how reproduction rates in the Moran model might be dependent on population composition (e.g. linear Nowak et al. (2004), exponential Traulsen et al. (2008)). These all share the common feature that they become functionally similar in the limit of weak selection and by construction also to the replicator equation in the limit of infinite population size. However, this equivalence between game theoretic formulations of the Moran model and the replicator equations has been shown only to hold in the limit of weak selection Traulsen et al. (2005). When selection is strong, alternative update processes, such as an imitation process, map more cleanly to deterministic replicator dynamics Traulsen et al. (2005).

More recently, efforts have been made to understand the role of demographic noise in the LV model by developing stochastic LV models Khasminskii and Klebaner (2001); Spagnolo et al. (2003); Zelnik et al. (2015). Quite naturally with this has come a desire to understand how these probabilistic systems might be related to other canonical models. In Noble et al. (2011) it was shown that a similar frequency-dependent Moran model to that used in the game theory literature shares deterministic fixed points with the LV model. This is perhaps unsurprising, as this formulation of the Moran model by construction maps to the replicator equation in the deterministic limit, which in turn is known to map to the LV model. It is then suggested in Noble et al. (2011) that this equivalence in the deterministic limit might be used to try to understand the stochastic dynamics of the LV model. However this misses the crucial point that two systems that share the same deterministic dynamics need not have the same stochastic dynamics. This was demonstrated in the present context in Ref. Parsons et al. (2010), where it was shown that while increasing the longevity of a type in an SLVC model (while keeping its lifetime fecundity fixed) had no effect on the deterministic dynamics, it could affect the stochastic dynamics, increasing the fixation probability of the type. In that paper the emergent differences between the Moran model with constant selection and a formulation of the SLVC model with homogeneous competition were discussed, primarily for the two-allele case. Finally, the SLVC model in two variables has been mapped to the Moran model in a single variable Constable and McKane (2015).

As we have indicated, in this paper we will be using timescale separation arguments in developing a mapping between the SLVC model and three formulations of the Moran model; the neutral model, the model with constant selection and the replicator model. Unlike traditional deterministic approaches that map $M$ LV equations to $M+1$ replicator equations Hofbauer and Sigmund (1998), we will map $M$ SLVC equations to $M-1$ stochastic replicator equations of the Moran type by elimination of the fast transient associated with a fast approach of the system to carrying capacity. When dealing with both the Moran models and the SLVC model, we will exclusively look at the limit of weak selection. We will therefore not have to be overly concerned about the breakdown of equivalence of the Moran formulation of the replicator equations and the replicator equations in large selection strength regimes Traulsen et al. (2005). In applying our dimensional reduction, we will be careful to deal correctly with the noise terms in the SLVC model. This allows us to determine a full stochastic mapping between the SLVC model and the Moran models, rather than simply inferring the mapping based on a deterministic equivalence Noble et al. (2011). We will not consider noise-induced selection effects of the type identified in Parsons et al. (2010), as these lie outside the scope of this paper. However unlike in Ref. Parsons et al. (2010), where competition rates were taken to be symmetric,we will analyze the effect of varying competition rates. This will allow for the extension of the mapping from the SLVC model to the replicator equations, as well as giving the conditions that the competition matrix must satisfy in order for the Moran model with constant selection to be a valid approximation. Finally, this paper will provide a multivariate extension to the work reported in Ref. Constable and McKane (2015). Not only does this extend the treatment given there, but the multivariate analysis also provides a deeper insight into the mapping between the models.

II Model definitions and the mesoscopic formulation

We define the models at the very basic level of individuals which are born and die and where changes occur from one type to another due to the process of competition. In these individual-level models or individual-based models, the state of the system at a given time is specified by how many individuals of the different types are present at that time. The models are essentially defined by giving functional forms for the rates at which transitions from one state to another occur. We will now describe these for each of the models in turn, starting with the SLVC model.

II.1 The SLVC model

As discussed in the Introduction, the system is a population of $n_{1}$ haploid individuals each of which carries an allele of type $1$, $n_{2}$ haploid individuals each of which carries an allele of type $2$, …, $n_{M}$ haploid individuals each of which carries an allele of type $M$. We denote the state of the system by the vector $\bm{n}=(n_{1},\ldots,n_{M})$. Individuals of type $\alpha$ reproduce at a rate $b_{\alpha}$ and die at a rate $d_{\alpha}$, $\alpha=1,\ldots,M$. The total number of individuals, $\sum^{M}_{\alpha=1}n_{\alpha}$ is not fixed, instead it is regulated by the process of competition, which occurs between individuals of type $\alpha$ and $\beta$ at a rate $c_{\alpha\beta}$.

The transition rates from state $\bm{n}$ to a new state $\bm{n}^{\prime}$ are generalizations of those given for the case of two alleles in Ref. Constable and McKane (2015):

[TABLE]

where $\alpha=1,\ldots,M$ and where only the alleles which change in number have been given as arguments of the transition rates (the original state is to the right and the new state to the left). The parameter $V$ is a measure of the size of the system, such as the volume. It will be shortly be used to make a transition to a mesoscopic description, via the diffusion approximation.

Since the transition rates in Eq. (1) only depend on the current state of the system, the process is Markovian, and can be described by a master equation for the probability, $P_{\bm{n}}(t)$, of finding the system in state $\bm{n}$ at time $t$ van Kampen (2007). It is given by

[TABLE]

where $\bm{\nu}_{\alpha+}$ specifies the number of individuals of type $\alpha$ which increase during the reaction $\alpha\,+$ and $\bm{\nu}_{\alpha-}$ specifies the number of individuals of type $\alpha$ which decrease during the reaction $\alpha\,-$. So, $\bm{\nu}_{\alpha+}=(0,\ldots,0,1,0,\ldots,0)$ and $\bm{\nu}_{\alpha-}=(0,\ldots,0,-1,0,\ldots,0)$, the nonzero entries being in the $\alpha^{\rm th}$ position. Equations (1) and (2), together with an initial condition for $P_{\bm{n}}$, completely specify the stochastic dynamics, so that we can, in principle, find $P_{\bm{n}}(t)$ for all $t$.

We now make the diffusion approximation, mentioned above, that is, $V$ is assumed sufficiently large that $x_{\alpha}\equiv n_{\alpha}/V$ is approximately continuous Crow and Kimura (1970). The other aspect of the approximation involves expanding out the master equation as a power series in $V^{-1}$, and neglecting powers of $V^{-3}$ and higher. The master equation for $P_{\bm{n}}(t)$ then becomes a Fokker-Planck equation (FPE) for $P(\bm{x},t)$ Gardiner (2009):

[TABLE]

The precise form of the functions $A_{\alpha}(\bm{x})$ and $B_{\alpha\beta}(\bm{x})$ are found by carrying out the expansion, but explicit expressions for them also exist in terms of the $\bm{\nu}_{\alpha\pm}$ and the transition rates McKane et al. (2014). One finds that

[TABLE]

and $B_{\alpha\beta}=0$, for $\alpha\neq\beta$.

The FPE (3) is useful for systems with one degree of freedom. However for those with more than one degree of freedom, it is very difficult to analyse, and just as importantly, it is difficult to understand intuitively. For this reason we move over to the completely equivalent, but very different, formulation in terms of a set of stochastic differential equations (SDEs). For the FPE (3) these take the form Gardiner (2009)

[TABLE]

where $\tau=t/V$, $\eta_{\alpha}(\tau)$ is a Gaussian white noise with zero mean and with a correlator

[TABLE]

and where the SDE is to be interpreted in the sense of Itō. Equations (5) and (6) together give the mesoscopic description of the system. The familiar, deterministic, Lotka-Volterra equations form the macroscopic description, and can be found by taking the $V\to\infty$ limit of Eq. (5).

II.2 The Moran model

If we are to discuss the relationship of this model to the Moran model, we need to carry out a similar derivation to that given above, but for the Moran model, since we are not aware that the master equation or the Fokker-Planck equation appears in the literature for the case of $M$ alleles and selection, or at least not in the form that we require here. The derivation itself looks more complicated than the one for the SLVC model, due mainly to the fact that we have to implement the fixed $N$ constraint by expressing one variable in terms of the other $M-1$. Therefore we will only give the definition of the model in terms of the transition rates, and the final form for the FPE and the SDEs here, leaving the intermediate steps to Appendix A.

The states of the system will be labelled by $n_{1},\ldots,n_{M-1}$, since $n_{M}$ can be expressed in terms of the other $(M-1)$ through $n_{M}=N-\sum^{M-1}_{a=1}n_{a}$. We will also use the notation $\underline{n}=(n_{1},\ldots,n_{M-1})$. If we write $n_{M}$, then it should be understood as being equal to $n_{M}=N-\sum^{M-1}_{a=1}n_{a}$. In what follows Greek indices $\alpha,\beta,\gamma,\ldots$ always run from $1$ to $M$ and Roman indices $a,b,c,\ldots$ always run from $1$ to $(M-1)$.

II.2.1 The neutral Moran model

First of all, suppose there is no selection. This case is discussed in the literature; the transition rates are given by Blythe and McKane (2007)

[TABLE]

with $a\neq b$, and

[TABLE]

if either allele $a$ increases at the expense of allele $M$, or allele $M$ increases at the expense of allele $a$, respectively.

II.2.2 The Moran model with frequency-independent selection

We can now add constant selection, that is, selection for each allele that does not depend on population composition. Suppose that $W_{\alpha}$ is the fitness weighting of allele $\alpha$, $\alpha=1,\ldots,M$. Then Eqs. (7) and (8) become

[TABLE]

if $a\neq b$;

[TABLE]

and

[TABLE]

where

[TABLE]

We can simplify these expressions somewhat by taking the limit of weak selection. To do this we express $W_{\alpha}$ as equal to unity plus a small deviation of order $s$:

[TABLE]

where $\rho_{\alpha}$ is of order one but can be positive or negative. However this still leads to rather cumbersome expressions and the details are given in Appendix A. There we also show that going over to new (continuous) variables $x_{a}=n_{a}/N$, one finds that the system can be described by the FPE

[TABLE]

with

[TABLE]

to first order in $s$ and

[TABLE]

As with our treatment of the SLVC model, we note that the above FPE is equivalent to an Itō SDE

[TABLE]

where $\tau=t/N$ and $\eta_{a}(\tau)$ is a Gaussian white noise with zero mean and with a correlator

[TABLE]

This is very similar to Eqs. (5) and (6), but with indices $a$ and $b$ replacing $\alpha$ and $\beta$ and $N$ replacing $V$, and with the functions $A_{a}$ and $B_{ab}$ taken from Eq. (15) and Eq. (16).

The result for $B_{ab}(\underline{x})$ is just that of the neutral case, and has been known for a long time Kimura (1955). The result for $A_{a}(\underline{x})$ can be checked by directly calculating $\mathrm{d}x_{a}/\mathrm{d}\tau$ from the master equation, as described in Appendix A. As we are explicitly considering the mapping between the models in the limit of weak selection, there is no need to go to higher order in $s$.

II.2.3 The Moran model with frequency-dependent selection

We now consider the dynamics of the Moran model if the selective advantage experienced by an allele is dependent on the composition of the population. In this scenario the fitness of an allele $a$ is now denoted by $W_{a}(\underline{n})$, the inclusion of the explicit $\underline{n}$ argument indicating the dependence of the fitness on the nature of the population. The equations for the transition rates then take a similar form to those in the case when selection was constant, Eqs. (9)-(12), but with $W_{a}$ replaced with $W_{a}(\underline{n})$.

We are now left with a choice about how the fitness function $W_{a}(\underline{n})$ depends on the population composition. A common approach is to set $W_{a}(\underline{n})$ to a constant reproductive rate, moderated by a payoff from a game that each allele “plays” with every other allele in the population Nowak (2006). There are many distinct ways to implement this, however, in line with Antal et al. (2009) we make the specific choice

[TABLE]

where $g_{\alpha\beta}$ is the payoff to allele $\alpha$ from interacting with type $\beta$.

As in the case of the Moran model with constant selection (addressed in Section II.2.2), we can expand the master equation in terms of $1/N$ and $s$, and assuming that $N$ is large and $s$ small (formally $s\approx N^{-1}$), obtaining an approximation for the system dynamics in terms of an FPE of form Eq. (14). This is discussed in Appendix A, where it is shown that in this case $\underline{A}(\underline{x})$ is given by

[TABLE]

to first order in $s$, while the form of $B(\underline{x})$ remains unchanged from that given in Eq. (16). Here we have defined the quantities

[TABLE]

Once again, this FPE is equivalent to an SDE of the form Eq. (17), but with $\underline{A}(\underline{x})$ taken from Eq. (20).

It is interesting that it is the quantities $\mathcal{G}_{aM}$ and $G_{ab}$ which appear in the final expression for $A_{a}(\underline{x})$, and not simply $g_{ab}$. The quantity $\mathcal{G}_{a\beta}$ can be interpreted as a relative fitness, namely the payoff to allele $a$ against an opponent $\beta$ relative to the payoff to allele $M$ against the same opponent. Similarly, $G_{ab}$ is a relative relative fitness, namely the relative payoff to allele $a$ against an opponent with an allele $b$ relative to the relative payoff against an opponent with an allele $M$.

At this order in $s$ the dynamics of the system in the deterministic $N\rightarrow\infty$ limit are equivalent to replicator dynamics Nowak et al. (2004); Hofbauer and Sigmund (1998). As we discuss in the Introduction this equivalence does not hold at higher orders in $s$ Traulsen et al. (2005). However, since we will work for the remainder of the paper in the weak selection limit, the mapping that we will develop between the SLVC model and the Moran model with frequency-dependent selection can also be interpreted as a mapping between the SLVC model and a stochastic version of the replicator dynamics.

III Reduction of the Lotka-Volterra Model

In this section we will show that at medium to long-times the LV model with $M$ degrees of freedom, reduces to an $(M-1)$-dimensional model. We can then ask if there are any similarities between this reduced SLVC model and the Moran model. The reduction is accomplished by the systematic elimination of a fast mode, using techniques that we developed previously Constable and McKane (2014a, b, 2015). These require that we first understand the broad features of the deterministic ($V\to\infty$) dynamics, before we go on to study the stochastic dynamics. We begin with the neutral ($s=0$) model; effects due to selection will be introduced as perturbative corrections to the neutral case, since, as usual, we expect $s$ to be small.

III.1 The neutral model

The assumption that individuals of type $\alpha$, $\alpha=1,\ldots,M$, have equal fitness, that is, the theory is neutral, implies that they all have equal birth, death and competition rates: $b_{\alpha}\equiv b_{0},d_{\alpha}\equiv d_{0},c_{\alpha\beta}\equiv c_{0}$. In this case Eq. (4) reduces to

[TABLE]

To characterize the deterministic dynamics, we first determine the fixed points of the dynamics by setting $A_{\alpha}(\bm{x})=0$ for all $\alpha$. It is clear that there are two classes of fixed points depending on whether $\sum^{m}_{\beta=1}x_{\beta}$ is or is not equal to $(b_{0}-d_{0})c^{-1}_{0}$. If it is not, then $x_{\alpha}=0$ for all $\alpha$. So there is a fixed point at the origin and an $(M-1)$-dimensional hyperplane of fixed points given by

[TABLE]

It is useful at this stage to rescale the $x_{\alpha}$ variables and time in order to eliminate the constants $b_{0},c_{0}$ and $d_{0}$ from as much of the calculation as possible. To do so we introduce the new variables $y_{\alpha}=c_{0}x_{\alpha}/(b_{0}-d_{0})$, $\alpha=1,\ldots,M$ and a new timescale $\tilde{\tau}=(b_{0}-d_{0})\tau$. Then the deterministic dynamics becomes

[TABLE]

using Eq. (22). The fixed points are now the origin and the $(M-1)$-dimensional hyperplane $\sum^{M}_{\beta=1}y_{\beta}=1$.

Further insight can be gained by calculating the Jacobian

[TABLE]

at points on this hyperplane. This is a highly degenerate matrix, with all columns identical to each other and equal to the matrix with entries $-y_{\alpha}$. It follows that there is a degenerate set of $(M-1)$ eigenvalues equal to zero, reflecting the existence of the hyperplane of fixed points. The remaining eigenvalue is equal to the trace of the Jacobian: $-\sum^{M}_{\beta=1}y_{\beta}=-1$. We will label the zero eigenvalues $\lambda^{(a)}=0$, $a=1,\ldots,M-1$ and the non-zero eigenvalue $\lambda^{(M)}$. In addition, we will denote the left- and right-eigenvectors of the Jacobian corresponding to the eigenvalue $\lambda^{(\alpha)}$ by $\bm{u}^{(\alpha)}$ and $\bm{v}^{(\alpha)}$ respectively. They will be normalised so that $\sum^{M}_{\gamma=1}\,u^{(\alpha)}_{\gamma}v^{(\beta)}_{\gamma}=\delta_{\alpha\beta}$. For example, the left- and right-eigenvectors corresponding to the non-zero eigenvalue $\lambda^{(M)}=-1$ are

[TABLE]

respectively.

It is now possible to describe the neutral deterministic dynamics rather simply. If the system starts away from the hyperplane $\sum_{\alpha}y_{\alpha}=1$, it will move towards the hyperplane at a rate governed by the non-zero eigenvalue of the Jacobian, that is, unity in the rescaled time and $(b_{0}-d_{0})^{-1}$ in the original time variable. Once it reaches the hyperplane, it remains at this point, since all points of the hyperplane are fixed points. The hyperplane is thus a center-manifold (CM) of the dynamics Wiggins (2003). This dynamics is, of course, so simple that it is of limited interest, except that it forms the basis of the method that we will use to investigate the stochastic dynamics and of the dynamics with selection. For example, under the influence of weak noise, we would expect the system to similarly collapse onto the CM, but with a noisy trajectory, and once on the CM to move around purely stochastically. For this reason the direction perpendicular to the CM, $\bm{v}^{(M)}$, is called the fast direction and the other directions, $\bm{v}^{(a)}$, $a=1,\ldots,(M-1)$, the slow directions.

In our subsequent analysis we will need to map the initial condition of the SLVC model, $\bm{y}^{(0)}$, to an initial condition on the CM $\underline{y}^{\rm CM}$. This is most easily achieved by first introducing new variables to describe the deterministic dynamics of the neutral model. These are $f_{a}=x_{a}/\sum^{M}_{\beta=1}x_{\beta}$, $a=1,\ldots,M-1$, and $\psi=\sum^{M}_{\beta=1}x_{\beta}$. By direct substitution into the dynamical equations $dx_{\alpha}/dt=A_{\alpha}(\bm{x})$, with $A_{\alpha}(\bm{x})$ given by Eq. (22), one finds that $df_{a}/dt=0$, that is the $f_{a}$ are constants of the motion. This implies that they do not change from their initial values: $f_{a}=f^{(0)}_{a}$ at all times $t$. The only dyamical variable is $\psi$, which satisfies the equation $d\psi/dt=\psi[(b_{0}-d_{0})-c_{0}\psi]$. Therefore the dynamics simply consists of $\psi(t)$ decreasing (increasing) if it is larger (smaller) than $(b_{0}-d_{0})/c_{0}$, until it reaches this value. At this point the deterministic dynamics ceases, as the system is now on the CM. Using $f_{a}=f^{(0)}_{a}$ this point, labelled with the superscript CM, is given by

[TABLE]

Using $\sum^{M}_{\beta=1}x^{\rm CM}_{\beta}=(b_{0}-d_{0})/c_{0}$, and going over to the scaled variables $y_{\alpha}$, one finds that the initial condition on the CM is given by

[TABLE]

A useful mathematical tool to separate the fast and slow dynamics is the projection operator defined by

[TABLE]

Suppose it is used to operate on the vector $\phi_{\beta}=\sum^{M}_{\kappa=1}C_{\kappa}v^{(\kappa)}_{\beta}$, where the $C_{\kappa}$ are arbitrary constants. Then

[TABLE]

that is, the fast term $C_{M}v^{(M)}_{\beta}$ has been wiped out.

We can use the projection operator directly as given by Eq. (35), or observe that

[TABLE]

by completeness. Since $u^{(M)}_{\beta}=1$ for all $\beta$, we may write this as

[TABLE]

We also note that $\sum^{M}_{\alpha=1}P_{\alpha\beta}=0$.

The projection operator can be used to determine the stochastic dynamics on the CM. In terms of the rescaled variables Eq. (5) becomes

[TABLE]

since from Eq. (22) we see that on the CM $A_{\alpha}(\bm{x})=0$. We also note that $B_{\alpha\beta}=2b_{0}x_{\alpha}\delta_{\alpha\beta}$ on the CM, and so

[TABLE]

Application of the projection operator in Eq. (38) to $\mathrm{d}y_{\beta}/\mathrm{d}\tilde{\tau}$, which we denote as $\dot{y}_{\beta}$, gives $\dot{y}_{\alpha}-\sum^{M}_{\beta=1}v^{(M)}_{\alpha}\dot{y}_{\beta}$. However on the CM, $\sum^{M}_{\beta=1}\dot{y}_{\beta}=0$, and we recover $\dot{y}_{\alpha}$. Therefore, defining a projected noise $\zeta_{\alpha}=\sum^{M}_{\beta=1}P_{\alpha\beta}\tilde{\eta}_{\beta}$, the SDE (39) becomes

[TABLE]

Although, Eqs. (39) and (41) look similar, there are some significant differences. First, in Eq. (41), there are only $(M-1)$ independent variables (since $\sum_{\alpha}y_{\alpha}=1$) and noises (since $\sum_{\alpha}P_{\alpha\beta}=0$ implies that $\sum_{\alpha}\zeta_{\alpha}=0$). Second, the noises, $\zeta_{\alpha}(\tilde{\tau})$, now have a different correlation function:

[TABLE]

using Eqs. (38) and (40).

In summary, the noisy dynamics on the CM is governed by the SDEs in Eq. (41), where the noise correlation is given by Eq. (LABEL:correlator_new). However, as pointed out above, only $(M-1)$ of the $y_{\alpha}$ and the noises are independent. We therefore choose the dynamical variables to be the first $(M-1)$ $y_{\alpha}$ and denote these by $z_{a}$.

The SDEs then become

[TABLE]

where $\zeta_{a}(\tilde{\tau})$ is a Gaussian noise with zero mean and with a correlator

[TABLE]

The choice of independent variables is chosen to mirror those in the Moran model: the first $(M-1)$ being independent, with the final one being determined through the condition $y_{M}=1-\sum^{M-1}_{a=1}y_{a}$. In the Moran model this condition comes from the constraint $\sum^{M}_{\alpha=1}n_{\alpha}=N$, whereas in the SLVC model it comes from the equation of the CM.

III.2 The model with selection

Introducing selection into the neutral model implies that the birth, death and competition rates can now be different for different alleles. We write

[TABLE]

where $\epsilon$ is a small selection constant, which will be taken to be proportional to the selection constant, $s$, in the Moran model, when we compare both models.

The SLVC model to first order in $\epsilon$ can now be constructed. There is no need to modify the $B_{\alpha\beta}(\bm{x})$ beyond their neutral form, as working in the weak selection limit, $\epsilon$ is small, and terms of order $\epsilon/N^{2}$ in the FPE can be neglected. The functions $A_{\alpha}(\bm{x})$ are given by (after rescaling as in Eq.(24))

[TABLE]

If we now search for fixed points of the dynamics: values of $\bm{y}$ such that $\tilde{A}_{\alpha}(\bm{y})=0$, it is found that generically there is at most one fixed point which is not on the boundary of the allowed region of $\bm{y}$ variables. This result is known to be true for the LV model Hofbauer and Sigmund (1998); Noble et al. (2011); Gokhale and Traulsen (2010); we infer the result for the reduced model by exploiting the mapping to the Moran model with frequency-dependent selection, which we will demonstrate in Sec. IV.3. We therefore deduce that a CM does not exist for nonzero $\epsilon$.

Although a CM no longer exists, we would still expect the time-scale separation argument used previously for the neutral case to apply. The time-scale for the fast mode to collapse will only be changed by small terms of order $\epsilon$, but now the collapse will not be onto a CM, but onto an $(M-1)$-dimensional subspace on which there is a weak deterministic dynamics (of strength $\epsilon$) in addition to the same noise as was found in the neutral case. This subspace will not be planar, but we will nevertheless determine it by asking that there is no deterministic dynamics in the fast direction as defined by the eigenvector $\bm{v}^{(M)}$ found in the neutral case, i.e., we ask that $\bm{u}^{(M)}\cdot\tilde{\bm{A}}(\bm{y})=0$. An equivalent condition, using the explicit form for $\bm{u}^{(M)}$ given in Eq. (32) is $\sum^{M}_{\alpha=1}\tilde{A}_{\alpha}=0$. This is an approximation, but in previous work where it has been used Constable and McKane (2014a, b, 2015) it has been found to be a very good one. The condition is also consistent with the other terms which appear in the SDE, given in the neutral case by Eq. (41), that is, $\sum_{\alpha}\dot{y}_{\alpha}=0$ and $\sum_{\alpha}\zeta_{\alpha}=0$. In a similar fashion we will also continue to use Eq. (34) as an approximation for the initial condition of the system on the slow subspace.

This condition on $\tilde{\bm{A}}(\bm{y})$ determines the equation of the slow manifold. To zeroth order in $\epsilon$ it gives $y_{M}=1-\sum^{M-1}_{a=1}y_{a}$, the equation of the CM. Therefore to determine the order $\epsilon$ correction to this equation which gives the slow manifold we write

[TABLE]

where $f$ is a function to be determined. Substituting Eq. (47) into the condition $\sum^{M}_{\alpha=1}\tilde{A}_{\alpha}=0$ gives

[TABLE]

where we have introduced the combinations of constants

[TABLE]

This is the same combination of $\gamma_{\alpha\beta}$ as appears in the definition of $G_{ab}$ given by Eq. (21), and the same interpretation in terms of relative quantities holds. In fact, using Eq. (45) we may introduce a similar quantity for the full competition rates $c_{\alpha\beta}$:

[TABLE]

We can now eliminate $y_{M}$ from the function $\tilde{A}_{\alpha}(\bm{y})$ given by Eq. (46) to find

[TABLE]

where SS indicates that this is $\tilde{A}_{\alpha}(\bm{y})$ evaluated on the slow-subspace and where, for clarity, we have introduced

[TABLE]

Introducing the reduced variables $z_{a}$ as in the neutral case, we may add the term in Eq. (51) to the SDE given in Eq. (43) to give the reduced SDE in the case with selection:

[TABLE]

where the noise is as in the neutral case, that is, with the correlator given by Eq. (44) and where $\left.\tilde{A}_{a}(\underline{z})\right|_{\textrm{SS}}$ is evaluated at first order in $\epsilon$.

In this section we have shown that even though the $M$-allele SLVC model begins with one more degree of freedom than the $M$-allele Moran model, after some time the extra degree of freedom decays away and the models begin to resemble each other. In the next section we seek to determine the precise conditions under which their dynamics are equivalent.

IV Comparison of the reduced SLVC model and Moran models

IV.1 The neutral models

It is clear that the reduced SLVC model written in the form of Eq. (43) is precisely the neutral Moran model (Eqs. (15) and (16) with $s=0$) up to a constant in the correlation function. The identification can be taken a little bit further by noting that the average population size of the SLVC model on the CM, Eq. (23), in terms of the number of individuals in the system is

[TABLE]

We can then transform Eq. (43) back into a FPE equation in the natural units $t$ of the SLVC process, to obtain

[TABLE]

where

[TABLE]

Comparing Eq. (14) (with $\underline{A}(\underline{x})=0$ and $B(\underline{x})$ taken from Eq. (16)) with Eq. (55), we see that in units of $t$, the natural timescale of the underlying stochastic processes, the timescale of the neutral SLVC model, is related to that of the Moran model by $t^{\mathrm{(M)}}=[b_{0}(b_{0}-d_{0})/c_{0}]t^{\mathrm{(LV)}}$ for fixed $N=\langle N\rangle$.

Some further algebra also allows us to calculate how the neutral Moran model and the neutral SLVC model are related in the SDE representation given Eqs. (17) and (5). Along with the scaling $x_{a}^{(\mathrm{M})}=c_{0}x_{a}^{\mathrm{(LV)}}/(b_{0}-d_{0})$, we find $\tau^{\mathrm{(M)}}=b_{0}\tau^{\mathrm{(LV)}}$. These results are summarized in Table 1.

IV.2 The models with frequency-independent selection

In Sec. IV.1 we showed that the neutral reduced SLVC model and neutral Moran model are identical up to a rescaling of time. In this section we ask under what conditions the reduced SLVC model behaves identically to the Moran model with frequency-independent selection, defined by Eqs. (14)-(16).

This is carried out by comparing the constant terms, the terms linear in $y_{b}$ and the quadratic terms, in the bracket multiplying $\epsilon y_{a}$ in Eq. (51) with the corresponding bracket in Eq. (15). Comparing the constant terms gives $\rho_{a}-\rho_{M}=(\Phi_{a}-\gamma_{aM})-(\Phi_{M}-\gamma_{MM})$. The linear terms now match if $\sum_{b}\Gamma_{ab}y_{b}=0$ for all $y_{b}$, that is, if $\Gamma_{ab}=0$ for all $a$ and $b$. This last condition implies that there is no quadratic term, as is clear from Eq. (15). Thus the condition for the mapping to exist is $\Gamma_{ab}=0$ and the relationship between the $\rho_{a}$ and the $\gamma_{aM}$ that is required in order to carry out the mapping is $\rho_{a}=\Phi_{a}-\gamma_{aM}+$ constant. If we also identify $\rho_{M}=\Phi_{M}-\gamma_{MM}$, then this constant is zero, and we may then state that the reduced SLVC model behaves identically to the Moran model with frequency-independent selection if $\Gamma_{ab}\equiv 0$ and if we make the identification

[TABLE]

Assuming that this mapping can be carried out, we now transform Eq. (51) back into a FPE equation in the natural units $t$ of the SLVC process to obtain

[TABLE]

where $\mathcal{A}_{a}(\underline{z})$ is a rescaled form of $\tilde{A}_{a}(\underline{y})|_{\textrm{SS}}$ — reversing the rescaling carried out in Eq. (46):

[TABLE]

and where $\mathcal{B}(\underline{z})$ retains the form given in Eq. (56). Comparing Eqs. (14)-(16) with Eqs. (58), (59) and (56) we see that the two equations are identical if we rescale time in the SLVC such that $t^{\mathrm{(M)}}=[b_{0}(b_{0}-d_{0})/c_{0}]t^{\mathrm{(LV)}}$ and simultaneously make the identification

[TABLE]

along with the condition $\Gamma_{ab}=0$ for all $a$ and $b$. The equivalence of the reduced SLVC model and Moran model in the SDE setting is obtained by again enforcing the condition on $\Gamma_{ab}$ and (60), but rescaling time instead by $\tau^{\mathrm{(M)}}=b_{0}\tau^{\mathrm{(LV)}}$. We note that is different to the rescaling adopted in the analysis of the $M=2$ version of this correspondence Constable and McKane (2015), where factors of $\gamma_{\alpha\beta}$ were included.

One can also ask about the nature of the fixed points in the Moran model and the reduced SLVC model under the conditions outlined above. From Eq. (15) the fixed points of the model are given by solutions to

[TABLE]

where $\hat{\rho}_{a}\equiv\rho_{a}-\rho_{M}$. An analysis of this equation shows that the only fixed points are on a boundary, unless all the $\rho_{\alpha}$ are equal. However, if all the $\rho_{\alpha}$ are equal there is no selection, so in this case of frequency independent selection there are no interior fixed points.

IV.3 The models with frequency-dependent selection

We can now repeat a similar analysis to that described in Sec. IV.2 to ask under what conditions the reduced SLVC model behaves identically to the Moran model with frequency-dependent selection, defined by Eqs. (14), (20) and (16).

Once again, we compare the constant terms, the terms linear in $y_{b}$ and the quadratic terms, in the bracket multiplying $\epsilon y_{a}$ in Eq. (51) with the corresponding bracket in Eq. (20). Comparing the constant terms gives $g_{aM}-g_{MM}=(\Phi_{a}-\gamma_{aM})-(\Phi_{M}-\gamma_{MM})$. The linear terms now match if $\sum_{b}G_{ab}y_{b}=-\sum_{b}\Gamma_{ab}y_{b}$ for all $y_{b}$, that is, if $G_{ab}=-\Gamma_{ab}$ for all $a$ and $b$. This last condition implies that the quadratic terms match. Thus the parameters of the two models are related by $\Gamma_{ab}=-G_{ab}$ and also by $g_{aM}=\Phi_{a}-\gamma_{aM}+$ constant. If we also identify $g_{MM}=\Phi_{M}-\gamma_{MM}$, then this constant is zero, and we may then state that the reduced SLVC model behaves identically to the Moran model with frequency-independent selection if $G_{ab}=-\Gamma_{ab}$ and if we make the identification

[TABLE]

We may once again carry out rescalings and transformations as in Sec. IV.2 which lead to Eqs. (58), (59) and (60).

Under the conditions outlined above it can be shown that the deterministic frequency-dependent Moran model (equivalently the replicator equations) admits at most one stable fixed point on the interior region Gokhale and Traulsen (2010). This is clearly also true for the reduced version of the SLVC model, since we have shown that it maps to the frequency-dependent Moran model. A more detailed discussion of these deterministic considerations can be found in Noble et al. (2011).

V Utilization of the Mappings

A useful feature of the mappings that have been derived above, is that we can now use analytic results obtained for the Moran model to make predictions about the SLVC model. In this section we illustrate this by mapping various results for fixation probabilities and times in the Moran model to those in the SLVC model, providing predictions for the extinction probabilities and mean time to extinction of species in certain scenarios.

V.1 The neutral models

In the neutral Moran system, it it possible to separate out the dynamics of the different alleles Baxter et al. (2007); rather than consider the dynamics of the many different interacting alleles, we can instead split the population into two subpopulations, one containing an allele of interest and the other containing all of the remaining alleles. Using this approach it is possible to calculate the probability of any series of extinctions and the time until the $r^{\rm th}$ extinction. These results, together with the mappings provided by Table 1, allow us to calculate the same quantities for the SLVC model. The full calculation is detailed in Appendix B, however here we state the main results for the SLVC model.

Let $\underline{S}=(\alpha_{M},\alpha_{M-1},\ldots,\alpha_{2})$ be a series of species extinctions in the neutral $M$-species SLVC model, such that species $\alpha_{M}$ goes extinct first, followed by $\alpha_{M-1}$ finally leaving only species $\alpha_{1}$ in the population. The probability of this series of extinctions, $P_{\mathrm{Fix}}^{\underline{S}}$, given initial conditions $\bm{x}^{(0)}$ is given by

[TABLE]

We can also use this result to calculate the probability that species $\alpha$ goes extinct first, $P_{\mathrm{Ext}}^{\left\{\alpha\right\}}$; this is simply the sum of Eq. (63) over all $\underline{S}$ that do not contain $\alpha_{\alpha}$ as an element (i.e. species $\alpha$ never fixates but rather is the first to reach extinction). We find excellent agreement between this result and those obtained from simulations (see Figure 2).

We now move to considering the mean unconditional time until the extinction of the $r^{\rm th}$ species, $\langle T_{(r)}^{\mathrm{(LV)}}(\bm{x}^{(0)})\rangle$. Given an initial distribution of species $\bm{x}^{(0)}$ in the SLVC model, this is given by

[TABLE]

where $\langle N^{(\mathrm{LV})}\rangle$ is the average number of individuals in the SLVC model at carrying capacity (see Eq. (54)). Here the summation over $\bm{\alpha}=\{\alpha_{1},\alpha_{2},\ldots,\alpha_{a}\}$ denotes summation over all possible subsets of the set of positive integers $\{1,2,\ldots,M\}$ with $a$ entries. Again, this matches very well the results obtained from simulating the underlying stochastic models (see Figure 2).

Note that in the Moran model, extinctions will occur on a time-scale proportional to $N^{2}$ when measured in units of $t^{\mathrm{(M)}}$ (see Appendix B), while in the SLVC model we predict fixation to occur on a timescale $c_{0}\langle N\rangle^{2}/\left[b_{0}(b_{0}-d_{0})\right]$ (see Table 1). From this we see that increasing the birthrate, $b_{0}$, in the SLVC model increases the rate of ‘genetic drift’ in the SLVC model relative to the Moran model. However a perhaps less intuitive result is that increasing the average genotype lifetime $b_{0}-d_{0}$ also increases the relative rate of genetic drift, while conversely increasing competition rate, $c_{0}$, slows down the rate of genetic drift in the SLVC model relative to the Moran model.

V.2 The models with frequency-independent selection

It is very difficult to obtain results on the fixation probability of an allele in a multi-allele Moran model when each of the alleles is under a different selection pressure, and to our knowledge no analytic results for this problem have been obtained. Progress can be made however if one assumes that all of the alleles are under one or other of two selective pressures. In this case the entire system dynamics can be decomposed into two processes; a Moran process with selection between the two subpopulations, and neutral drift within each population. The fixation probability of any allele can then be expressed as the product of the fixation probability of its subpopulation, multiplied by its fixation probability within the subpopulation.

The full calculation is detailed in Appendix B, while here we simply give the key results in the context of the SLVC model. Suppose that the $M$ species in the SLVC model are interacting in such a way that their dynamics can be described by frequency-independent selection (see Eq. (57) and Table 1). Now further suppose that although their birth, death and interaction parameters may all be distinct, they are such that they only give rise to two distinct selection pressures when mapped to the Moran model with frequency-dependent selection; that is

[TABLE]

The fixation probability of any species in the SLVC model is then described by

[TABLE]

where

[TABLE]

and where we have made use of Eqs. (34) and (60). We find excellent agreement between these results and results obtained from simulations, as illustrated in Figure 3. Determining the probability of first extinction is not possible however; this requires knowing the probability of time ordering of extinctions (whether fixation in the subpopulation occurs before fixation in the population as a whole).

More generally, we find good agreement between the Moran model and the SLVC model with conditions and mappings taken from Table 1, even when analytic results are not available (see Figure 5).

V.3 The models with frequency-dependent selection

The inclusion of cubic terms in the multiallelic model model with frequency-dependent selection makes obtaining analytic results for the fixation probability and time even more challenging than in the case with frequency-independent selection. However analytic progress can again be made under the condition that the elements of the payoff matrix can be partitioned such that only two distinct strategies exists within the population (see Appendix B). Here we describe the results of this calculation in relation to the SLVC model.

We first make the assumption that the competition matrix $\gamma_{\alpha\beta}$ can be partitioned such that

[TABLE]

while the birth and death terms can be partitioned such that

[TABLE]

Then the fixation probability of any species in the population is shown in Appendix B to be given by

[TABLE]

where

[TABLE]

Here erfi and erfc are respectively the imaginary and complimentary error functions Abramowitz and Stegun (1965); Erdélyi (1953), and the function $l(\bm{x}^{(0)})$ is defined as

[TABLE]

In addition, the notation $l(1)$ and $l(0)$ means the value of $l$ when $\left(\sum_{b=1}^{\theta}x^{(0)}_{b}\right)/\left(\sum_{\beta=1}^{M}x^{(0)}_{\beta}\right)=1$ and $\left(\sum_{b=1}^{\theta}x^{(0)}_{b}\right)/\left(\sum_{\beta=1}^{M}x^{(0)}_{\beta}\right)=0$, respectively. Using this analytic result we find good agreement between our theory and results obtained from simulations, as illustrated in Figure 4.

We note however that in Figure 4 there begins to be disagreement for the time to fixation. This difference is likely caused by the existence of an interior stable fixed point in the system. In general, our method of fast-variable elimination is expected to perform well when $s$ is small and $N$ is relatively large Constable and McKane (2014a). However, when an interior fixed point is present, another factor comes into play; the validity of the diffusion approximation itself. As the distribution of stochastic trajectories about the fixed point becomes increasingly stable (for instance, by increasing $N$), the FPE obtained via the diffusion approximation is known to become less accurate near the boundaries Doering et al. (2006). Therefore, although our approximation of the FPE becomes more accurate with increasing $N$, the FPE itself becomes less reliable in this regime for increasing $N$. The disagreement hinted at here between the Moran-replicator and SLVC models will then only be calculable outside the diffusion limit.

Solving for fixation times in this scenario becomes increasingly difficult; as we have discussed they are not straightforward to obtain even in a two-type system, as the emergence of very stable fixed points creates deviations between the underlying stochastic process and the diffusion approximation that have to be carefully corrected for Doering et al. (2006). Adding more types compounds this difficulty. However, comparing results of simulations of the Moran model with the SLVC model, we can see that the mapping between the models remains qualitatively intact (see Figure. 4) and that the mapping also holds under a broad range of parameter regimes in which there is no interior fixed point (see Figure 5).

VI Conclusions

In this paper we have defined a mapping between various incarnations of the Moran model in $M-1$ variables and the SLVC model in $M$ variables, which accounts for demographic noise. In particular we have concentrated on the mapping between three forms of the Moran model: the neutral Moran model, the Moran model with frequency-independent selection and the Moran model with frequency-dependent selection. While the Moran model is formulated in terms of a population of fixed size $N$, the populations in the SLVC model do not strictly have a fixed size. However, at long times they approach a carrying capacity around which they stochastically fluctuate. The key to making an analytic bridge between these models has been in noting the following: if selection is weak and the population size is large, the SLVC model approaches its carrying capacity (and is confined in its vicinity) on a much faster timescale than that on which the composition of the population changes. By working in an SDE (or equivalent FPE) setting that approximates the stochastic dynamics, the fast timescales can be clearly identified. On removing these fast transient dynamics, the SLVC model can be approximated by a reduced model in $M-1$ variables that is of similar form to the Moran model, and thus a precise mapping can be determined.

Our analysis begins with a consideration of the neutral SLVC model, which we define as that in which all birth, death and competition rates are equal for all species. This setting allows us to determine a CM for the system and in turn calculate the reduced form of the neutral SLVC model. The neutral setting also serves as a useful reference case for determining the reduced form of the non-neutral SLVC model; the reduced non-neutral dynamics can be calculated using a perturbation theory around the neutral system. We find that the SLVC model maps to the Moran model with frequency-dependent selection with no conditions other than those stated in our approximation, that selection is weak and the population size is large (see Table 1). We note that for any given set of parameters in the SLVC model, a unique choice of payoff matrix in the Moran model does not exist. This is because it is the relative (rather than absolute) values of the payoff matrix that are consequential for the dynamics (see Eq. (21)). The SLVC model maps to the Moran model with frequency-independent selection with conditions on the competition matrix that ensures that there is no frequency-dependent selection in the SLVC model. Finally, and as already stated, the SLVC model maps to the neutral Moran model with the largest number of conditions; that birth, death and competition rates are the same for each species.

Although we have used the neutral SLVC model as a reference case for our fast-variable elimination procedure, other choices of reference case are possible. In fact a similar approach to that outlined in this paper can be taken by implementing a perturbation theory around any choice of parameters that generates a CM in the LV dynamics. More generally, CMs in the LV system can also be achieved by setting $c_{\alpha\kappa}=c_{\beta\kappa}$ and $b_{\alpha}-d_{\alpha}=b_{\beta}-d_{\beta}$ for all $\alpha$ and $\beta$. Although the species are no longer identical (even in the absence of selection), a reduced form of the dynamics can still be calculated. However, breaking the symmetry between the species in this way can give rise to noise-induced selection (see, for instance, Constable et al. (2016), where the competition terms are varied and Parsons et al. (2010); Lin et al. (2015), where the birth and death rates are varied in two-species systems). In this paper we have not explored these possibilities, instead focusing on developing a clean mapping between models in which neutrality is defined in the usual way. However considering the emergence of such noise-induced effects for the $M$-species SLVC model will be of great interest for further work.

In general we have stated that we expect the mapping that we have developed between the SDEs for each model to hold under the same range of validity as the approximations we have employed. From the diffusion approximation that gave rise to the SDEs we require that the number of individuals in both the SLVC model and Moran models is large (large $V$ and $N$ respectively) while from the fast-variable elimination procedure we additionally require that selection is weak (small $\epsilon$ and $s$ respectively). However, as noted in Section V.3, this picture becomes more complicated in the case of frequency-dependent selection when a stable coexistence fixed point is present. In this case, as the system size increases, the diffusion approximation can become an increasingly poor predictor of the fixation properties of the system Doering et al. (2006). Thus although the mapping between the SDEs becomes increasingly accurate with increasing $N$ and $V$, the underlying stochastic dynamics of the systems may not converge in this limit. It is likely that this is the origin of the discrepancy between the fixation times in the SLVC model and Moran model in Figure 4. Further work is required using other approximations of the underlying master-equation, such as the WKB approximation Assaf and Meerson (2010), to determine to what extent the mapping established here remains valid with respect to fixation times. However it should be stressed that although the fixation time becomes poorly matched in this case, the mapping in terms of fixation probabilities still works very well (see Figure 4).

The mappings defined in this paper have a clear utility. Namely, the problem of calculating many stochastic quantities relating to the SLVC model (such as fixation probabilities and times) is reduced to finding a related result in the relevant Moran model together with an implementation of the mapping that we have defined in Table 1. While stochasticity in the Lotka-Volterra model has received relatively little attention, there is a vast literature of results pertaining to the Moran model and its variants. In Section V of this paper, we have illustrated how fixation properties for the SLVC model can be derived in a straightforward way by employing the mapping to the Moran process. No doubt there are more rich and interesting behaviors that can be uncovered in a similar way. Since the SLVC model maps with the least stringent conditions to the Moran model with frequency-dependent selection, the increasing number of analytic results related to multi-strategy game theory Gokhale and Traulsen (2010); Antal et al. (2009) are of special note. While consideration of these numerous possible extensions is beyond the scope of this current paper, the inclusion of mutation in the SLVC model may be of particular interest. This is in part because of the large body of analytic work that exists in the population genetic Vogl and Clemente (2012) and game theory literature Antal et al. (2009) that relies on small mutation rates to gain analytic traction.

The form of the Moran model familiar to most population geneticists is that with frequency-independent selection. By construction this ignores interactions between alleles that alter each other’s fitness. In terms of deriving analytic results on population genetics, understanding the frequency-independent case is clearly the first point of order. However, it is in some sense naive to imagine that nature would conform to this scenario. Frequency-dependent selection is in some sense the hallmark of an ecological system and arises continuously in biological systems. Indeed, despite conscious attempts to remove ecology from microbial experiments, frequency-dependent selection is often seen to emerge Maddamsetti et al. (2015). Just as these experimental findings have started to motivate theoretical studies Huang et al. (2012), from the perspective of theoretical ecology there is an increasing awareness that demographic stochasticity can have important consequences Gokhale et al. (2013). It is our hope that the work presented here will prove to be of use to researchers in ecology, population genetics and game theory in seeing concrete parallels, distinctions and applications in each other’s work.

Acknowledgements.

GWAC thanks the Finnish Center for Excellence in Biological Interactions for funding.

Appendix A Derivation of the Fokker-Planck equation for the Moran model

In Sec. II of the main text we defined the $M$-allele Moran model with selection through the transition rates in Eqs. (9)-(11). We considered two types of selection. In the first type, the fitness weightings, $W_{a}$, were independent of the number of individuals carrying a particular allele, that is, independent of $\underline{n}$. For the second type the fitness weightings depended on the population composition in a way which was given by Eq. (19). In both cases, when these transition rates are substituted into the master equation

[TABLE]

they give the stochastic dynamics of the population. In this Appendix we will derive the FPE for the $M$-allele Moran model with selection, by applying the diffusion approximation to the master equation with these particular transition rates.

To make the diffusion approximation we write $x_{\alpha}=n_{\alpha}/N$ and introduce the notation $F_{\alpha,\beta}(\underline{x})$ for the above transition rates for moving from state $x_{\alpha}$ to state $x_{\alpha}+N^{-1}$ and from state $x_{\beta}$ to state $x_{\beta}-N^{-1}$. In terms of the $\underline{x}$ variables, the master equation (71) becomes, after expanding in powers of $N^{-1}$, relabelling and combining,

[TABLE]

This is the FPE for the Moran model:

[TABLE]

where

[TABLE]

This is as far as one can go without using the specific forms for the transition rates, so now we consider the two cases of frequency-independent selection and frequency-dependent selection in turn.

A.1 The Fokker-Planck equation for the Moran model with frequency-independent selection

Before beginning the derivation, we simplify the expressions in Eqs. (9)-(11) by using Eq. (13) and expanding in $s$. Since $\sum^{M}_{\gamma=1}W_{\gamma}n_{\gamma}=N+s\,\sum^{M}_{\gamma=1}\rho_{\gamma}n_{\gamma}$, we have

[TABLE]

Therefore,

[TABLE]

if $a\neq b$,

[TABLE]

and

[TABLE]

In terms of the $F_{\alpha\beta}(\underline{x})$ these are, omitting terms of order $s^{2}$ and higher,

[TABLE]

Using these specific forms, one finds from Eq. (74) that

[TABLE]

After the introduction of the new time scale $\tau=t/N$, the FPE may be written in terms of the set of $(M-1)$ SDEs (17). The function $A_{a}(\underline{x})$ may be slightly rewritten to give Eq. (15). Similarly the diagonal elements of the functions $B_{ab}(\underline{x})$ may be simplified to $B_{aa}(\underline{x})=2x_{a}-2x^{2}_{a}=2x_{a}(1-x_{a})$, giving Eq. (16).

The result for $B_{ab}(\underline{x})$ is known Kimura (1955). One check on the form of $A_{a}(\underline{x})$ is to take $M=2$, so that $a$ and $b$ take only the value $1$. Then, dropping the index on $A$ and on $x$, we have that $A(x)=s\rho_{1}x(1-x)-sx\rho_{2}(1-x)=s(\rho_{1}-\rho_{2})x(1-x)$, which is the known result for two alleles, up to a constant Crow and Kimura (1970). As a second, and more substantial, check we can start from the master equation with transition rates (9)-(11) and work out the equation for $d\langle n_{k}\rangle/d\tau$. One finds

[TABLE]

There are a number of points to be made about this equation. First, we have replaced $n_{a}$ by $Nx_{a}$. Since $n_{a}$ always appears in the combination $n_{a}/N$, there was no need to take $N\to\infty$ to eliminate extra factors of $N$. Second, there should be angle brackets around all the terms on the right-hand side; these have been omitted so that the expression does not look so cluttered. This is permitted, since in the limit $N\to\infty$, the average of products is the product of the averages, that is, $\langle x_{a}x_{b}\rangle=\langle x_{a}\rangle\langle x_{b}\rangle$. In the following we will also omit the angle brackets on the left-hand side, since we are attempting to derive the macroscopic equation, which is an equation for the macroscopic variable $x_{a}$, written without angle brackets. Finally, we also note that the factor $\sum^{M}_{\beta=1}W_{\gamma}x_{\gamma}$ is common throughout in the denominator, and so we may multiply through by it, to find

[TABLE]

We now write $W_{\alpha}=1+s\rho_{\alpha}$, as in Eq. (13). The terms of order $1$ on the right-hand side are seen to cancel, and so the whole of the right-hand side is of order $s$, and $W_{\alpha}$ may be replaced everywhere by $s\rho_{\alpha}$. In addition the factor in brackets on the left-hand side is $\sum^{M}_{\gamma=1}x_{\gamma}+\mathcal{O}(s)$, and so to the order we are working at is simply equal to $1$. Therefore the macroscopic equation reads

[TABLE]

which has the form $dx_{k}/d\tau=A_{k}(\underline{x})$, where $A_{k}(\underline{x})$ is given by Eq. (15).

A.2 The Fokker-Planck equation for the Moran model with frequency-dependent selection

In this case the average fitness of the population (see Eq. (12)) can be expressed as

[TABLE]

This then leads to the following results for the combinations of transition rates which are of interest to us:

[TABLE]

and

[TABLE]

Substituting these expressions into Eq. (74), one finds that

[TABLE]

with the form of $B_{(}\underline{x})$ being the same as in the frequency-independent case, given in Eq. (79). Now interchanging the $b$ and $c$ labels in the sum involving $g_{Mb}x_{b}x_{c}$ in Eq. (86), and introducing $\mathcal{G}_{aM}=g_{aM}-g_{MM}$ and $G_{ab}=g_{ab}-g_{aM}-g_{Mb}+g_{MM}$ one obtains Eq. (20) in the main text.

As for the case of frequency-independent selection, we can check this result for $A_{a}(\underline{x})$ by calculating $d\langle n_{k}\rangle/d\tau$ directly from the master equation. The steps leading to Eq. (81) still hold — as long as $W_{\alpha}$ is replaced by $\langle W_{\alpha}(\underline{x})\rangle$, but as in the frequency-independent case, since we are deriving what is the macroscopic equation this is simply $W(\langle\underline{x}\rangle)$, that is, $W(\underline{x})$. The comments below Eq. (81) also hold here, and substitution of the order $s$ term in Eq. (19) does indeed give $dx_{k}/d\tau=A_{k}(\underline{x})$, where $A_{k}(\underline{x})$ is given by Eq. (86).

Appendix B Using results from the Moran model to calculate fixation quantities in the SLVC

In this appendix we will use results on fixation probabilities and times in the Moran model to calculate equivalent quantities in the SLVC model. In the first section, we will consider the neutral model. Here we will calculate the probability of fixation of any species and the probability of any series of extinctions, as well as the unconditional mean time until each successive extinction. In the following two sections we will derive results for the fixation of alleles in the frequency-independent model and frequency-dependent model respectively. In particular, in these later two cases we shall calculate the fixation probability of alleles in degenerate scenarios in which there are $M$ species but only two selection strengths (frequency-independent selection) or two competition regimes (frequency-dependent selection).

B.1 The neutral models

We first calculate the probability of a particular sequence of extinctions. Consider an $M$-allele single locus neutral haploid Moran model with the frequency of each allele denoted $z_{a}$, $a=1,\ldots M-1$. We wish to calculate the fixation probability of a particular allele, say allele $b$, with frequency $z_{b}$. We begin by noting that if we are only interested in the dynamics of a single allele, we can group together the remaining alleles and treat these as a single type. The frequency of the remaining type, which we shall denote $z^{\prime}_{r}$ is then simply given by $z^{\prime}_{r}=1-z_{b}$. In a two allele system, the fixation probability of an allele is simply equal to its relative initial frequency in the population Crow and Kimura (1970);

[TABLE]

This is the probability that allele $b$ fixates, or equivalently, that $b$ does not go to extinction.

We next consider the dynamics within the subpopulation of frequency $z^{\prime}_{r}$. We ask what is the fixation probability of an allele $c$ (of frequency $z^{(0)}_{c}$ in the global population) within the subpopulation. This is given by

[TABLE]

The probability that $c$ fixates first within the subpopulation, followed by $b$ fixating in the global population is then

[TABLE]

Equivalently, this is the probability that allele $c$ is the final allele to become extinct, while allele $b$ does not go to extinction.

Iterating this argument over successive subpopulations, the probability of a given sequence $\underline{S}=(\alpha_{M},\alpha_{M-1},\ldots,\alpha_{2})$ of extinctions is Baxter et al. (2007)

[TABLE]

where $\alpha_{M}$ is the first to go extinct, $\alpha_{M-1}$ the second,…, until only allele $\alpha_{1}$ remains. Our final step is to simply transform this into a function of the original $\bm{x}^{(0)}$ initial condition variables using Eq. (34). This leads to Eq. (63) of the main text.

We next calculate the unconditional mean time to fixation of the $r^{\rm th}$ allele. In the Moran model with the frequency of each allele denoted by $z_{a}$, this is given by $\langle T_{(r)}^{\mathrm{(M)}}(\bm{z}^{(0)})\rangle$ Baxter et al. (2007)

[TABLE]

Here the summation over $\bm{\alpha}=\{\alpha_{1},\alpha_{2},\ldots,\alpha_{a}\}$ denotes summation over all possible subsets of the set of positive integers $\{1,2,\ldots,M\}$ with $a$ entries. For instance, for $M=3$ and $a=2$

[TABLE]

We also note that the result Eq. (87) differs slightly from that given in Constable and McKane (2015) in that we have stated it in natural time units of the Moran model (that is in $t^{(\mathrm{M})}$, see Table 1), rather than rescaled time, which accounts for the $N^{2}$ prefactor in our description.

Finally we require that Eq. (87) is given in units and variables appropriate for the untransformed SLVC model. Using Table 1, we can transform each of the initial conditions in Eq. (87) into their equivalent values in the SLVC model formulation (see also Eq. (34)) as well as rescaling time into the natural units of the Moran model $t^{(\mathrm{LV})}$. Recalling that the average population size in the SLVC model at carrying capacity is denoted $\langle N^{(\mathrm{LV})}\rangle$, we find that the time to fixation in the SLVC model in its natural units is given by Eq. (LABEL:eq_fixtime_SLVC_neutral).

B.2 The models with frequency-independent selection

We begin by recalling the dynamics of the two-allele Moran model, in which allele $1$ has fitness $W_{1}=1+s\rho_{1}$ and allele $2$ has fitness $W_{2}=1+s\rho_{2}$, that is, the $M=2$ version of the model described in Sec. II.2.2. The drift and diffusion terms for this system are respectively

[TABLE]

to the order we are working in $s$. The fixation probability of allele $1$ can be obtained by solving the backward FPE Gardiner (2009); Constable (2015), and is given by

[TABLE]

The fixation probability of allele $2$ is then $P_{\mathrm{Fix}}^{\{2\}}=1-P_{\mathrm{Fix}}^{\{1\}}$.

We now move on to the $M$-allele Moran model, but where the $M$ alleles are acted on by only two distinct selection pressures, $\bar{\rho}_{1}$ and $\bar{\rho}_{2}$. Suppose that $\theta$ of the alleles are acted upon by selection pressure $\bar{\rho}_{1}$. We choose to label these $1,\ldots,\theta$. Then those labelled $\theta+1,\ldots,M$ are acted under selection pressure $\bar{\rho}_{2}$, that is,

[TABLE]

Therefore the first $\theta$ alleles can be said to constitute subpopulation $1$, while final $M-\theta$ alleles can be designated subpopulation $2$. Since the dynamics within each subpopulation are neutral, we can say that the probability of each subpopulation fixating is simply given by $P_{\mathrm{Fix}}^{\{1\}}$ and $P_{\mathrm{Fix}}^{\{2\}}$ (see Eq (90)). The reason that this can be done, is that different alleles in each subpopulation only differ by the labels given to them; we can therefore choose to label them as only belonging to a particular subpopulation, without changing the dynamics. Meanwhile, again since the subpopulations are neutral, the probability that each allele fixates within its respective subpopulation is simply equal to its initial frequency within the subpopulation. Therefore, the probability that any allele $\alpha$ fixates is equal to the product of the probability that it fixates within its subpopulation and the probability that its subpopulation fixates;

[TABLE]

where the initial frequency of the $M^{\rm th}$ allele, $z^{(0)}_{M}$ is understood to be given by $1-\sum^{M-1}_{b=1}z^{(0)}_{b}$. Notice that unlike in the neutral case, we cannot determine the probability of allele extinctions in a particular order. This is because we can only break this non-neutral system down into two distinct subpopulations, whereas in the neutral case we could break the system down into any focal allele, plus the remainder of the population.

The final task is to write the above equations in terms of the original SLVC model variables and parameters. Note that the mapping for the parameter $\rho_{a}$ depends on a combination of the parameters in terms of the original SLVC model variables (see Table 1). In order to make use of the calculation above, we require that each $\rho_{\alpha}$ for $\alpha\leq\theta$ is identical, however this does not mean that the parameters in the SLVC model need be identical; we merely require that $b_{0}\beta_{\alpha}-d_{0}\delta_{\alpha}-(b_{0}-d_{0})\gamma_{\alpha M}=\bar{\rho}_{1}$ for $\alpha\leq\theta$ and $b_{0}\beta_{\alpha}-d_{0}\delta_{\alpha}-(b_{0}-d_{0})\gamma_{\alpha M}=\bar{\rho}_{2}$ for $\alpha>\theta$. This leads to the equivalent expressions given by Eq. (65) for the model for the model under consideration, which can be seen clearly in the parameters used to generate figure 4 (see Appendix C). The fixation probabilities in SLVC model notation are then given by Eq. (66) of the main text.

B.3 The models with frequency-dependent selection

Let us begin by recalling the dynamics of the two-strategy model, that is, the $M=2$ version of the model described in Sec. II.2.3:

[TABLE]

to the order we are working in $s$. The fixation probability of allele $1$ can be obtained by solving the backward FPE Gardiner (2009); Constable (2015). Defining the function

[TABLE]

one finds

[TABLE]

where

[TABLE]

Here erfc and erfi are the complementary and imaginary error functions respectively Abramowitz and Stegun (1965); Erdélyi (1953). Again, the fixation probability of strategy $2$ is $P_{\mathrm{Fix}}^{\{2\}}=1-P_{\mathrm{Fix}}^{\{1\}}$.

In a similar fashion to the frequency-independent case, we now envisage a scenario in which we have multiple strategies, but the payoff matrix can be partitioned such that;

[TABLE]

Let the first $\theta$ alleles be subpopulation $1$, while final $M-\theta$ alleles can be designated population two. Since the dynamics within each subpopulation are neutral, we can say that the probability of each subpopulation fixating is simply given by $P_{\mathrm{Fix}}^{\{1\}}$ (see Eq (95)) and $P_{\mathrm{Fix}}^{\{2\}}=1-P_{\mathrm{Fix}}^{\{1\}}$. Meanwhile the probability that each allele fixates within its respective subpopulation is simply equal to its initial frequency within the subpopulation. Therefore, the probability that any allele $\alpha$ fixates is equal to the product of the probability that it fixates within its subpopulation and the probability that its subpopulation fixates:

[TABLE]

where

[TABLE]

and where

[TABLE]

Here, by $l(1)$ and $l(0)$ we mean the value of $l$ when $\sum_{b=1}^{\theta}z^{(0)}_{b}=1$ and $\sum_{b=1}^{\theta}z^{(0)}_{b}=0$, respectively. Finally, writing this as a function of the SLVC variables and parameters, we obtain Eq. (LABEL:eq_fixprob_SLVC_dep) of the main text.

Appendix C Parameters used in figures

In Figure 1, the parameters used for the SLVC are

[TABLE]

In Figure 2, the parameters used for the SLVC are

[TABLE]

and the parameters used for the Moran model are

[TABLE]

The initial conditions used in the SLVC model and Moran model are respectively

[TABLE]

In Figure 3, the parameters used for the SLVC are

[TABLE]

[TABLE]

[TABLE]

and the parameters used for the Moran model are

[TABLE]

[TABLE]

Note that these choices of parameters make the models mappable to each other (see Table 1). Also note that although the values of the parameters $\gamma_{Ma}$ are needed to specify the dynamics of the SLVC simulation, they are not required for the mapping, and thus the same mapping between the SLVC model and the Moran model with frequency-independent selection holds for any order one choice for the parameters $\gamma_{Ma}$. The initial conditions used in the SLVC model and Moran model are respectively

[TABLE]

In Figure 4, the parameters used for the SLVC are

[TABLE]

[TABLE]

[TABLE]

and the parameters used for the Moran model are

[TABLE]

[TABLE]

Note that these choices of parameters make the models mappable to each other (see Table 1). The initial conditions used in the SLVC model and Moran model are respectively

[TABLE]

In Figure 5, the parameters used for the SLVC are

[TABLE]

[TABLE]

[TABLE]

and the parameters used for the Moran model are

[TABLE]

[TABLE]

Note that these choices of parameters make the models mappable to each other (see Table 1). The initial conditions used in the SLVC model and Moran model are respectively

[TABLE]