Fixation probabilities in populations under demographic fluctuations

We study the fixation probability of a mutant type when introduced into a resident population. We implement a stochastic competitive Lotka–Volterra model with two types and intra- and interspecific competition. The model further allows for stochastically varying population sizes. The competition coefficients are interpreted in terms of inverse payoffs emerging from an evolutionary game. Since our study focuses on the impact of the competition values, we assume the same net growth rate for both types. In this general framework, we derive a formula for the fixation probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document}φ of the mutant type under weak selection. We find that the most important parameter deciding over the invasion success of the mutant is its death rate due to competition with the resident. Furthermore, we compare our approximation to results obtained by implementing population size changes deterministically in order to explore the parameter regime of validity of our method. Finally, we put our formula in the context of classical evolutionary game theory and observe similarities and differences to the results obtained in that constant population size setting.


Introduction
The evolutionary dynamics of a mutant strain in a resident population is a well-studied topic in the field of population dynamics. Results concerning the fixation probability, the average fixation time or coexistence behavior can be applied in various biological fields, e.g. population genetics, bacterial evolution, viral dynamics or cancer initiation (Nowak 2006;Altrock et al. 2015;Ewens 2004). While the first theoretical analysis of such processes relied on deterministic differential equations, over the course of time more detailed models were studied describing the stochasticity of microscopic processes on the individual level. These individual based models can be approximated and studied by the replicator equation (in the large population size limit), by branching processes (in case of a small number of invading mutants) or by multi-type birth-death processes (Nowak 2006;Hofbauer and Sigmund 1998;Sandholm 2010;Haccou et al. 2005;Ewens 2004). However, the ecological evolution of the entire population is mostly neglected in these kinds of models and a constant population size is assumed instead. On the other hand, in population genetics and theoretical ecology, studies focused more on the effect that population dynamics have on the fixation probability rather than the concrete interaction mechanisms between the mutant and wild-type individuals (Ewens 1967;Kimura and Ohta 1974;Otto and Whitlock 1997). More recently, researchers started investigating models connecting the stochastic interaction between individuals and stochastic population dynamics from a theoretical point of view (Lambert 2005(Lambert , 2006Champagnat and Lambert 2007;Parsons and Quince 2007a, b;Melbinger et al. 2010;Cremer et al. 2011;Gabel et al. 2013;Constable et al. 2016;Chotibut and Nelson 2017). For a historical overview on the calculation of fixation probabilities, see Patwa and Wahl (2008).
To our knowledge, the first analytical approximation of fixation probabilities under stochastically varying population sizes is due to Lambert (2005Lambert ( , 2006. In these papers, the author analyzes models of interacting species by considering the corresponding diffusion equations under the constraint of weak selection. Going one step back on the descriptive scale and analyzing the Kolmogorov forward equation instead of its diffusion approximation, Champagnat and Lambert study the effect of various model parameters on the fixation probabilities and extend the previous results (Champagnat and Lambert 2007). In parallel to these studies, Parsons and Quince examined the effect of variable growth rates on the fixation probability and mean fixation time in a two species system with stochastically varying population size (Parsons and Quince 2007a, b). These results were later complemented and refined in Parsons et al. (2010). Instead of focusing on variable growth rates, in this paper we concentrate on the effect of variable competition coefficients on the fixation probability.
The model we will work with was introduced in Huang et al. (2015). It is a generalized stochastic Lotka-Volterra-model which connects an evolutionary game with the competition coefficients of the model. Individuals reproduce at constant rates and die spontaneously or based on competition within and between species. This leads to stochastically induced demographic fluctuations driven by interactions within the population. We focus on two types only and our goal is to calculate the probability that a mutant takes over such a population of changing size.
Recently, further models have been studied which connected game theoretical dynamics with exogenous population growth (Melbinger et al. 2010;Tao and Cressman 2007;Traulsen et al. 2012;Cremer et al. 2011;McAvoy et al. 2018). For instance, Ashcroft et al. (2017) consider a model with deterministic cell growth defined by a power law and stochastic species interactions derived from an evolutionary game. The authors rely on simulation results suggesting that the evolutionary outcome not only depends on the game played by the species, but also on the growth exponent of the power law governing the population growth.
In the context of public goods dynamics (Constable et al. 2016) study the invasion probability of producers and non-producers of the public good under varying population sizes. Using a time-scale separation under a weak selection approximation, they find that producers can successfully invade a colony of non-producers even though they have a lower fitness than the resident type.
The present paper is structured as follows: In Sect. 2 we describe the generalized Lotka-Volterra-model and restate some basic properties of the system, which were already described in Huang et al. (2015). In Sect. 3 we apply tools developed by Lambert (2006) in order to derive a formula for the fixation probability in the weak selection limit. This allows us to interpret the impact of the competition coefficients separately. Furthermore, we compare the results for various competition matrices induced by different games with each other, i.e. the differences between coordination, coexistence and dominance games. Finally, in Sect. 4 we compare our results to similar models derived in the physics and population genetics literature dealing with deterministic ecological evolution of the model. Furthermore, we examine the fixation probability of a single mutant in a wild-type population, which allows us to compare our findings with those obtained in classical evolutionary game theory, i.e. in finite but fixed population sizes.

Model
The model we consider is a competitive Lotka-Volterra system consisting of two types, the mutant X and the wild-type Y . We assume a well-mixed population, i.e. dynamics do not depend on the spatial arrangement of individuals, and a discrete state space describing the number of individuals of the two types, X and Y .
The evolution of the system is described by birth, death and competition processes, which we assume can be written in terms of chemical reactions. Each individual of the two types can reproduce or die independently of the other individuals. This leads to four reactions for the birth-death-processes, Here, β X , γ X and β Y , γ Y denote the birth and death rates of the mutant and the wildtype, respectively. Additionally, each individual competes with other individuals and might die due to this process. These reactions occur at the rates where M controls the total population size in stationarity. Later on, we interpret the competition rates as inverse payoffs of an evolutionary two-player game with payoff matrix This interpretation of the competition processes and a descriptive study of the stochastic competitive Lotka-Volterra system as well as a stability analysis of the stationary points of the corresponding deterministic system was performed by Huang et al. (2015). This setup has the advantage that the average size of a monomorphic population reflects the payoffs. For example, focus on a prisoner's dilemma, where a population of cooperators would be better off than a population of defectors, but selection acts in favor of defectors. A population of cooperators would be larger than a population of defectors, reflecting the fitness values within the population. Letting X (t) and Y (t) be the number of mutant and wild-type individuals at time t, respectively, the differential equations of the deterministic model that emerges in the large population limit read For a > c and d > b as well as for a < c and d < b, these equations have an internal stationary point where both species exist. It is given by Its stability depends on whether a coordination (a > c and b < d) or a coexistence game (a < c and b > d) is played. Additionally, we see that M indeed characterizes the scale of the total population size in equilibrium. In the following, we will work with the rescaled variables x(t) = X (t) M and y(t) = Y (t) M while the fraction of mutants in the population is given by p = x x+y . We denote the steady state of this value by Our goal is to extend the analysis of this particular system by approximating the fixation probability of the mutant type X in a population of Y individuals. The techniques we use rely on the theory of stochastic diffusions, see e.g. Ewens (2004). Hence, we will work with the diffusion approximation of the above system; for a detailed derivation see "Appendix A". We find the stochastic differential equations where W 1 and W 2 are two independent, one-dimensional Brownian motions. The stochastic integrals are interpreted in the sense of Itô (van Kampen 1997;Gardiner 2004). The solution of this system of differential equations is a two-dimensional Markov processes with infinitesimal generator given by (see "Appendix A" or Kallenberg 2002, Chapter 21) We now proceed in deriving the fixation probability of the mutant type X .

Fixation probabilities
The main result of this paper is the approximation of the probability of a mutant strain to reach fixation in a resident population of randomly fluctuating size under weak selection. Note first that due to the competition processes neither of the two species is able to go to ∞ and hence each of them will die out at a (finite) random time (Lambert 2006). We define fixation of the mutant X as follows (see also Lambert 2006, Definition 1.1 for a similar definition): Definition 1 (Fixation) Species X becomes fixed if for some t ≥ 0 we have y(t) = 0 and x(t) > 0.
In order to quantify the fixation probability we make use of the generator description of the model. Let ϕ(x 0 , y 0 ) be the fixation probability of species X if the initial type-frequencies are x 0 and y 0 . Then diffusion theory, see also Ewens (2004) and Gardiner (2004) or "Appendix B", implies that ϕ solves In order to solve this partial differential equation we first do a parameter transformation to the coordinates p = x x+y and z = x + y = X +Y M , the fraction of X -individuals in the population and the whole rescaled population size, respectively. Given the same birth and death rates for both species, i.e. β X = β Y = β and γ X = γ Y = γ , and noting that the generator transforms to (the detailed calculations are given in "Appendix C") From now on, we drop the indices of p 0 and z 0 since the fixation probability always depends on the corresponding initial values. Our goal is to approximate the solution of Eq. (8). Therefore, we start with the neutral setting which forms the basis of the subsequent calculations.

Neutral model
In formal terms, a neutral setting is given when individuals are exchangeable under labelling which in our case is equivalent to choosing a constant competition matrix, i.e. a = b = c = d. In this scenario the generator in Eq. (7) simplifies to Solving Gϕ neu ( p, z) = 0 with boundary conditions ϕ neu (0, z) = 0 and ϕ neu (1, z) = 1 for z > 0 we obtain ϕ neu ( p, z) = p, the standard fixation probability of a mutant in an evolutionary process without selection.

Fixation probability under weak selection
Based on the result of the neutral setting we approximate the fixation probability ϕ( p, z) in the case of weak selection. In our model, this translates to the coefficients of the competition matrix being similar. To be more concrete we need the following conditions 1. Condition (i) demands that the parameters d and b are close together such that their squared difference is negligible. Due to conditions (ii) and (iii) it follows that also a and c are close to d and thus b which then implies that the whole competition matrix is just a small perturbation of the neutral matrix in which all entries are the same. Finally, condition (iv) is needed to ensure that there indeed is an internal fixed point in the deterministic system which is close to 1/2. For the case that there exists no internal fixed point we refer to Sect. 3.4.
We note again, that we claim that the product of two values which are already small is negligible, the single terms separately however are important for our approximation formula. We see in the derivation of the following theorem that the condition on the equilibrium p * is necessary for the separation of the variables p and z reflecting a separation of the ecological and the evolutionary variable. For a extended discussion on the choice of fixed points in these kind of models see also Czuppon and Gokhale (2018). In the following we will make use of asymptotic notation, i.e.
We now state our main result.

Theorem 1 (Fixation probability under weak selection) Under conditions (i)-(iv) the solution of Eq. (8) can be written as
where ψ(z) is independent of the initial frequency of mutants p and solves Remark 1 We note that this is basically a linearization around the neutral fixation probability and reduces to the neutral model if all payoff coefficients are equal, i.e. ϕ( p, z) = p due to 1 − d b = 0.
The proof of the Theorem is given in "Appendix D". The idea is to apply the infinitesimal generator G from Eq. (7) to the formula stated in Eq. (9). Inserting conditions (i)-(iv) then gives the result. It seems remarkable that the initial population size does not affect the qualitative behaviour of the fixation probability, i.e. ψ(z) does not change signs as we will see later. However, the initial frequency of the mutant compared to the internal steady state and the payoffs b and d can change the sign of the first order effect under weak selection. An interesting application is to consider fixation out of the neighbourhood of the internal steady state. Precisely at that point, we find ϕ( p * , z) = p * under weak selection. This is also implicit in the formulas derived in Lambert (2006), Champagnat and Lambert (2007), and Constable and McKane (2015) but has not been stated explicitly. For p = p * + ε, we have which implies that the fixation probability out of a neighborhood of a stable steady state p * (d < b) is smaller than neutral for positive deviations in p and larger than neutral for negative deviations in p. On the other hand, for an unstable steady state p * (d > b), the fixation probability out of the neighborhood is larger than neutral for positive deviations in p and smaller than neutral for negative deviations in p. For a detailed study of fixation probabilities when leaving the deterministic steady state see also Park and Traulsen (2017). Next, we investigate different competition parameter constellations, i.e. conditions on the evolutionary game. We consider the following cases (a) coexistence game a < c and b > d; The cases (a) and (b) allow for a mixed steady state in the deterministic model given in Eq.
(3). For coexistence games, this internal equilibrium is stable whereas for coordination games it is unstable, see for instance (Huang et al. 2015).
A qualitatively different picture arises in case (c). Here, either type X or type Y strictly dominates the other species in a game theoretic sense. This implies that the deterministic model only allows for single species equilibria where the stationary point of the dominant (inferior) type is stable (unstable). Thus, Theorem 1 does not hold in this case since p * does not tend to 1 2 . In fact, p * does not even exist. Instead we will replace condition (iv) by an adapted version which then gives a similar approximation, see Eq. (14) in Sect. 3.4 below.

Coexistence and coordination games
In this section, we compare the resulting fixation probabilities in a coexistence and coordination game with the neutral fixation probability ϕ neu ( p, z) = p. In order to do so we need a Lemma characterizing the impact of the initial population size which we prove in "Appendix E".
Remark 2 In fact the function ψ is a growing function in z as can be seen in Fig. 1. This can be explained as follows: since the only difference between the strains are the corresponding competition parameters, regions where these parameters play a dominant role should have a larger impact on the fixation probability. Translating this to initial population sizes yields that larger initial population sizes, being affected mostly by competition processes, also have a larger influence on the evolutionary outcome and thus a larger value for ψ. Now, we can state some immediate consequences of the fixation probability which follow from Eq. (9).
Corollary 1 (Impact of competition parameters) Given the assumptions of Theorem 1 we find the following: 1. For arbitrary a, c > 0 and p < p * we have that ϕ > ϕ neu iff b > d.

The probability of fixation is an increasing function in the competition parameter a.
Proof Part 1. follows immediately by comparing ϕ neu with ϕ from Eq. (9) and Lemma 1. For part 2., we differentiate the representation of the fixation probability from Eq. (9) with respect to a which gives The last inequality holds for all choices of the parameter values which finishes the proof.
Remark 3 The first statement of the Corollary has the obvious implication that for a mutant to invade a resident population it is important to perform well against the wildtype. This also implies that species with lower single species equilibria (i.e. a < d) can have a higher chance of fixating than neutral. This can end up in a decrease of the overall population size. But still, as the second part of the Corollary shows, species with higher single-species equilibria also have a higher chance to reach fixation.
Before turning to dominance games we take a brief look at some special cases in the context of coexistence and coordination games.

Symmetric and asymmetric games
Here, we assume that a = d and distinguish between symmetric games, i.e. b = c and asymmetric games, b = c. In the case of symmetric games we have a c = d b and thus p * = 1/2. Therefore the fixation probability in Eq. (9) simplifies to Additionally, note that in this case we do not need assumption (iv) for the solution of the generator equation in (8). For an illustration of the fixation probability with some simulated data points, see Fig. 2. For coexistence games (b > d) the fixation probability lies above the neutral line while for coordination games (b < d) the fixation probability is lower. Choosing b closer to d improves the analytical prediction due to the weak selection approximation in conditions (i)-(iii). For asymmetric games, i.e. we still assume a = d but now b = c, we obtain similar results. In this case, the fixation probability is given by The condition for invasion, i.e. for ϕ asym > ϕ neu to hold, is again given by b > d. Thus, for a mutant to reach fixation in the resident population, it is primarily important to have a high payoff when playing against the resident, i.e. the more abundant type. However, when comparing ϕ sym and ϕ asym we see that here the parameter c does play a role. To be more precise, whenever b > c we have ϕ asym > ϕ sym . This condition resembles the shifting of the internal equilibrium of the deterministic system towards the mutant-axis, i.e.

Dominance game
In contrast to coexistence and coordination games the dominance game does not have an internal equilibrium in its deterministic counterpart. This is due to one strategy strictly dominating the other. In terms of parameters this means that either a > c and b > d or a < c and b < d hold. As already mentioned, for the analysis of the fixation probability it does not make sense to assume condition (iv) which states that the internal equilibrium should be close to 1 2 . Hence, we can already infer that the analytical solution will not intersect with the neutral line due to one strategy being favored independently of its frequency (Fig. 4). For the calculation of ϕ dom we still assume conditions (i)-(iii) but instead of condition (iv) we need the following: This is plausible since the fraction in the last term should approximate −1 when considering a dominance game under weak selection, i.e. either c > a and d > b or c < a and d < b.

Theorem 2 Under conditions (i)-(iii) and (v), we find
where ψ(z) satisfies The proof is an imitation of the proof of Theorem 1 and therefore spared out. We see that indeed ϕ dom is always larger or smaller than ϕ neu dependent on b being larger or smaller than d, respectively. This finding is rather trivial, since we consider a dominance game and b > d ensures the mutant being advantageous. More importantly, Eq. (14) allows to calculate the first order approximation of the neutral result.

Discussion
Now that we have derived the fixation probabilities in all possible scenarios captured by our model, we can compare them to previously derived results. Since our model takes into account stochastic fluctuations in population size, it is especially interesting to see how that affects the fixation probability in comparison to deterministically varying or constant population sizes. Concerning constant population sizes we also compare our results to the 1 3 -rule which deals with the invasion of a single mutant into a wild-type population.

Comparison to models with deterministically varying population size
Fixation probabilities under deterministically changing population sizes have been extensively studied in the population genetics literature (Uecker and Hermisson 2011;Otto and Whitlock 1997;Kimura and Ohta 1974;Waxman 2011;Wahl and Gerrish 2001). It is possible to, at least partially, compare our results to these models. In order to do so, we consider the case of a dominating trait with frequency independent selection, i.e. we are in the context of Theorem 2 with parameters being equal along the rows of the competition matrix, i.e. a = b and c = d. This transforms the competition being only dependent on the population density and not on the population composition. Then typically s = a − c represents the fitness advantage of a mutant when compared to a wild-type resulting in a higher replacement (or reproduction) rate of wild-type individuals. This effect cannot easily be translated to our mechanistic model in (1) and (2). In our model selection happens due to competition and thus affects the death rates of individuals, see (Baar et al. 2017) for a similar model in the adaptive dynamics framework. To account for that we should set s = 1 d − 1 b which is already part of the first order term in our derived approximations in equations (9) and (11), see also Lambert (2006) for a similar definition of selection in a similar model.
Using the selection coefficient s and the population size M, a classical assumption in population genetics is Ms 1 which is violated in our setting since our approximation is applicable when Ms < 1 which emerges from the selection interaction modeled by the microscopic competition reactions stated in (2). Furthermore, once an individual dies due to a competition event it does not get replaced by another one, allowing for fluctuating population sizes. This is, as far as we can judge, different to models with varying population size in existing population genetics literature where the population size dynamics are often separated from the evolutionary dynamics. These two differences prevent a rigorous comparison between results obtained in the population genetics context and our microscopic model.
Employing a mechanistic individual based model as we do has both advantages and disadvantages. As a drawback, one could argue that externally driven fluctuations in population size, for example seasonal fluctuations or bottleneck experiments as analyzed in Uecker andHermisson (2011), andWahl andGerrish (2001), cannot be dealt with-at least not in our studied model. On the other hand the microscopic model allows for concrete interpretation of the selection dynamics and allows us to explore the stochastic effects in small populations in more detail. That is, as opposed to a general selection value s affecting the birth rate, we can study the specific impact of varying the different parameters corresponding to single reactions much like it was performed in Lambert (2006), Champagnat and Lambert (2007), Parsons and Quince (2007a, b), Parsons et al. (2010), and Czuppon and Gokhale (2018). This microscopic approach has recently received some attention, see Doebeli et al. (2017).
We now proceed to compare our approximation with those obtained in the literature considering similar scenarios. Otto and Whitlock (1997), Kimura and Ohta (1974) In Otto and Whitlock (1997) the authors derive a stochastic diffusion approximation for populations consisting of two types, one strictly dominant, under frequencyindependent selection and deterministically varying population sizes. In case of logistic growth they find in their Eq. (11) for the fixation probability of a single mutant

Deterministic logistic growth-
where s is the selection advantage, K the carrying capacity and r the logistic growth parameter. These parameters can be translated to our setting. As the selection strength we set s = 1 d − 1 b which is the selective advantage of the wild-type in the frequency-independent case, i.e. a = b and c = d. The carrying capacity is set to K = (β −γ )a M and r = β − γ is the effective growth. The equation derived by Otto and Whitlock is valid for s, r < 0.1 and large population sizes z M. Since neither r < 0.1 nor z M 1 are met in our simulations it is not surprising to see deviations between our approximation and (16), see also Fig. 5.
Further, the authors mention that when the conditions s, r < 0.1 do not hold, one can use the approximation derived by Kimura and Ohta (1974). There, the Kolmogorov backward equation is used and the result they obtain is given by (see equation (5) in Kimura and Ohta 1974 or (A1) in Otto and Whitlock 1997) where Q is given by (see Eq. (A8) in Otto and Whitlock 1997) However, for the derivation of the equations to hold, one needs Qs 1 which is not satisfied in our setting such that again it is not unexpected to see differences between the two approximations, see Fig. 5.
Due to the underlying assumptions, these approximations only work in case of an advantageous mutant, i.e. s > 0. Therefore, in the figures we only compute the corresponding values in case of a = b > c = d.

Deterministic logistic growth-Uecker and Hermisson (2011)
The model developed by Uecker and Hermisson (2011) is different to the previous one since it strictly distinguishes between ecological and evolutionary steps. On the ecological scale both types are the same while due to evolutionary updating the fitter type on average replaces the other type more often. Since selective advantage is again frequency-independent we are again in the case of dominance. Their result for logistic growth and constant selective advantage, equation (24) in Uecker and Hermisson (2011) reads where s is the selective advantage of the mutant, r the logistic growth rate, b the birth rate of individuals, ξ the neutral evolutionary replacement rate and K the carrying capacity. Translating our modeling parameters into this setting we obtain s = 1 d − 1 b , r = β −γ , b = β, ξ = 0 and K = ra M. We note, that the choice of ξ is to some extent arbitrary. We chose to set it to 0 since we do not have any comparable mechanism in our model. range covered by our approximation such that we find a disagreement between the approximations, see Fig. 5.
In contrast to the derived approximations before, we can use the method from Uecker and Hermisson (2011) to solve for varying initial frequencies. Even though the branching process assumption only makes sense for an invading and initially rare mutant, the formula can be applied to larger values of p. Their equation (14) together with (24) then give A visual comparison of this result with our approximation is given in Fig. 8. As we see there, the branching process approach overestimates the simulation results. This discrepancy can again be attributed to ϕ U H (1/(z M))M 10 not being satisfied.
However, the branching process approach is, as the approximations by Otto and Whitlock and Kimura and Ohta, limited to the case of beneficial mutations, i.e. s > 0. We therefore compute only the case a = b > c = d for this approximation in the figures below.
Remark 4 It might seem odd that both branching process estimates obtained by Otto/Whitlock and Uecker/Hermisson overestimate the actual fixation probability. Usually, branching process approaches tend to underestimate fixation probabilities since they do not account for the absorbing state at p = 1. However, the overestimation in these cases comes from the distinct implementation of the selection processes. While in our model selection acts implicitly through the competition processes, in the mentioned papers it is explicitly added to the birth rate of the mutant. The complicated intertwining of selection and population ecology in our model causes difficulties in computing the fixation probabilities in the same vein as suggested in the branching process approaches. On the other hand the added selective advantage to the birth rate causes the mutant to have a higher reproductive rate independent of the population composition and size which does not hold either in our case. Hence the overestimation of the results obtained in Otto and Whitlock (1997) and Uecker and Hermisson (2011).

Nearly constant population size-Constable and McKane (2015)
We now want to compare our results to a model with nearly constant population size, i.e. basically constant population size with small fluctuations around the ecological equilibrium. Therefore, we use a time-scale separation argument, developed by Constable and McKane (2015) in the case of Lotka-Volterra dynamics.
In order to do so, we distinguish between the fast ecological scale and the slow evolutionary process. This is reasonable since due to our choice of reaction rates in Eqs. (1) and (2) we see that the population quickly equilibrates. More precisely, looking at Eq. (3) we see that for a initial population size of order O(1) the process first grows until it reaches size O(M) since until then the competition terms, which are of order O(M −1 ), can be neglected. A similar reasoning shows that for initial population sizes of order O(M α ) for α > 1 the population decreases to order O(M) since the linear birth and death terms are negligible, at least when compared to the competition rates. Thus, competition is the key driver in this case. This reasoning is true as long as M is sufficiently large such that the birth-death and competition processes can be separated. As we see in Fig. 6, choosing M = 100 is already large enough to observe this time-scale separation.
For the formal approximation of the population equilibrium we use the projection technique developed in McKane (2015, 2017) which maps the twodimensional dynamics onto a one-dimensional manifold. Briefly, the method works in the context of weak selection, i.e. in our case small differences of the competition parameters. Additionally, the neutral model needs to exhibit a stable center manifold which holds for general Lotka-Volterra dynamics. Under these assumptions one ends up with a stochastic diffusion along this manifold which can then be analyzed by stochastic diffusion techniques, i.e. making use of the scale function of the diffusion process (for details see Ewens 2004, Chapter 4.3). Translating the results obtained in Constable and McKane (2015) to our model we obtain the following diffusion (the derivation is given in "Appendix G"-see also equation (14) in Constable and McKane (2015))  . 6 A stochastic simulation of the dynamics in a coexistence game. One can see that the trajectories quickly enter the neighborhood of the deterministic equilibria of the system, displayed by the dashed lines. Further, once type Y goes extinct, the blue trajectory approximates the monomorphic equilibrium of X visualized by the blue dotted line. This shows that a time-scale separation is indeed a reasonable assumption in a system under Lotka-Volterra dynamics. The initial condition is X 0 = 100 and Y 0 = 100. We now compare our results for the fixation probability with those obtained by using (20). As we can see in Fig. 7 the two results are very similar in the cases of coexistence and coordination game dynamics. Looking at the subfigures of Fig. 7, we only observe small differences between the two approaches. The largest deviation between the two models arises in the bistable case in subfigure (a) where our result is slightly more accurate. In that case we consider initial population sizes smaller (z = 0.2) than the carrying capacity on the center manifold (z = 0.5). Interestingly this effect is reversed when starting above the center manifold, see subfigure (b). A possible explanation might be an overestimation of stochastic fluctuations by our approach in (b) and an underestimation of stochastic effects by Eq. (20) in (a). The overestimation of stochasticity would also explain our approximation often being smaller than the observed fixation probabilities for small initial mutant populations (∼ p smaller than 0.3). Fluctuations can drive the rare mutant type to extinction while deterministic dynamics pull it to the center manifold where the reasoning derived in Constable and McKane (2015), our Eq. (20), is valid. If deterministic dynamics are strong enough (or stochasticity weak enough) initial fluctuations can be neglected which corresponds to the approach taken in Constable and McKane (2015).
In case of a dominant strain, we observe some differences. As we can see in Fig. 8a our result captures the simulation results slightly better than the approximation from Constable and McKane (2015). This difference emerges since we do not start close to the center manifold but at z = 0.1. Hence, the initial value on the center manifold does not agree or is not close enough to the initial frequency of mutants of the simulations. For techniques of how to estimate the initial value on the slow sub-system in a more sophisticated way we refer to Constable and McKane (2018) and Roberts (1989). This difference however vanishes when we increase the population size as done in subfigure (b) which also means that we decrease the impact of stochastic fluctuations. This can be seen by our result being slightly off the data points whereas the model derived in Constable and McKane (2015) neglecting the demographic fluctuations closely follows the simulation results.
Summarizing the comparison between the different models we observe that for lower population sizes (M ∼ 100) where stochastic fluctuations have a strong impact on the model dynamics, our approximation which explicitly accounts for these fluctuations, provides the best fit. For intermediate population sizes (M ∼ 1000) our approximation seems to overestimate the effect of demographic stochasticity and thus methods like the projection of Constable and McKane describe the evolutionary trajectories best. However both, our approximation as well as the time-scale separation work with very weak selection strengths (s = 1 − d b ) such that Ms = O(1). Other scalings of Ms yield different limiting models such that other techniques are necessary for analyzing the fixation probability. Therefore, e.g. Uecker and Hermisson use branching processes to study a system with a larger population size (or selection strength) shifting Ms to larger values. Their results work best for parameter choices which yield Ms > 10 and provide better fits than the previously mentioned approximations.
Remark 5 In Fig. 8b we chose M = 1, 000 which is in population genetics terms still a small population. Unfortunately, for values of M larger than that we ran into numerical problems (unstable convergence of the method) solving the differential Eq. (15) for ψ.

Fixation of a single mutant under frequency-dependent selection
In this section, we consider the case that p = 1 zM = 1 1+Y 0 , i.e. initially there is exactly one mutant present in the population. This is probably the most interesting scenario as seen from a biologist's perspective. In contrast to the previous sections we fix the initial fraction of mutants and instead vary the initial population size z M. As can be seen in Fig. 9, the probability of fixation is a decreasing function of the initial population sizes. This translates to the already observed fact that fixation of a mutant strain is more likely in a smaller population and therefore initially growing population than in a larger one which is initially decreasing, cf. Kimura and Ohta (1974). We can also infer this from the formula describing the fixation probability: since we are working in the weak selection limit the governing part of ϕ( p, z) is the initial frequency of mutants, here 1 zM , which is decreasing for increasing z which means that ϕ is a decreasing function in z.

Comparison to models with constant population sizes
In models with a constant population size N , the fixation probability of a mutant strain under weak selection can be translated to the location of the mixed equilibrium of the corresponding replicator equation. More specifically, in the case of bistability (coordination game) a mutant has a higher probability than neutral, i.e. 1 N , to invade a resident population if the basin of attraction of the wild-type is smaller than 1 3 . This is referred to as the 1 3 -rule first derived in Nowak et al. (2004) and later generalized in scope in Lessard and Ladret (2007); Pfaffelhuber and Wakolbinger (2018). An equivalent approximation has also been made for multiplayer games (Kurokawa and Ihara 2009;Gokhale and Traulsen 2010;Lessard 2011). However, as noted in Traulsen  (2017) the order of limits is relevant for the rule to hold. One first needs to take the weak selection limit and then perform an approximation for large population sizes N . Just by looking at our conditions for the approximation to hold, especially condition (iv) stating that p * ≈ 1 2 , we would not expect the 1 3 -law to hold. Still, we are able to draw some more general conclusions from our approximation and those obtained in Constable and McKane (2015) which suggests that a similar rule may not hold in populations with fluctuating size.
Our analytical approach relies on a stochastic diffusion equation which depends on a large population size approximation, i.e. N 1 or in our case M 1. Just after that approximation we let selection be weak which translates to the competition coefficients being similar. Thus again, it should not be surprising that the 1 3 -rule is violated in our model. Instead, the invasion probability only depends on the competition rates b and d describing the competition pressure of the species due to the resident type. However, these do not imply any properties of equilibria in the corresponding deterministic system. Even more, the choice of the parameters a and c does not affect the qualitative behavior of the fixation probability when compared to the neutral case as was already pointed out in Corollary 1. Thus, in the present model the difference between the neutral fixation probability and the probability of fixation under weak selection can be entirely described by the parameter d b and is independent of a c as can be seen in Fig. 10.
Still, our approximation of the fixation probability is in some sense similar to the approximation obtained by first taking the weak selection limit and then the population size approximation. In Nowak et al. (2004), and Lessard and Ladret (2007) where ω is the selection intensity. This has the same form as our result obtained in Theorem 1 with ω = 1 − d b . The reason for this is that even though we perform a diffusion approximation, we still let N (or rather M) be finite which results in the first term 1 zM being non-zero in contrast to being 0 for M → ∞. Additionally, our first order term also involves the deterministic equilibrium p * which is again reminiscent of the formulas derived under weak selection. Hence, our result can be seen as a connection between the two limit orderings since it involves both, ecological and evolutionary terms in the first order approximation. Thus, our approach does not belong to the class of "first population limit, then selection limit" models as for instance analyzed in Sample and Allen (2017) nor does it belong to the "first weak selection, then large population limit" class but is rather somewhere in between the two limit orderings (even though formally according to Definition 2 provided in Sample and Allen (2017) our limit belongs to the category of ωN -limits). This can also be related to the fact that our approximation holds neither in case M(1 − d/b) → ∞ (typically associated with the N ω-limit) nor for M(1 − d/b) → 0 (associated with the ωN -limit).
Moreover, this sort of breakdown of the 1 3 -rule under demographic fluctuations has also been observed by Ashcroft et al. (2017). Instead, at least in our case, the fixation behavior simply depends on the stability of the internal fixed point of the deterministic system. For a asymptotically stable fixed point the fixation probability of a single mutant is larger than neutral, while for an unstable fixed point it is smaller. This similarity between the general replicator equation, i.e. the deterministic evolutionary behavior in the fixed population size model and the fixation probability holds for both, our stochastically varying population size model, as well as for the model with constant population size derived in Constable and McKane (2015), see also Fig. 7 in Sect. 4.1. Furthermore, the similarity between the fixation probability and its corresponding replicator dynamics is generalizable to models with multiple player interactions under demographic fluctuations (Czuppon and Gokhale 2018). These observations lead us to our conjecture that the qualitative behavior of the fixation probability under weak selection and varying population sizes, at least in a mechanistic modeling approach, is largely influenced by the deterministic replicator dynamics as opposed to the stochastic birth-death dynamics which led to the 1 3 -rule. The main reason for it being the conditions on the fixed point under which ecological and evolutionary processes can be separated.

Conclusion
The goal of this manuscript is the analysis of the invasion of a mutant strain when introduced into a wild-type population. Assuming a constant population size or deterministically varying population sizes this quantity has already been studied extensively. Here, we extend the analysis to systems including stochastic demographic fluctuations. Dealing with a competitive Lotka-Volterra model we are able to approximate the fixation probability under the weak selection assumption, i.e. the interaction rates between individuals just differ slightly. This mechanistic implementation of the model allows for a straightforward interpretation of the selection effects in terms of competition rates, much in the spirit of Doebeli et al. (2017).
In order to obtain an expression for the fixation probability, we approximate the model by stochastic diffusions and apply tools from stochastic differential equation theory. We observe that the evolutionary success of a mutant mainly depends on its interspecific competition rate b. This is due to the resident type being more frequent initially yielding a higher probability for a mutant individual to interact with a resident type. This implies that a larger payoff for the mutant interacting with the wild-type ensures an enhanced fixation probability. This can be seen explicitly by the factor (1 − d b ) occurring in Eq. (9) which is the only term in the formula that can switch signs given an initially rare mutant, i.e. p < p * .
Still, the values a and c play a role in the overall evolutionary picture. While lowering c increases p * and thus the region where the invading type has a selective advantage, the parameter a has an impact on the overall population size after fixation of the mutant strain. This might end in a decrease of the total number of individuals if a < d, even though the mutant has a selective advantage over the wild-type due to b > d. Furthermore, we studied the fixation probability of exactly one mutant in the initial population. Non-surprisingly and as already observed in systems with deterministic population growth/decrease, cf. Kimura and Ohta (1974), we see that the fixation probability monotonically decreases for increasing initial population sizes. Additionally, we find that in our system due to the varying population size the famous 1 3 -rule for fixed population sizes, see Nowak et al. (2004), does not hold anymore. However, we can relate the deterministic internal equilibrium to the intersection of the neutral fixation probability and its counterpart including selection, i.e. ϕ neu ( p * , z) = ϕ( p * , z) = p * . This has already been observed in Constable and McKane (2015) and Lambert (2006) but not been stated explicitly We also explored the parameter range under which our modeling assumptions are valid. As it turns out the regime is different when compared to results from population genetics. There, usually Ms 1 is assumed while in our framework the approximation holds for Ms < 1. This is due to the mechanistic implementation of the population dynamics while in population genetics the former condition is usually assumed in order to derive a reasonable stochastic diffusion limit.
The evolutionary result of populations under stochastically fluctuating population sizes has been studied in various scenarios over the last few years (Melbinger et al. 2010;Chotibut and Nelson 2015;Constable et al. 2016;Chotibut and Nelson 2017). The stochasticity of the system as opposed to a deterministic modeling approach allows for different asymptotic behaviors and especially can reverse the deterministic behavior. This triggers the question for calculating fixation probabilities. We added new insight on the impact of the competition parameters on the fixation probability. Furthermore, we showcase a method from stochastic diffusion theory and developed in Lambert (2006) for approximating this quantity at least in the weak selection limit. Even though it is limited to the study of two interacting species, it is adaptable to many other models (and not only Lotka-Volterra-type systems) which include stochastic variation on the population size level.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Derivation of the diffusion approximation
We derive the Fokker-Planck equation corresponding to our model. The birth-and death-processes are given by with β X , γ X and β Y , γ Y being the birth and death rates. The competition processes read with M scaling the population size in stationarity. We follow the derivation of the Fokker-Planck equation as done in Huang et al. (2015) for the very same model (see also van Kampen 1997;Gardiner 2004). We set as the transition rates of the system and derive the master equation. That is, for p t (X t , Y t ) being the probability for the system to be in state (X, Y ) at time t we have Rescaling the parameters, i.e. setting x = X M and y = Y M , we get In the following we neglect the time subscript of the variables x and y. Then, doing a Taylor expansion of the function f and the transition rates T around (x t , y t ) up to the second order yields Inserting the terms for the transition rates T +/− x/y gives the Fokker-Planck equation which corresponds to the diffusion equation given in Eq. (4), see for example (van Kampen 1997;Gardiner 2004). ries resulting in the extinction of type Y individuals under the condition that type X individuals are still present.
On the other hand we have where the operator G is given in Eq. (7). Hence, we need to solve Gϕ = 0 with boundary conditions ϕ(0, y) = 0 and ϕ(x, 0) = 1 for x, y > 0 which follow immediately from the model dynamics which then gives the partial differential equation with boundary conditions as stated in Eq. (6).
Hence, the generator changes to + Next, applying the weak selection limit, i.e. using conditions (i)-(iii) we obtain is positive. Since the following section is very technical the reader who is satisfied with a numerical argument should skip this section and continue reading at "Appendix F".
(iv) w decreases in the neighborhood of 0+, we get = 1 z 4 h(z) − Next, we look for a combination of these derivatives such that these give the inhomop t+δt = p t + with Lastly, we need to identify the deterministic part of the system with weak selection. Therefore, we need to compute Eq. 19 in the reference. In order to do so, we need to transform their variables into our parameter framework. Transforming the birthand death rates is straightforward while for competition rates we need to do some arithmetics. The competition coefficients in Constable and McKane (2015), denoted by c i j , can be written as follows, exemplarily done for c 12 : 1 b = c 12 = c 0 (1 + εγ 12 ) ⇔ γ 12 = 1 ε 1 bc 0 − 1 .
As the reference case we set c 0 = 1 such that we can now transform between the two systems. Their Eq. 19 now transforms to In summation we now have expressions for B(q) and A(q) which not yield the stochastic differential equation which is also stated in Eq. (20) in the main text.