Effects of niche overlap on coexistence, fixation and invasion in a population of two interacting species

Synergistic and antagonistic interactions in multi-species populations—such as resource sharing and competition—result in remarkably diverse behaviours in populations of interacting cells, such as in soil or human microbiomes, or clonal competition in cancer. The degree of inter- and intra-specific interaction can often be quantified through the notion of an ecological ‘niche’. Typically, weakly interacting species that occupy largely distinct niches result in stable mixed populations, while strong interactions and competition for the same niche result in rapid extinctions of some species and fixations of others. We investigate the transition of a deterministically stable mixed population to a stochasticity-induced fixation as a function of the niche overlap between the two species. We also investigate the effect of the niche overlap on the population stability with respect to external invasions. Our results have important implications for a number of experimental systems.

One common theory, known as the Gause's rule or the competitive exclusion principle, postulates that due to abiotic constraints, resource usage, inter-species interactions, and other factors, ecosystems can be divided into ecological niches, with each niche supporting only one species in steady state; that species is termed to have fixated within the niche [18][19][20][21].Commonly, the number of ecological niches can be related to the number of limiting factors that affect growth and death rates, such as metabolic resources or secreted molecular signals like growth factors or toxins, or other regulatory molecules [22][23][24][25].The strength of interspecies can be commonly characterized by a quantity known as niche overlap [30,33,[77][78][79].Observed biodiversity can also arise from the turnover of transient mutants or immigrants that appear and go extinct in the population [26][27][28].However, the exact definition of an ecological niche is still the subject of debate and varies among different works [4,[29][30][31][32][33][34][35].The maintenance of the biodiversity of species that occupy similar niches is still not fully understood [25,36,37].
Deterministically, ecological dynamics of mixed populations has been commonly described as a dynamical system of equations that governs the dynamics of the numbers of individuals of each species and the concentrations of the limiting factors [22][23][24].Steady state co-existence of multiple species commonly corresponds to a stable fixed point in such a dynamical system, and the number of stably co-existing species is typically constrained by the number of limiting factors [22][23][24].In some cases, deterministic models allow co-existence of more species than the number of limiting factors, for instance when the attractor of the system dynamics is a limit cycle rather than a fixed point [24,38].Particularly pertinent for this paper is the case when the interactions of the limiting factors and the target species have a redundancy that results in the transformation of a stable fixed point into a marginally stable manifold of fixed points.This situation underlies many neutral models of population dynamics and evolution [19,20,26,[39][40][41][42].However, in this situation the stochastic fluctuations in the species numbers become important [24,25,[43][44][45][46].
Stochastic effects, arising either from the extrinsic fluctuations of the environment [47,48], or the intrinsic stochasticity of the individual birth and death events within the population [35,40,[49][50][51][52][53][54][55], modify the deterministic picture.The latter type of stochasticity, known as demographic noise, is the focus of this paper.Demographic noise causes fluctuations of the populations abundances around the deterministic steady state until a rare large fluctuation leads to an extinction of one of the species [20,40,54].In systems with a deterministically stable co-existence point, the mean time to extinction is typically exponential in the population size [51,[56][57][58][59] and is commonly considered to imply stable long term co-existence for typical biological examples with relatively large numbers of individuals [59,60].
By contrast, in systems with a stable neutral manifold that restore fluctuations out of the manifold but not along it, mean extinction times scale as a power law with the population size, indicating that the co-existence fails in such systems on biologically relevant timescales [39,40,51,61,62].This type of stochastic dynamics parallels the stochastic fixation in the Moran-Fisher-Wright model that describes strongly competing populations with fixed overall population size, and is the classical example of the class of models known as neutral models [9,14,20,34,39,[63][64][65].
A broad class of dynamical models of multi-species populations interacting through limiting factors can be mapped onto the class of models known as generalized Lotka-Volterra (LV) models, which allow one to conveniently distinguish between various interaction regimes, such as competition or mutualism, and which have served as paradigmatic models for the study of the behavior of interacting species [33, 35, 40, 43-46, 51-55, 66, 67].Remarkably, the stochastic dynamics of LV type models is still incompletely understood, and has recently received renewed attention motivated by problems in bacterial ecology and cancer progression [9,25,34,35,54,[68][69][70].A number of approximations have been used to understand the stochastic dynamics of the LV model, predominantly in the weak competition (independent) regime, or the opposite neutral limit [35,40,[51][52][53][54]71].However, the approximations tend to break down for the strong competition near the neutral limit, which is of central importance for understanding the deviations of ecological and physiological systems from neutrality [35,67,[72][73][74].
In this paper, we analyse a Lotka-Volterra model of two competing species using the master equation and first passage formalism that enables us to obtain numerically exact solution to arbitrary accuracy in all regimes, avoiding inaccuracies of various approximate descriptions of the stochastic dynamics of the system.The emphasis of this paper is on the effect of the ecological niche overlap on the transition from the deterministic co-existence to the stochastic fixation.In particular, we focus on the balance between the stochastic and deterministic effects in the strong competition regime near complete niche overlap, as expressed in the interplay between the exponential and algebraic terms in the extinction time as a function of the population size.Furthermore, we study the complementary question of the population stability with respect to the invasion of an immigrant or a mutant into a stable ecological niche of an already established species.This analysis informs many theories of transient diversity, such as Hubbell's neutral theory, and serves as a stepping stone toward multi-species models in non-neutral regime.
The paper is structured as follows.Section 2 discusses the definition of an ecological niche and briefly examines the regimes of deterministic stability of the system.In Section 3 we introduce the stochastic description of the LV model.In Section 3.1 we analyse the fixation/extinction times as a function of the niche overlap between the two species.We find that for any incomplete niche overlap the extinction time contains both terms that scale exponentially and algebraically with the system size, only reaching the known algebraic scaling at complete niche overlap; we calculate the exact functional form of these dependencies.Section 3.2 discusses the dynamic underpinning of the results of Section 3.1.In Section 3.3 we consider the invasion of an immigrant or a mutant into a stable ecological niche of an already established species, and calculate the invasion probabilities as a function of the system size and the niche overlap.The probability of success of an invasion attempt depends strongly on the competition strength, decreasing linearly with niche overlap, while only weakly dependent on the system size.All invasion attempts occur rapidly regardless of the outcome, or of the niche overlap between the two species or system size.We conclude with the Discussion section of these results in the context of previous works, and their implications for a number of experimental systems.
2 Deterministic Behavior of the Lotka-Volterra Model: Overview at the indicated values of the niche overlap a for a 12 = a 21 ≡ a.The fixed point is stable for a < 1.At a = 0 the two species evolve independently.As a increases, the deterministically stable fixed point moves toward the origin.At a = 1 the fixed point degenerates into a line of marginally stable fixed points, corresponding to the Moran model.The dashed lines illustrate the deterministic flow of the system: blue for a = 0.5, and orange for a = 1.2.The zoom inset illustrates the stochastic transitions between the discrete states of the system.Fixation occurs when the system reaches either of the axes.See text for details.
The model being employed in this paper is based on the generalized two dimensional Lotka-Volterra (LV) equations, that have served as a paradigmatic system to study the behavior of interacting populations [33, 35, 40, 43-46, 51-55, 66, 67].An example of the derivation of the model, from the underlying mechanistic model describing the exchange of control factors, is shown in the Supplementary Information.The dynamical equations of the numbers of individuals of species 1 (denoted as x 1 ) and species 2 (denoted as x 2 ) are [54,75,76]: The turnover rates r i set the timescales of the birth and death for each species, and K i are known as the carrying capacities.The interaction parameters a ij provide a mathematical representation of the intuitive notion of the niche overlap between the species, commonly used to characterize the competition of multispecies systems.These niche overlap parameters, which represent the competition between species arising from their shared usage of resources and other factors, dictate the long-term co-existence, viability and the diversity of the population, as well as the apportionment of the different species among different niches [30,33,[77][78][79].The niche overlap parameters can be directly derived from an underlying model describing the production-consumption dynamics of the limiting factors that control the birth and death rates in the system (see Supplementary Information).When a ij = 0, species j does not affect the species i, and they occupy separate ecological niches.At the other limit, a ij = 1, the species j compete just as strongly with species i as species i does within itself, and both species occupy same niche.The deterministic behavior of the generalized Lotka-Volterra model is well known, and can result in various stability regimes, including co-existence, competitive exclusion, mutualism, and bi-stability [54,75,76], as reviewed in the Supplementary Information.Here, we briefly review the behavior of the deterministic Equations (2.1), which have four fixed points [54,75,76]: The origin O is the fixed point corresponding to both species being extinct, and is unstable with positive eigenvalues equal to r 1 and r 2 along the corresponding on-axis eigendirections.The single species fixed points A and B are stable on-axis (with eigenvalues −r 1 and −r 2 , respectively), but are unstable with respect to invasion if point C is stable, reflected in the positive second eigenvalue equal to r 2 (1 − a 21 K 1 /K 2 ) and r 1 (1 − a 12 K 2 /K 1 ), respectively.Fixed point C corresponds to the co-existence of the two species and is stable when both a ij < 1 for K 1 = K 2 , also known as the weak competition regime ( [54,75,76] and Supplementary Information).
In the special case a 12 = a 21 = 1 the stable co-existence point degenerates into a neutral line of stable points, defined by x 2 = K − x 1 , as shown in Figure 1.Each point on the line is stable with respect to perturbations off line, but any perturbations along the line are not restored to their unperturbed position [23,80].This line correspond to the singular case of complete niche overlap where the two species are functionally identical with respect to their interactions with limiting factors like resources, space, toxins, predators, etc. (see also Supplementary Information).The stochastic dynamics along this line correspond to the classical Moran model [40,[53][54][55] as discussed below, and in the following we refer to this line as the Moran line.
Figure 1 shows the phase portrait of the system, in the symmetric case of K 1 = K 2 ≡ K, r 1 = r 2 ≡ r, and a 12 = a 21 ≡ a, where neither of the species has an explicit fitness advantage.This equality of the two species, also known as neutrality, serves as a null model against which systems with explicit fitness differences can be compared.
Stochastic effects, described in the next section, modify the deterministic stability picture and can lead to extinction or fixation of species even in deterministically stable case.In this paper, we focus on the effects of stochasticity in the deterministically stable weak competition regime (0 ≤ a ≤ 1), finding the scaling of the mean times to fixation and invasion as a function of the niche overlap a and comparing these to the known limits.The results in the asymmetric case are qualitatively similar and are relegated to the Supplementary Information.

Effects of Stochasticity
Stochasticity naturally arises in the dynamics of the system from the randomness in the birth and death times of the individuals -commonly known as the demographic noise [49,[81][82][83].Competitive interactions between the species can affect either the birth rates (such as competition for nutrients) or the death rates (such as toxins or metabolic waste), and in general may result in different stochastic descriptions [84,85].In this paper, we follow others [40,52,53] in considering the case where the inter-species competition affects the death rates so that the per capita birth and death rates b i and d i of species i are: The per capita death rate of a member of species i is increased by the presence of other members of the same species, and to a lesser effect (in proportion to their niche overlap) also by members of the other species.In the deterministic limit of negligible fluctuations the model recovers the mean field competitive Lotka-Volterra Equations (2.1) [40].The effects of an intrinsic death rate and competition modifying birth rates will be studied elsewhere.
The system is characterized by the vector of probabilities P (s, t|s 0 ) to be in a state s = {x 1 , x 2 } at time t, given the initial conditions s 0 = (x 0 1 , x 0 2 ): P (t) ≡ . . ., P (s, t|s 0 ), . . .[86].The forward master equation describing the time evolution of this probability distribution is [81] where M is the (semi-infinite) transition matrix, with more details given in the following paragraph.
Because the approximate analytical and semi-analytical solutions of the master Equation (3.2) often do not provide correct scaling in all regimes ( [85,[87][88][89]; see also the Supplementary Information), we analyse the master equation numerically in order to recover both the exponential and polynomial aspects of the mean time to fixation.To enable numerical manipulations, we introduce a reflecting boundary condition at a cutoff population size C K > K for both species to make the transition matrix finite [86,90,91] and enumerate the states of the system with a single index [86] via the mapping of the two species populations (x 1 , x 2 ) to state s as where s serves as the index for our concatenated probability vector.
In this representation, the non-zero elements of the sparse matrix M are Ms,s = −b 1 (s)−b 2 (s)−d 1 (s)−d 2 (s) along the diagonal, giving the rate of transition out of state (x 1 , x 2 ) to the adjacent states (x 1 + 1, x 2 ), (x 1 , x 2 + 1), (x 1 − 1, x 2 ), or (x 1 , x 2 − 1) respectively, corresponding to a species 1 or 2 birth or death event.Off from the diagonal, the elements Ms,s+1 = d 2 (s + 1) give the transition rate from the state of populations (x 1 , x 2 ) of the two species to the state (x 1 , x 2 − 1) organisms, corresponding to the death of a member of species 2. Off-diagonal elements Ms+1,s = b 2 (s) are the transition from state (x 1 , x 2 ) to state (x 1 , x 2 + 1), similarly corresponding to the birth of an organism from species 2. The off-diagonal elements at ±C K are the remaining two transitions: the death of species 1 member is given by Ms,s+CK = d 1 (s + C K ), and its birth is Ms+CK,s = b 1 (s).Some diagonal elements are modified to ensure the reflecting boundary at x i = C K .We have found that the choice C K = 5K is more than sufficient to calculate the mean fixation times to at least three significant digits of accuracy.

Fixation time as a function of the niche overlap
In this section we calculate the first passage times to the extinction of one of the species and the corresponding fixation of the other, induced by demographic fluctuations, starting from an initial condition of the deterministically stable co-existence point.The stochastic system tends to fluctuate near the deterministic fixed point (or Moran line), but rare fluctuations can take the system far from equilibrium.Should the fluctuations bring the system to either of the axes, then one species has gone extinct.In this model each axis acts as an absorbing manifold, from which the system cannot recover.Biologically, once a species has gone extinct from the system it can no longer reproduce.The system then proceeds to evolve as a one species logistic model with demographic noise, a process that is well-studied in the literature [56,59,92].
The time the system takes to fixation is a random variable, the distribution of which can be calculated in a number of ways.The master Equation (3.2) has a formal solution obtained by the exponentiation of the matrix: P (t) = e M t P (0).Then ( P ) s=(x1,0) and ( P ) s=(0,x2) give the cumulative distribution of the fixation of species 1 and 2 respectively.However, direct matrix exponentiation, as well as direct sampling of the master equation using the Gillespie algorithm [93,94], are impractical since the fixation time grows exponentially with the system size; nevertheless, we used Gillespie tau-leaping simulations to verify our results up to moderate system size (see Supplementary Information).More amenably, the moments of the first passage time distribution can be calculated directly without explicitly solving the master equation [95].In an ensemble of trajectories, the fraction of trajectories in state s at time t is given by P (s, t|s 0 ).Thus, averaging over the ensemble, the mean residence time in state s during the system evolution is the integral of this quantity over all time [95]: For each trajectory, the time to reach fixation is simply the sum of the times it spends in each state on the way.Thus the mean time to fixation starting from a state s 0 is the sum of the residency times [96]: This expression can be also derived using the backward master equation formalism [96].The matrix inversion was performed using LU decomposition algorithm implemented with the C++ library Eigen [97], which has algorithimic complexity of the calculation scaling algebraically with K. Increasing the cutoff C K enables calculation of the mean fixation times to an arbitrary accuracy.
The left panel of Figure 2 shows the calculated fixation times with the initial condition at the deterministically stable co-existence fixed point as a function of the carrying capacity K for different values of the niche overlap a.While the transition has not previously been characterized, the independent and neutral limits are known, providing comparisons at the a = 0 and a = 1 extremes.
In the limit of non-interacting species (a = 0), each species evolves according to an independent stochastic logistic model, and the probability distribution of the fixation times is a convolution of the extinction time distributions of a single species, which are dominated by a single exponential tail [56,59,92] (see also the Supplementary Information).Mean extinction time of a single species can be calculated exactly and is well known, and asymptotically for K 1 it varies as 1 K e K [98,99].Denoting the probability distribution of the extinction times for either of the independent species as p(t) = αe −αt and its cumulative as f (t) = t s=0 p(s)ds, the probability that either of the species goes extinct in the time interval , where is the mean time of the single species extinction.This analytical limit τ 1 2K e K is shown in Figure 2 in a black dashed line and agrees well with the numerical results of Equation (3.5).From the biological perspective, for sufficiently large K, the exponential dependence of the fixation time on K implies that the fluctuations do not destroy the stable co-existence of the two species.
In the opposite limit of complete niche overlap, a = 1, any fluctuations along the line of neutrally stable points are not restored, and the system performs diffusion-like motion that closely parallels the random walk of the classic Moran model [35,40,46,51,53,54,71] (see also the Supplementary Information).The Moran model shows a relatively fast fixation time scaling algebraically with K [39,40], τ ln(2)K 2 ∆t.This known result can be obtained from the backward Fokker-Planck equation ∂ 2 τ ∂x 2 = −∆t K 2 x(K−x) of the Moran model, which directly gives the above scaling for equal population initial conditions [39,40].The fixation time predicted by the Moran model is shown in Figure 2 as a red dotted line and shows excellent agreement with our exact result.Note that the average time step ∆t in the corresponding Moran model is ∆t ≈ 1/K because the mean transition time in the stochastic LV model is proportional to 1/(rK) close to the Moran line [54]; see the Supplementary Information for more details.The relatively short fixation time in the complete niche overlap regime implies that the population can reach fixation on biologically realistic timescales.
The exponential scaling of the fixation time with K persists for incomplete niche overlap described by the intermediate values of 0 < a < 1.However, both the exponential and the algebraic prefactor depend on the niche overlap a.The exponential scaling is expected for systems with a deterministically stable fixed point [35,52,59,88,89], as indicated in [40,51,54] using Fokker-Planck approximation and in [52] using the WKB approximation.However, the Fokker-Planck and WKB approximations, while providing the qualitatively correct dominant scaling, do not correctly calculate the scaling of the polynomial prefactor and the numerical value of the exponent simultaneously [59,71,85].For large population sizes and timescales, effective species co-existence will be typically observed experimentally whenever the fixation time has a non-zero exponential component.
To quantitatively investigate the transition from the exponentially stable fixation times to the algebraic scaling in the complete niche overlap regime, we use the ansatz τ (a, K) = e h(a) K g(a) e f (a)K . (3.6) In the Moran limit, a = 1, we expect f (1) = 0, g(1) = 1.In the independent species limit with zero niche overlap, a = 0, we expect f (0) = 1 and g(0) = −1.The right panel of Figure 2 shows the ansatz functions f (a), g(a), and h(a), obtained by numerical fit to the fixation times as a function of K shown in the left panel.
The numerical results agree well with the expected approximate analytical results for a = 0 and a = 1 with small discrepancies attributable to the approximate nature of the limiting values.Notably, f (a), which quantifies the exponential dependence of the fixation time on the niche overlap a, smoothly decays to zero at a = 1: only when two species have complete niche overlap (a = 1) does one expect rapid fixation dominated by the algebraic dependence on K.In all other cases the mean time until fixation is exponentially long in the system size [59,92].Even two species that occupy almost the same niche (a 1) effectively co-exist for K 1, with small fluctuations around the deterministically stable fixed point.The dependence of the exponential function e f (a)K on a and K can be understood in the spirit of Kramers' theory, as discussed in the next section.
Interestingly, for small K the algebraic prefactors may dominate the exponential scaling such that in a region of a and K space, the fixation time from a deterministically stable fixed point (with a < 1) can be shorter than that of a Moran case with the same K. See the Supplementary Information for more details.Also explored in the Supplementary Information is a breaking of the parameter symmetry, which gives similar behavior to the above.In particular, the exponential dependence of the escape time from the fixed point on the system size also persists in the non-neutral case, although the dependence can be much weaker.

Route to fixation and the origin of the exponential scaling
To gain deeper insight into the fixation dynamics and to understand the origins and the limits of the exponential scaling with K, in this section we calculate the residency times in each state during the fixation process, given by Equation (3.4).The results are shown as a contour plot in Figure 3, for two different niche overlaps, one at and the other far from the Moran limit.The set of states lying along the steepest descent lines of the contour plot, shown as the black line, can be thought of as a "typical" trajectory [52,71,100].However, even for two species close to complete niche overlap the system trajectory visits many states far from this line.This departure is even greater for weakly competing species, where the system covers large areas around the fixed point before the rare fluctuation that leads to fixation occurs [50].These deviations Figure 3: The system samples multiple trajectories on its way to fixation.Contour plot shows the average residency times at different states of the system, with pink indicating longer residence time, deep green indicating rarely visited states.The colored line is a sample trajectory the system undergoes before fixation; color coding corresponds to the elapsed time with orange at early times, purple at the intermediate times and red at late stages of the trajectory.The red dot shows the deterministic co-existence point.See text for more details.Left: Complete niche overlap limit, a = 1, for K = 64.Right: Independent limit with a = 0 and K = 32.from a "typical" trajectory are related to the inaccuracy of the WKB approximation in calculating the scaling of the pre-exponential factor [50,89,99]; see also the Supplementary Information.This occupancy landscape can be qualitatively thought of as an effective Lyapunov function/effective potential of the system dynamics [101], although the LV system does not possess a true Lyapunov function.
A similar issue also arises in the Fokker-Planck approximation [54,101].Nevertheless, the pseudo-potential landscape provides an intuitive underpinning for the general exponential scaling in the incomplete niche overlap regime: the fixation process can be thought of as the Kramers-type escape from a pseudo-potential well [102].The Kramers' escape time is dominated by the exponential term τ exp(∆U ), where ∆U is the depth of the effective potential well, corresponding to the dominant scaling τ exp(f (a)K) of this paper.The depth of the pseudo-potential well can be estimated in the Gaussian approximation to the Fokker-Planck equation [81,103,104].For y i ≡ x i /K, the Fokker-Planck equation is , similar terms for F 2 and D 22 , and D 12 = D 21 = 0. Linearized about the deterministic fixed point y * the Fokker-Planck equation becomes The steady state solution to this approximation is a Gaussian distribution centered on the deterministic fixed point.It gives an effective well depth of ∆U = (1−a) 2(1+a) K, providing a qualitative basis for the numerical results in the right panel of Figure 2. As a increases and the species interact more strongly, the potential well becomes less steep, resulting in weaker exponential scaling.In the complete niche overlap limit, the pseudo-potential develops a soft direction along the Moran line that enables relatively fast escape towards fixation.The change in the topology of the pseudo-potential is also reflected in the Pearson correlation coefficient between the two species: The covariance ranges from zero in the independent limit (as expected) to anti-correlated in the Moran limit, reflecting the fact that the trajectories typically diffuse along the trough around the Moran line.

Invasion of a mutant/immigrant into a deterministically stable population
Transient co-existence during the fixation/extinction process of immigrants/mutants has also been proposed as a mechanism for observed biodiversity in a number of contexts [26,31,67,70,72,105,106].The extent of this biodiversity is constrained by the interplay between the residence times of these invaders and the rate at which they appear in a settled population.In the previous sections we calculated the fixation times in the two species system starting from the deterministically stable fixed point.In this section we investigate the complementary problem of robustness of a stable population of one species with respect to an invasion of another species, arising either through mutation or immigration, and investigate the effect of niche overlap and system size on the probability and mean times of successful and failed invasions.
To this end, we study the case where the system starts with K − 1 individuals of the established species and 1 invader.Each species' dynamics is described by the birth and death rates defined by Equations (3.1) above.We consider the invasion successful if the invader grows to be half of the total population without dying out first.Note that this is different from the bulk of invasion literature [40, 46, 51-55, 61, 67, 90, 105, 107-110] (but see [111]), which typically considers invasion successful when a mutant fixates in the system.However, this is only sensible in the Moran limit, since otherwise the system tends to the co-existence fixed point, after which there is an equal chance of fixation and extinction, with the MTE calculated above; hence our definition of invasion.We denote the probability of a invader success as P, the mean time to a successful invasion as τ s , and the mean time of a failed invasion attempt, where the invader dies out before establishing itself in the population, as τ f .More generally, invasion probability and the successful and failed times starting from an arbitrary state s 0 are denoted as P s 0 , τ s 0 s and τ s 0 f , respectively.Similar to Equation (3.5) above, the invasion probability can be obtained from [96,103] and the times from where α s is the transition rate from a state s directly to extinction or invasion of the invader and Φ s 0 = τ s 0 P s 0 is a product of the invasion or extinction time and probability.Similar equations describe τ f [96,103].
Figure 4 shows the calculated invasion probabilities as a function of the carrying capacity K and of the niche overlap a between the invader and the established species.In the complete niche overlap limit, a = 1, the dependence of the invasion probability on the carrying capacity K closely follows the results of the classical Moran model, P s 0 = 2/K [39,90], shown in the blue dotted line in the left panel, and tends to zero as K increases.In the other limit, a = 0, the problem is well approximated by the one-species stochastic logistic model starting with one individual and evolving to either 0 or K individuals; the exact result in this limit is shown in black dotted line, referred to as the independent limit [103].In the independent limit, a = 0, the invasion probability asymptotically approaches 1 for large K, reflecting the fact that the system is deterministically drawn towards the deterministic stable fixed point with equal numbers of both species.Interestingly, the invasion probability is a non-monotonic function of K and exhibits a minimum at an intermediate/low carrying capacity, which might be relevant for some biological systems, such as in early cancer development [13] or plasmid exchange in bacteria [17].
For the intermediate values of the niche overlap, 0 < a < 1, the invasion probability is a monotonically decreasing function of a, as shown in the right panel of Figure 4.For large K, the outcome of the invasion is typically determined after only a few steps: since the system is drawn deterministically to the mixed fixed point, the invasion is almost certain once the invader has reproduced several times.At early times, the invader birth and death rates (3.1) are roughly constant, and the invasion failure can be approximated by the extinction probability of a birth-death process with constant death d mut and birth b mut rates.The invasion probability is then P = 1 − d mut /b mut ≈ 1 − a.This heuristic estimate is in excellent agreement with the numerical predictions, shown in the right panel of Figure 4 as a purple dashed and the blue lines respectively.Perhaps unexpectedly, at large K the dependence of the invasion probability on a is much stronger than its dependence on K, such that increasing the carrying capacity has little effect of the eventual outcome of the bmut+dmut is an appropriate small carrying capacity limit.Successive lines are at larger system size, and approach the solid magenta line of system, but a small change in the niche overlap between the invader and the established species leads to a large change in the probability of invasion success.For small K either invasion or extinction typically occurs after only a small number of steps.The invasion probability in this limit is dominated by the probability that the lone mutant reproduces before it dies, namely bmut bmut+dmut = K K(1+a)+1−a , as shown in black dotted line in the right panel of Figure 4.
The top left panel of Figure 5 shows the dependence of the mean time to successful invasion, τ s , on K and a. Increasing K can have potentially contradictory effects on the invasion time, as it increases the number of births before a successful invasion on the one hand, while increasing the steepness of the potential landscape and therefore the bias towards invasion on the other.Nevertheless, the invasion time is a monotonically increasing function of K for all values of a.In the complete niche overlap limit a = 1 the invasion time scales linearly with the carrying capacity K, as expected from the predictions of the Moran model, with ∆t 1/K, as explained above.The quantitative discrepancy arises from the breakdown of the ∆t 1/K approximation off of the Moran line.In the opposite limit of non-interacting species, a = 0, the invading mutant follows the dynamics of a single logistic system with the carrying capacity K, resulting in the invasion time that grows approximately logarithmically with the system size, as shown in the left panels of Figure 5 as a black line (see also Supplementary Information).For all values 0 ≤ a < 1 the invasion time scales sub-linearly with the carrying capacity, indicating that successful invasions occur relatively quickly, even when close to complete niche overlap, where the invading mutant strongly competes against the stable species.By contrast, the failed invasion time, shown in the bottom left panel of Figure 5, is non-monotonic in K.The analytical approximations of the Moran model and the of two independent 1D stochastic logistic systems recover the qualitative dependence of the failed invasion time on K at high and low niche overlap, respectively.All failed invasion times are fast, with the greatest scaling being that of the Moran limit.For a < 1 these failed invasion attempts approach a constant timescale at large K.
The dependence of the time of an attempted invasion (both for successful and failed ones) on the niche overlap a is different for small and large K, as shown in the right panels of Figure 5.For small K both τ s and τ f are monotonically decreasing functions of a, with the Moran limit having the shortest conditional times.In this regime, the extinction or fixation already occurs after just a few steps, and its timescale is determined by the slowest steps, namely the mutant birth and death.Thus τ ≈ Figure 5: Mean time of a successful or failed invasion attempt.Upper Left: Dotted lines connect the numerical results of invasion times conditioned on success, from a = 0 at the bottom being mostly fastest to a = 1 being slowest.The solid green line shows for comparison the predictions of the Moran model in the complete niche overlap limit, a = 1; see text.The solid purple line correspond to the solution of an independent stochastic logistic species, a = 0, and overestimates the time at small K but fares better as K increases.Upper Right: The red line shows the results of successful invasion time for carrying capacity K = 4, and successive lines are at larger system size, up to K = 256.The cyan line is 1/(b mut + d mut ) and matches with small K. Lower Panels: Same as upper panels, but for the mean time conditioned on a failed invasion attempt.
as shown in the figures as the dashed blue line.By contrast, at large K, the invasion time is a non-monotonic function of the niche overlap, increasing at small a and decreasing at large a.This behavior stems from the conflicting effect of the increase in niche overlap: on the one hand, increasing a brings the fixed point closer to the initial condition of one invader, suggesting a shorter timescale; on the other hand, it also makes the two species more similar, increasing the competition that hinders the invasion.

Discussion
Maintenance of species biodiversity in many biological communities is still incompletely understood.The classical idea of competitive exclusion postulates that ultimately only one species should exist in an ecological niche, excluding all others.Although the notion of an ecological niche has eluded precise definition, it is commonly related to the limiting factors that constrain or affect the population growth and death.In the simplest case, one factor corresponds to one niche, which supports one species, although a combination of factors may also define a niche, as discussed above.The competitive exclusion picture has encountered long-standing challenges as exemplified by the classical "paradox of the plankton" [4,31] in which many species of plankton seem to co-habit the same niche; in many other ecosystems the biodiversity is also higher than appears to be possible from the apparent number of niches [26,31,36,112,113].
Competitive exclusion-like phenomena appear in a number of popular mathematical models, for instance in the competition regime of the deterministic Lotka-Volterra model, whose extensive use as a toy model enables a mathematical definition of the niche overlap between competing populations [30,[77][78][79].Another classical paradigm of fixation within an ecological niche is the Moran model (and the closely related Fisher-Wright, Kimura, and Hubbell models) that underlies a number of modern neutral theories of biodiversity [19,20,26,[39][40][41][42].Unlike the deterministic models, in the Moran model fixation does not rely on deterministic competition for space and limiting factors but arises from the stochastic demographic noise.Recently, the connection between deterministic models of the LV type and stochastic models of the Moran type has accrued renewed interest because of new focus on the stochastic dynamics of the microbiome, immune system, and cancer progression [13,25,40,46,53,54,70,89].
Much of the recent studies of these systems employed various approximations, such as the Fokker-Planck approximation [35,40,51,53,54], WKB approximation [52,71] or game theoretic [46] approach.The results of these approximations typically differ from the exact solution of the master equation, used in this paper, especially for small population sizes [59,71,85,88,89].In this paper, we have interrogated stochastic dynamics of a system of two competing species using a numerically arbitrarily accurate method based on the first passage formalism in the master equation description.The algorithmic complexity of this method scales algebraically with the population size rather than with the exponential scaling of the fixation time, (as is the case with the Gillespie algorithm [93]).Our approach captures both the exponential terms and the algebraic prefactors in the fixation/extinction times for all population sizes and enables us to rigorously investigate the effect of niche overlap between the species on the transition between the known limits of the effectively stable and the stochastically unstable regimes of population diversity.
Stochastic fluctuations allow the system to escape from the deterministic co-existence fixed point towards fixation.If the escape time is exponential in the (typically large) system size, in practice it implies effective co-existence of the two species around their deterministic co-existence point.If the escape time is algebraic in K, as in the degenerate niche overlap case (closely related to the classical Moran model), the system may fixate on biological timescales [39,105].Nevertheless, for biological systems with small characteristic population sizes, exponential scaling may not dominate the fixation time, and the algebraic prefactors dictate the behavior of the system (see the Supplementary Information for more details).
The transition that we have characterized from the exponential scaling of the effective co-existence time to the rapid stochastic fixation in the Moran limit is governed by the niche overlap parameter.We find that the fixation time is exponential in the system size unless the two species occupy exactly the same niche.The numerical factor in the exponential is highly sensitive to the value of the niche overlap, and smoothly decays to zero in the complete niche overlap case.The implication is that two species will effectively co-exist unless they have exactly the same niches.
Our fixation time results can be understood by noticing that the escape from a deterministically stable co-existence fixed point can be likened to Kramers' escape from a pseudo-potential well [51,59,92,114], where the mean transition time grows exponentially with well depth [59].Approximating the steady state probability with a Gaussian distribution shows that this well depth is proportional to K(1−a) and disappears at a = 1.With complete niche overlap the system develops a "soft" marginally stable direction along the Moran line that enables algebraically fast escape towards fixation [51,54,90].Similar to the exponential term, the exponent of the algebraic prefactor is also a function of the niche overlap, and smoothly varies from −1 in the independent regime of non-overlapping niches to +1 in the Moran limit.
The fixation times of two co-existing species, discussed above, determines the timescales over which the stability of the mixed populations can be destroyed by stochastic fluctuations.Similarly, the times of successful and failed invasions into a stable population set the timescales of the expected transient co-existence in the case of an influx of invaders, arising from mutation, speciation, or immigration.The probability of invasion strongly depends on the niche overlap, and for large K decreases monotonically as 1 − a with the increase in niche overlap.By contrast, the effect of the carrying capacity on the invasion probability is negligible in the large K limit.
The mean times of both successful and failed invasions are relatively fast in all regimes, and scale linearly or-sublinearly with the system size K.Even for high niche overlap, where the invasion probability is low due to the strong competition between the species, those invaders who do succeed invade in a reasonably fast time.In this high competition regime, the times of the failed invasions are particularly salient because they set the timescales for transient species diversity.If the influx of invaders is slower than the mean time of their failed invasion attempts, most of the time the system will contain only one settled species, with rare "blips" corresponding to the appearance and quick extinction of the invader.On the other hand, if individual invaders arrive faster than the typical times of extinction of the previous invasion attempt, the new system will exhibit transient co-existence between the settled species and multiple invader strains, determined by the balance of the mean failure time and the rate of invasion [26,31,77,106].
The weaker dependence of the invasion times on the population size and the niche overlap, as compared to the escape times of a stably co-existing system to fixation, implies that the transient co-existence is expected to be much less sensitive to the niche overlap and the population size than the steady state co-existence.Of note, both niche overlap and the population size can have contradictory effects on the invasion times (as discussed in section 3) resulting in a non-monotonic dependence of the times of both successful and failed invasions on these parameters.
Our results indicate that the niche overlap between two species reflecting the similarity (or the divergence) in how they interact with their shared environment, is of critical importance in determining stability of mixed populations.This has important implications for understanding the long term population diversity in many systems, such as human microbiota in health and disease [2,3,12], industrial microbiota used in fermented products [16], and evolutionary phylogeny inference algorithms [14,15].Our results serve as a neutral model base for problems such as maintenance of drug resistance plasmids in bacteria [17] or strain survival in cancer progression [13].The theoretical results can also be tested and extended based on experiments in more controlled environments, such as the gut microbiome of a c.elegans [70], or in microfluidic devices [115].
The results of the paper can also inform investigations of multi-species systems [24-27, 31-35, 66, 77, 90, 116, 117].It is helpful to focus on one of the species within a large community of different species, which experiences a distribution of niche overlaps with other species in the system.Since extinction times away from complete niche overlap scale exponentially in the system size, the extinction of the given species is likely to be dominated by the shortest timescale, imposed by whichever other species has the greatest niche overlap with it, excluding it most strongly.Therefore, species extinctions are expected to occur more frequently in communities where the probability of high niche overlap between the species is high -for instance for distributions with high mean overlap, or high variance in the niche overlap distribution.These qualitative predictions compare favorably with the work of Capitan et al [33], which shows that increasing the niche overlap in a multi-species system results on average in a smaller number of co-existing species, decreasing biodiversity.The probability of invasion into a multi-species system is more complicated than the co-existence considerations, because, unlike the fixation times, invasion times do not scale exponentially with the system size, so no single pairwise interaction will dominate the timescale.These directions are beyond of the scope of the present work and will be studied elsewhere.

Supplementary Information: Effects of niche overlap on co-existence, fixation and invasion in a population of two interacting species
Minimal model of two interacting species and the derivation of the Lotka-Volterra model As a minimal example, in this section we introduce a model of two interacting species whose growth is constrained by two secreted factors, inspired by the works of MacArthur and others [1][2][3][4].Each species x i has basal per capita birth rate β i , death rate µ i .Each species secretes soluble factors t j at rates g ji .Each factor t i is degraded at a rate λ i and affects the death rate of each bacterium linearly with the efficacy e ij .
Positive e ij may correspond to metabolic wastes, toxins or anti-proliferative signals [5][6][7][8][9][10], while negative e ij would describe growth factors or secondary metabolites [6,11,12].The model kinetics is encapsulated in the following equations for the turnover of the species numbers: and the equations for the production and the degradation of the secreted factors: It is useful to recast Equations ( 1), (2) defining vectors x = (x 1 , x 2 ) and t = (t 1 , t 2 ), so that where we have the matrices X = , and Ê = e 11 /r 1 e 12 /r 1 e 21 /r 2 e 22 /r 2 .
In many experimentally relevant systems, such as communities of microorganisms and cells, the timescale of production, diffusion, and degradation of secreted factors is on the order of minutes [13], whereas cell division and death occurs over hours [14,15], and the dynamics of the turnover of the secreted factors can be assumed to adiabatically reach a steady state t * given by t * = Ĝ • x [16][17][18].In this approximation the dynamical equations for the species number reduce to Written explicitly, this becomes the familiar generalized two-species competitive Lotka-Volterra system [2, [19][20][21][22][23][24][25]: where The number of deterministically viable species is typically constrained by the number of limiting factors [3] (see the following section).Namely, if both matrices Ê and Ĝ are non-singular and invertible, Equations (3) possess a mixed fixed point given by x * = (EG) −1 1.If this fixed point is stable and positive, it corresponds to the co-existence of the two species in two different (albeit potentially overlapping) ecological niches.
When the matrix ( Ê • Ĝ) is singular (a 12 a 21 = 1), the co-existence fixed point x * = (EG) −1 1 does not exist, and the Equations (3) are satisfied only if the population of one (or both) of the species is zero.Biologically, this condition corresponds to the complete niche overlap between two species, whereby only one species can survive in the niche.(Of note, exclusion of one species by the other can also occur in non-singular cases, as discussed in the next section.)Nevertheless, even in the complete niche overlap case, multiple species can deterministically coexist within one niche if the matrix ( Ê • Ĝ) possesses a further degeneracy, K 1 /K 2 = a 12 = 1/a 21 , corresponding to an additional symmetry in the interactions of the species with the constraining factors, as illustrated in the next paragraph.
These mathematical notions can be understood in a biologically illustrative example, when the matrix Ê is singular, so that det( Ê) = 0. Any singular 2 × 2 real matrix can be written in the general form , where α, β and γ are arbitrary real numbers [26].In this case Equation ( 1) becomes so that both secreted factors effectively act as one factor with concentration t ≡ t 1 + βt 2 .With γ = 1 this corresponds to the classic notion of two species and only one limiting factor.The two equations cannot be simultaneously satisfied and the only solution of Equations ( 6) is either x 1 = 0 or x 2 = 0 (or both).This is one example of competitive exclusion due to competition within a single niche.Finally, when γ = 1 (corresponding to a 12 = 1/a 21 = K 1 /K 2 ), both the species and the secreted factors are functionally identical, and the Equations ( 6) allow multiple solutions lying on the line in phase space defined by x * = Ĝ−1 t * and 1 = α (t * 1 + βt * 2 ) [21,27]; in this case many different mixtures of the two species can be deterministically stable, depending on the initial conditions.These derivations above provide a mathematical definition and a biological illustration of the niche overlap between two interacting species, and can be extended to a general case of N species interacting via M factors [27].

Stability analysis of the LV model
Different regions of the parameter space of the Lotka-Volterra model, shown in Figure 1, have different biological interpretations [19,[28][29][30][31]. Parasitism, or predation/antagonism, occurs when a 12 a 21 < 0, with one species gaining from a loss of the other.In the strong parasitism regime where the positive a ij is greater than one, the parasite/predator drives the prey to extinction deterministically, and the only stable point is the predator's fixed point (A or B).Conversely, weak parasitism allows co-existence of both species despite the detriment of one to the benefit of the other [19,30].
Both a ij < 0 corresponds to mutualistic/symbiotic interactions between the species [19,[28][29][30].Weak mutualism is mathematically similar to weak competition in that it results in stable co-existence.Strong mutualism with a 12 a 21 < −1 results in population explosion.Detailed study of this regime lies outside of the scope of the present work (but see [32]).
The quadrant with both a 12 > 0 and a 21 > 0 corresponds to the competition regime.At strong competition with either a 12 or a 21 greater than one, either the system is bistable and possesses two single-species stable fixed points A and B with separate basins of attraction (a 12 > 1 and a 21 > 1) or one of the species deterministically outcompetes the other.The complete niche overlap regime is given by the line a 12 a 21 = 1.These regimes correspond to the classical competitive exclusion theory, together with the strong parasitism  5) is stable in the green region and unstable in the blue region; in the white regions it is non-biological.Colored dots indicate the parameter range studied in the paper.The numbered regions correspond to different biological different regimes; see main text.Regions 4-6 correspond to competitive exclusion, with only single species fixed point A or B being stable (or both, in the bistable regime 5).In region 7 the populations experience unbounded growth.For the degenerate case a 12 = a 21 = 1, indicated by the red dot, the co-existence fixed point is replaced by a line of marginal stability, which we call the Moran line.
case.By contrast, weak competition where both 0 < a ij < 1 results in the stable co-existence at the mixed point C. In the special case a 12 = a 21 = 1 the stable fixed point degenerates into a neutral line of stable points, defined by x 2 = K − x 1 , as shown in Figure 1.

Long-term deterministic stability of interacting species
Quite generally, the dynamics of a system of N asexually reproducing species that interact with each other only through M limiting factors (such as food, soluble signaling and growth/death factors, toxins, metabolic waste) and experience no immigration, can be described by the following system of differential equations for the species x 1 , ..., x N and the limiting factor densities f 1 , ..., f M [3,27,33]: where f is the state of all factors that might affect the per capita birth rate β i f and the death rate µ i f of the species i.
The density of a factor j in the environment, f j , follows its own dynamical production-consumption equation ḟj = g j ( f , x) − λ j ( f , x)f j (8) where g j is a production-consumption rate that includes both the secretion and the consumption by the participating species as well any external sources of the factor f j , and λ j is its degradation rate.Alternatively, for some abiotic constrains such as physical space or amount of sunlight, the concentration of the factor f j can be set through a conservation equation of a form f j = c j ( f , x) [3,27].The fixed points of the N + M Equations ( 7) and (8) determine the steady state numbers of each of the N species and the corresponding concentrations of the M limiting factors.However, the structure of Equations (7) imposes additional constraints on the steady state solutions: at a fixed point β i f = µ i f for each of the N species, which determines the steady state concentrations of the M limiting factors f .However, if N > M , the system (7) of N equations is over-determined and typically does not have a consistent solution, unless the fixed point populations of N − M of the species are equal to zero [3,16,27,33,34].This reasoning provides a mathematical basis for the competitive exclusion principle, whereby the number of independent niches is determined by the number of limiting factors, and a system with M resources can sustain at most M species in steady state.
Nevertheless, as mentioned in the Introduction, the number of species at the steady state can exceed the number of limiting factors, when the N equations for the species are not independent and thus provide less than N constraints on the solutions.In this case, at steady state the populations of the non-independent species typically converge onto a marginally stable manifold on which each point is stable with respect to off-manifold perturbations but is neutral within the manifold [20,27,[35][36][37].

Mean fixation time in the classical Moran model
Here we re-derive the mean fixation time for the Moran model [38], to provide context for the results of the main text.The Moran model results come about as one limit of the stochastic Lotka-Volterra system, when niche overlap is complete.In the classical Moran model, at each time step, an individual is chosen at random to reproduce from the total population of size K.In order to keep the population constant, another one is chosen at random to die.The probabilities that the number of individuals of species 1 increases (b M (n)) or decreases (b M (n)) by one in one time step are [38]: where n is the number and f is the fraction of species 1 in the system (of total system size K).
The mean fixation time, τ (n), starting from some initial number n of species 1 is described by the following backward recursion equations [39]: where ∆t is the time step.Substituting the values of the 'birth' and 'death' probabilities of species 1 from Equation ( 9) we get At K 1, the Kramers-Moyal expansion in 1/K results in Integrating, using the boundary conditions τ (0) = τ (K) = 0, gives For the initial condition analogous to the co-existence point, n = K/2, this gives For direct comparison with the coupled logistic model we need to consider that in the Moran model one time step ∆t corresponds to one birth and one death event.Since the coupled logistic model spends most of its time near the WFM line x 1 + x 2 = K we assume that on average where b i and d i are the birth and death rates of the coupled logistic model.Since at the initial conditions the populations of each species are equal, and since b i (K/2, K/2) = d i (K/2, K/2) = K/2, we get ∆t ≈ 1/K.Therefore τ = ln(2) K, (12) as in the main text.The fixation time of the Moran model agrees well with the results of the coupled logistic model for complete niche overlap, as shown in Figure 2.

Exact and approximate mean extinction time for a single species stochastic logistic model
In this section, we re-derive the known results for the extinction time of a single species logistic model.Exact calculation.The mean extinction times for different initial states n 0 obey the usual backward recursion relation [39] τ where b(n) = r n and rate d(n) = r n n K are the birth and death rates of the process, respectively.The equation above can be rewritten as where and The logistic process becomes extinct in finite time, and the extinction time can be thus written as follows [39]: Evaluating this sum with b(n) = rn, d(n) = rn 2 /K [40] and the initial condition n 0 = K 1 gives the asymptotic limit to leading order [41].This result agrees with the more general approach of [42].
Fixation time of the coupled logistic model in the independent limit Here we calculate the mean fixation time in the independent limit of the coupled logistic model.The fixation occurs when either of the species goes extinct.Denoting the probability distribution of the extinction times for either of the independent species as p(t) and its cumulative as f (t) = t s=0 p(s)ds, the probability that either of the species goes extinct in the time interval [t, t + dt], is The mean time to fixation t is As shown in Figure 3, the probability distribution of the fixation times of a single species is dominated by the exponential tail.It can be approximated as with 1 α 1 K e K from the previous section.Finally, we obtain for the mean time to fixation of the two-species logistic model which leads to the equation τ  The bulk of the probability density is well described by an exponential distribution with the same mean, shown in the red dotted line.Data are generated using using the Gillespie algorithm for K = 16.For higher carrying capacities the assumption of exponentially distributed times becomes even more accurate.

Fokker-Planck approximation and the pseudo-potential
In this section we derive an analytical approximation for the results in Figure 2 in the main text, using Fokker-Planck approximation.The Fokker-Planck approximation to the coupled logistic system studied herein takes its traditional form [43]: where F is the drift (or force) vector and D is the diffusion matrix (in this case diagonal).Here, under symmetric conditions and non-dimensionalization by r, F i = ni K (K −n i −a ij n j ) and D ii = ni K (K +n i +a ij n j ).In general, Equation ( 23) cannot be reduced to diffusion in a potential U ( n) with an equilibrium distribution function P ( n) ∼ exp(U ( n)).The condition of zero flux at equilibrium, J i = F i P − 1/2 j ∂ j D ij P = 0, would require [39,43] ∂ However, for consistency it also requires [39,43].It is easy to show that this is not upheld for the two directions unless a 12 = a 21 = 0 and the system can be decomposed into two onedimensional logistic systems.Instead, we define the pseudo-potential as: where P ss (n 1 , n 2 ) is a quasi-stationary probability distribution function [44].We calculate P ss (n 1 , n 2 ) in the approximation to the Fokker-Planck Equation ( 23) linearized about the deterministic co-existence fixed point [39,43,45].The linearized equation is [39,43] where The quasi-equilibrium solution to Equation ( 25) is , a Gaussian centered on the co-existence point and with a variance given by the covariance matrix C = B • A −1 /2 in the symmetric case a 12 = a 21 = a, K 1 = K 2 = K [43].In this case the diagonal term of C is 1 1−a 2 K and gives the variance of a species about its mean value.The off-diagonal, which corresponds to the covariance between the two species, is − a 1−a 2 K. Thus the Pearson correlation coefficient between the two species is −a.
For the initial condition at the co-existence fixed point and assuming that the system escapes towards fixation once it reaches one of the axial fixed points (0, K) or (K, 0), from Equation ( 24) the well depth is proportional to carrying capacity K, being In a Kramers' type approximation, the escape time from the pseudo-potential well scales as ∼ exp(∆U ) [46], reproducing the exponential scaling of the extinction time with K, observed numerically.Moreover, the Fokker-Planck approximation also shows that the exponential scaling disappears as niche overlap a approaches unity, in accord with the numerical results in the main text.The correlation between the two species goes to negative one in this parameter limit, such that they are entirely anti-correlated.Whereas the well has a single lowest point at the co-existence fixed point for partial niche overlap, at a = 1 the potential shows a trough of equal depth going between the two axial fixed points.This is the Moran line, along which diffusion is unbiased; diffusion away from the Moran line is restored, as the system is drawn toward the bottom of the trough.Because everywhere along the Moran line is equally likely, the probability cannot be normalized, and the linearization approximation fails.This is to be expected, as it is an expansion about a fixed point, but the fixed point is replaced by the Moran line in the Moran limit of a = 1.
Breaking the parameter symmetries Left: Dashed lines denote the functions in the τ = e h(a21) K g(a21) e f (a21)K ansatz obtained from the fit to the numerical data.come from fitting the ansatz to generated data.Right: The right panel compares the results of the ansatz fit with Kramers'/Fokker-Planck estimate of the fixation time.
The main text treats the symmetric case of K 1 = K 2 ≡ K and a 12 = a 21 ≡ a.Here, we extend our results to the asymmetric case, where the symmetry between the parameters is broken.Although the exponential scaling of the fixation time with the system size (except at the Moran line), persists also in the asymmetric case although the exponential dependence can be much weaker in some of the asymmetric cases.
Figure 4 shows the dependence of the fixation time on the niche overlap a 21 while keeping a 12 = 0.5 for K 1 = K 2 ≡ K, using the similar ansatz we apply the same τ = e h(a21) K g(a21) e f (a21)K .As the niche overlap a 21 changes from 0 to 1, the location of the co-existence fixed point shifts from (K/2, K/4) to (K, 0).Accordingly, the fixation time starting from the fixed point maintains its exponential scaling with carrying capacity up until a 21 = 1, where the fixed time is equal to zero, as reflected in the shape of the of h(a 21 ).Notably, in the asymmetric case the exponential scaling function f (a 21 ) is much weaker compared to the symmetric case, partially because the fixed point is located closer to an axis that in the symmetric case even at a 21 = 0.The right panel of Figure 4 shows the comparison of the results of the ansatz fit with the estimates of the exponential part of the fixation time using Kramers'/Fokker-Planck pseudo-potential described in the previous section that explains the observed trends of f (a 21 ).Next let us consider breaking the symmetry such that the Moran line can still be recovered.The carrying capacity symmetry is broken, such that K 2 = 2K 1 .The two species are still independent when a 12 = a 21 = 0, but in this case the Moran line exists when a 12 = 2 and a 21 = 1/2.Figure 5 shows the results when the symmetry is broken both in the carrying capacity and the niche overlap.It shows the change in the fixation time as a function of the niche overlap a 21 for K 2 = 2K 1 ≡ K while the niche overlaps change along the line where a 12 = 4a 21 , starting from the independent case a 12 = a 21 = 0 to a 12 = 2 and a 21 = 1/2 where the system reaches its corresponding Moran line.The observed behaviour is very similar to that shown in the symmetric case, with the exponential dependence transitioning smoothly to zero at the Moran line.
I uphold the conclusion that only at the Moran line will fixation be fast; when the system parameters are even slightly off those niche overlap values which balance the carrying capacities and allow for the Moran line to exist, the fixation is exponential in the carrying capacity, to the point that the two species effectively co-exist.

Comparison with the Gillespie algorithm
In this section, we verify the numerical results for the mean fixation times obtained using the master equation approach with the truncated transition matrix, agrees with the exact sampling of the fixation process using Gillespie algorithm [47], as shown in Figure 6.fixating.To ensure accuracy of the mean times to 0.1% or better we have chosen the cutoff value C K = 5K although this is largely excessive and even C K = 2K is sufficient for accurate calculation of the fixation times all but the smallest carrying capacities.

Map of a, K showing Co-existence versus Fixation
In the main text we state that since biological system sizes are typically large, a fixation time that scales exponentially with carrying capacity effectively implies co-existence.However, some systems have only a few competing members, as in nascent cancers or plasmids in a single cell.We want to get a better sense of when the exponential scaling is relevant, especially since for those systems with almost complete niche overlap the exponential scaling is slow.To this end we compare the expected mean fixation time with that of the Moran model.The ansatz e h(a) K g(a) e f (a)K is fit to the data and then used to estimate the fixation time at a variety of parameter values.This time is compared to the fixation time of a Moran model with the same carrying capacity.In figure 7 the shaded region represents those parameter combinations for which the estimated fixation is faster than the corresponding Moran model.As is evident, a carrying capacity of forty is fully sufficient to allow for effective co-existence of two species which are not identical in their niches.Even for systems with a smaller carrying capacity, unless the two species are similar they are expected to co-exist for long times before fixation.The odd curvature at K = 5 comes from an extrapolation of the ansatz to low numbers; for a system with such a small carrying capacity, the simplifying assumptions underlying the model are expected to break down.

Invasion into a one-dimensional deterministic logistic model
In this section we estimate the time a deterministic single logistic system would take to grow from one organism close to its carrying capacity.
x o (K − x f ) .
If we assume x 0 = 1, K 1, and x f = (1 − )K then this becomes If we further assume that 1 we can conclude and so expect the invasion time to grow logarithmically with carrying capacity.These results are well-known in the literature [41,48].

Figure 2 :
Figure 2: Left: Dependence of the fixation time on carrying capacity and niche overlap.The lowest line, a = 1, recovers the Moran model results with the fixation time algebraically dependent on K for K 1.For all other values of a, the fixation time is exponential in K for K 1. Right: Niche overlap controls the transition from co-existence to fixation.Blue line: f (a) from the ansatz of Equation (3.6) characterizes the exponential dependence of the fixation time on K; it smoothly approaches zero as the niche overlap reaches its Moran line value a = 1.Green line: g(a) quantifies the scaling of the pre-exponential prefactor K g(a) with K. Yellow line: h(a) is the multiplicative constant.The bars represent a 95% confidence interval, and for f are thinner than the line.The dots at the extremes a = 0 and a = 1 are the expected asymptotic values.Red line: Gaussian approximation to Fokker-Planck equation (see next subsection).

Figure 4 :
Figure 4: Probability of a successful invasion.Left: Solid lines show the numerical results, from a = 0 at the top to a = 1 at the bottom.The purple solid line is the expected analytical solution in the independent limit.The green solid line is the prediction of the Moran model in the complete niche overlap case.Data come from Equation (3.7) and are connected with dotted lines to guide the eye.Right: The red data show the results for carrying capacity K = 4, and suggest the solid black line bmut

Figure 1 :
Figure 1: Stability phase diagram of the co-existence fixed point for K 1 = K 2 = K.The co-existence fixed point C of Equations (5) is stable in the green region and unstable in the blue region; in the white regions it is non-biological.Colored dots indicate the parameter range studied in the paper.The numbered regions correspond to different biological different regimes; see main text.Regions 4-6 correspond to competitive exclusion, with only single species fixed point A or B being stable (or both, in the bistable regime 5).In region 7 the populations experience unbounded growth.For the degenerate case a 12 = a 21 = 1, indicated by the red dot, the co-existence fixed point is replaced by a line of marginal stability, which we call the Moran line.

Figure 2 :
Figure 2: The coupled logistic model agrees with the Moran model in the limit of complete niche overlap, a = 1.The fixation time of the Moran model is shown in red; black are the numerical results of the the coupled logistic model for a = 1 using master equation approach.The population size of the Moran model is set equal to the carrying capacity of the corresponding coupled logistic model; K = 64.

1 2Ke
K in the main text.

Figure 3 :
Figure 3: Extinction time distribution of the logistic model is dominated by a single exponential tail.The bulk of the probability density is well described by an exponential distribution with the same mean, shown in the red dotted line.Data are generated using using the Gillespie algorithm for K = 16.For higher carrying capacities the assumption of exponentially distributed times becomes even more accurate.

Figure 4 :
Figure 4: Breaking the symmetry in a.Left: Dashed lines denote the functions in the τ = e h(a21) K g(a21) e f (a21)K ansatz obtained from the fit to the numerical data.come from fitting the ansatz to generated data.Right: The right panel compares the results of the ansatz fit with Kramers'/Fokker-Planck estimate of the fixation time.

Figure 5 :
Figure 5: Breaking the symmetry in K.As in Figure 2 in the main text, lines come from fitting the ansatz to generated data.The exponential dependence is non-zero except at the appearance of the Moran line at a 21 = 1/2.The extreme points are the expected asymptotic values.

Figure 6 :
Figure6: Directly solving the (truncated) master equation agrees with Gillespie simulations.Solid lines come from directly solving the backwards master equation by inverting the transition matrix, after a cutoff has been applied to the matrix to make it finite.Dashed lines are each an average of a hundred realizations of the stochastic process, as simulated using the Gillespie algorithm.

Figure 7 :
Figure 7: Parameter space in which fixation is fast.The white area shows where two species are expected to effectively co-exist, while the black shading identifies the regime where fixation is faster than a similar Moran model.Fixation is estimated by extrapolating the ansatz parameter fits to the a, K parameter space.