Does deterministic coexistence theory matter in a finite world?

Contemporary studies of species coexistence are underpinned by deterministic models that assume that competing species have continuous (i.e. non-integer) densities, live in infinitely large landscapes, and coexist over infinite time horizons. By contrast, in nature species are composed of discrete individuals subject to demographic stochasticity, and occur in habitats of finite size where extinctions occur in finite time. One consequence of these discrepancies is that metrics of species coexistence derived from deterministic theory may be unreliable predictors of the duration of species coexistence in nature. These coexistence metrics include invasion growth rates and niche and fitness differences, which are now commonly applied in theoretical and empirical studies of species coexistence. Here we test the efficacy of deterministic coexistence metrics on the duration of species coexistence in a finite world. We introduce new theoretical and computational methods to estimate coexistence times in stochastic counterparts of classic deterministic models of competition. Importantly, we parameterized this model using experimental field data for 90 pairwise combinations of 18 species of annual plants, allowing us to derive biologically-informed estimates of coexistence times for a natural system. Strikingly, we find that for species expected to deterministically coexist, habitat sizes containing only tens of individuals have predicted coexistence times of greater than 1, 000 years. We also find that invasion growth rates explain 60% of the variation in intrinsic coexistence times, reinforcing their general usefulness in studies of coexistence. However, only by integrating information on both invasion growth rates and species’ equilibrium population sizes could most (> 99%) of the variation in species coexistence times be explained. This integration is achieved with demographically uncoupled single species models solely determined by the invasion growth rates and equilibrium population sizes. Moreover, because of a complex relationship between niche overlap/fitness differences and equilibrium population sizes, increasing niche overlap and increasing fitness differences did not always result in decreasing coexistence times as deterministic theory would predict. Nevertheless, our results tend to support the informed use of deterministic theory for understanding the duration of species coexistence, while highlighting the need to incorporate information on species’ equilibrium population sizes in addition to invasion growth rates.


Introduction
Understanding how competing species coexist is a central problem in ecology [Hutchinson, 1961, Chesson, 2000a. Recent theoretical progress on this problem has replaced vague conceptualizations of coexistence requirements with tools that allow ecologists to quantify the niche differences that allow coexistence, and the fitness differences that drive competitive exclusion [Chesson, 1990, 2000a, Adler et al., 2007. This progress has, in turn, motivated a large number of empirical studies to apply tools derived from theory to quantify the drivers of species coexistence in the field [e.g., Hille Ris Lambers, 2009, Narwani et al., 2013]. Importantly, however, the deterministic theory on which these advances are based assumes that competition occurs between populations whose densities vary continuously on landscapes of infinite size [Faure andSchreiber, 2014, Schreiber, 2017]. Under these assumptions, the influence of processes occurring as a consequence of the discrete nature of individuals-such as demographic stochasticity-are excluded [Hart et al., 2016, Pande et al., 2020b. This generates a fundamental disconnect between theory and reality because while theory predicts that coexisting species will coexist indefinitely, in nature coexistence can only occur over finite periods of time. How well metrics derived from deterministic coexistence theory predict the duration of coexistence in the discrete, finite systems of nature -the ultimate object of study -remains largely unknown.
One of the most widely used metrics in contemporary theoretical and empirical studies of species coexistence are invasion growth rates [Hofbauer and Sigmund, 1998, Chesson, 2000a, Schreiber, 2000, Grainger et al., 2019. A metric from deterministic models, the invasion growth rate of a species is its per-capita growth rate at vanishingly low densities when its competitors' densities are at equilibrium. For two competing species, deterministic models predict that coexistence occurs when each species in a competing pair has a positive long-term invasion growth rate [Macarthur and Levins, 1967]. Provided there are no Allee effects or positive frequency-dependence at low densities [Schreiber et al., 2019], meeting this "mutual-invasibility" criterion implies that coexistence occurs indefinitely. Invasion analyses have been particularly powerful in studies of species coexistence because modelspecific expressions for the invasion growth rate can be used to derive expressions that quantify the magnitude of niche and fitness differences between species [Chesson, 1990, 2013, Hart et al., 2018. Thus, the mutual invasibility criterion has become central to much of our current understanding of species coexistence. However, when finite populations with positive deterministic invasion growth rates are depressed to low numbers of discrete individuals, they may still fail to persist because the negative effects of demographic stochasticity lead to extinction [Hart et al., 2018, Pande et al., 2020b. More generally, even when species have higher densities away from the invasion boundary, demographic stochasticity operating on competing populations of finite size ensures extinction in finite time [Reuter, 1961, Lande et al., 2003, Jagers, 2010, Schreiber, 2017.
Despite the dominance of coexistence theory that imposes continuously varying population densities and infinitely large landscapes on a discrete, finite world, the effect of demographic stochasticity on coexistence is receiving greater attention [Vellend, 2010, Shoemaker et al., 2020. Perhaps the most prominent example is neutral theory, in which demographic stochasticity is the sole driver of competitive dynamics between species [Hubbell andFoster, 1986, Hubbell, 2001]. Like fixation times for neutral alleles [Ewens, 2012], neutral theory predicts that competing species can coexist for long periods of time when the community size is large relative to the number of species [Hubbell, 2001]. The well-known problem with neutral theory, however, is that in emphasizing the role of demographic stochasticity, the theory simultaneously excludes all deterministic processes governing species interactions [Vellend, 2010]. This exclusion includes species-level fitness differences that drive competitive exclusion, and niche differences that promote coexistence [Adler et al., 2007]. Thus, depending on the relative strength of these deterministic processes, coexistence times predicted by neutral theory (i.e. by the action of demographic stochasticity alone), will be either significant over-or under-estimates.
More recently, population theoretic studies that incorporate demographic stochasticity into traditionally deterministic models of competition have begun to emerge [Adler and Drake, 2008, Gómez-Corral and López García, 2012, Gabel et al., 2013, Kramer and Drake, 2014. These studies first make the important point that interspecific competitive interactions cause extinction dynamics to be different to that predicted by single-species models, and they also demonstrate that stochasticity causes the identity of the winner in competition to be less-than-perfectly predicted by deterministic competition-model parameters [Gómez-Corral and López García, 2012, Gabel et al., 2013, Kramer and Drake, 2014. In addition, and together with models of extinction processes applied to single-species dynamics, these studies also highlight the likely importance of variables not traditionally considered in contemporary studies of species coexistence [Gómez-Corral andLópez García, 2012, Kramer andDrake, 2014]. For example, in single species models, time to extinction depends critically on equilibrium population size [Boyce, 1992, Grimm and Wissel, 2004, Ovaskainen and Meerson, 2010. If equilibrium population size is similarly important for the duration of species coexistence, then the sole focus on invasion growth rates as the primary arbiter of coexistence may be problematic.
Despite recent progress in this area, existing studies of the effects of demographic stochasticity on the duration of coexistence concentrate only on cases where competitive exclusion is deterministically ensured [Gómez-Corral andLópez García, 2012, Kramer andDrake, 2014]. What is missing, therefore, are assessments of the duration of coexistence when coexistence rather than exclusion is deterministically ensured, and how these durations relate to existing deterministic coexistence metrics (e.g. invasion growth rates, niche and fitness differences). The latter point is particularly important as long as the field continues to rely on deterministic theory to identify coexistence mechanisms, and to interpret empirical results [Ellner et al., 2019]. Indeed, Pande et al. [2020b] identify inadequacies in invasion growth rates derived from deterministic theory for predicting coexistence times of finite populations in fluctuating environments. They show that the environmentally-dependent distribution of species' growth rates in fluctuating environmentsnot just the average 'invasion growth rate' -influences coexistence times. While this finding is important, it leaves the independent effects of environmental and demographic stochasticity on coexistence times unresolved. For example, for coexistence mechanisms not reliant on environmental fluctuations there is no distribution of environmentally-dependent population growth rates, and so it remains unclear how demographic stochasticity alone influences coexistence times. Moreover, there is considerable empirical attention on, and support for, the contribution to coexistence of fluctuation-independent coexistence mechanisms [Adler et al., 2010, Chu and Adler, 2015, Ellner et al., 2016, Mordecai et al., 2016, Letten et al., 2017, Li et al., 2019, Muller-Landau and Visser, 2019, Wainwright et al., 2019, Zepeda and Martorell, 2019, Spaak and De Laender, 2020. For example, Zepeda and Martorell [2019] found that coexistence of 17 grassland species is mostly due to large fluctuation-independent coexistence mechanisms. Hence, it is particularly important to understand the applicability of these fluctuation-independent mechanisms for communities influenced by demographic stochasticity in nature.
In this paper, we explore the duration of species coexistence in discrete, finite natural systems experiencing demographic stochasticity. In particular, we determine the relationship between the duration of coexistence in the presence of demographic stochasticity and commonly used metrics of coexistence derived from deterministic theory. To determine these relationships, we introduce the concept of the intrinsic coexistence time, a multispecies analogue of Grimm and Wissel [2004]'s intrinsic extinction time for single-species stochastic population models. This metric corresponds to the mean time to losing one or more of these species after the species were coexisting sufficiently long to exhibit relatively stationary population dynamics i.e. quasi-stationarity. Our assessment of this metric is based on novel mathematical and computational methods that allow us to derive explicit relationships between the quasi-stationary behavior of a stochastic model and the dynamics of its deterministic model counterpart [Faure and Schreiber, 2014]. Importantly, we ground our analytical approach using estimates of competition model parameters from 90 pairs of annual plant species competing on serpentine soils in the field .

Models and Methods
Models and Methods Overview Our models and methods are composed of two broad parts. First, we describe the model and its empirical parameterization, we introduce metrics of deterministic coexistence, and we describe new analytical methods for calculating coexistence times in the presence of demographic stochasticity. Second, we describe our methods for answering a series of questions about the relationship between deterministic coexistence and stochastic coexistence times using the empirically-parameterized model.

Model
We base our analysis on an annual plant competition model [Beverton and Holt, 1957, Leslie and Gower, 1958, Watkinson, 1980, which is well studied analytically [Cushing et al., 2004], and does a good job of describing competitive population dynamics in plant communities in the field . In the deterministic model with no demographic stochasticity, the dynamics of species i (i = 1 or 2) can be expressed in terms of its density n i,t and its competitor's density n j,t at time t: where λ i describes offspring production in the absence of competition, and α ii and α ij are the competition coefficients, which describe the rate of decline in offspring production as conspecific and heterospecific competitor densities increase, respectively. Including the multiplicative factor λ i − 1 in the denominator ensures that the competition coefficients are in units of proportional reductions in λ i , such that when α ii n i + α ij n j = 1 the population stops growing. Parameterizing our model in this way ensures that the conditions for coexistence in our model correspond to the classical conditions for coexistence in the continuous-time Lotka-Volterra competition model. In the deterministic model (1), offspring production and population densities take on values in the non-negative, real numbers. The stochastic model takes the same functional form, but in contrast to the deterministic model, offspring production and population sizes N i,t take on non-negative integer values. Population densities N i,t /S are determined by the community size S and take on non-negative rational values. Specifically, in the stochastic model discrete individuals produce random numbers of discrete S= 10 D Figure 1: In panels A-C, the deterministic dynamics of the species densities n i,t for (1) are shown in light gray and light blue. Sample trajectories of the stochastic dynamics of the species densities N i,t /S for (2) are shown in blue and black. As the habitat size S increases, there is a closer correspondence between these dynamics, but eventually species go extinct in the stochastic model. In panel D, a sample simulation of the Aldous et al. [1988] algorithm to estimate the quasi-stationary distribution and the intrinsic coexistence time. Red circle correspond to times at which one or both species are going extinct; states at these times are replaced by states randomly sampled from the past. The frequency p of these red dots, for sufficiently long runs, determines the intrinsic coexistence time 1/p. offspring according to a Poisson distribution. These random individual-level reproductive events generate demographic stochasticity in our model. Assuming the Poisson-distributed offspring production of individuals are independent with the same mean as the deterministic model, the sum of these events is also Poisson distributed and the dynamics of our stochastic model at the population level can be expressed as: where Pois(µ) denotes a Poisson random variable with mean µ. The community size parameter S allows us to quantify the effects of demographic stochasticity on coexistence in landscapes of finite size. Importantly, when community size S is sufficiently large, the dynamics of the densities N i,t /S of our stochastic model (2) are well approximated by the deterministic model (1) (Fig. 1A-C). This justifies our analysis of the effects of metrics of deterministic coexistence on stochastic coexistence times. Specifically, for any prescribed time interval, say [0, T ], N i,t /S is highly likely to remain arbitrarily close to n i,t provided that S is sufficiently large and N i,0 /S =n i,0 (Appendix S2). Intuitively, for fixed initial densities of both species, larger community sizes correspond to greater population sizes and, consequently, smaller stochastic fluctuations in their densities. However, in contrast to the deterministic model, both species eventually go extinct in finite time in the stochastic model (Appendix S2).
Model parameterization In our analyses we use parameter values for our competition model that were estimated for 90 pairwise combinations of 18 annual plant species competing in experimental field plots on serpentine soils in California . We used established methods to estimate these parameters [Hart et al., 2018]. Briefly, prior to the growing season we sowed focal individuals of each species into a density gradient of competitors. By recording the fecundity of each focal individual near the end of the growing season within a competitive neighborhood, we were able to estimate per germinant seed production in the absence of competition, and the rate of decline in seed production as competitor density increases [Hart et al., 2018]. We also measured the germination rate of each species and seed survival rate in the seedbank. See  for full details of our empirical and statistical methods.
These experiments were used to parameterize a model of the form where n i,t+1 is the density of seeds in the soil after seed production but prior to germination, r i is the per germinant seed production in the absence of competition, a ii and a ij are the competition coefficients for germinated individuals, g i is the fraction of germinating seeds and s i is seed survival. As including the seed bank greatly complicates our analysis, yet ignoring it would unfairly bias competitive outcomes in the system, we assume that seeds that ultimately germinate do so in the year after they are produced. This implies that higher seed survival rates ultimately increase the average annual germination rate,g i (Appendix S1). After making this adjustment, model (3) is equivalent to model (1) after setting s = 0, λ i = r igi and α ij =g j a ij /(λ i − 1). In the absence of a fluctuating environment (where seeds survival may buffer populations from extinction), this change is expected to have no effect on the duration of species coexistence Metrics of deterministic coexistence Ultimately, we wish to relate coexistence times to metrics of deterministic coexistence commonly used in contemporary studies of species coexistence.
Here, we focus on the invasion growth rate, and quantitative definitions of niche and fitness differences, which are themselves derivatives of the invasion growth rate. For (1), the invasion growth rate for species i is given by .
where 1/α jj is the single-species equilibrium density of species j. In the deterministic model, coexistence occurs if the invasion growth rates of both species are greater than one. Equivalently, this occurs when the minimum of the invasion growth rates, min{I 1 , I 2 }, is greater than one. Conversely, if min{I 1 , I 2 } < 1, then one of the species' densities will approach zero over the infinite time horizon. Each species' invasion growth rate is determined by: 1) frequency-dependent processes that provide growth-rate advantages to both species when they are at low relative density, and 2) frequencyindependent processes that always favor the growth of one species over another regardless of relative density. These contributions to invasion growth rates have been quantitatively formalized as the niche overlap (ρ) and the average fitness ratio κ j κ i of the two species, respectively (see [Chesson, 2013] for mathematical details): Niche overlap, ρ, decreases as the strength of intraspecific competition increases relative to the strength of interspecific competition. Low niche overlap causes species at high relative density to limit their own growth through intraspecific competition more than they limit the growth of species at low relative density through interspecific competition, leading to higher invasion growth rates. The fitness ratio κ j κ i increases as one species becomes more sensitive to the total amount of competition in the system. For example, when κ j κ i is greater than one, species i is more impacted by competition and will be displaced by species j when there is perfect niche overlap. The fitness ratio is used to quantify fitness differences: the greater the magnitude of the log fitness ratio, the greater the fitness difference. Ultimately, both species have positive invasion growth rates when ρ < κ j /κ i and ρ < κ i /κ j , demonstrating that deterministic coexistence occurs when niche overlap is less than 1 and small relative to the fitness difference.
If the species deterministically coexist they will approach a globally stable equilibrium, (n * 1 , n * 2 ). Solving for the densities at which both species' fitness equals 1 in the deterministic model (1) gives Invasion growth rates (I i ) and equilibrium densities (n * i ) can also be expressed in terms of the niche overlap and the fitness ratio: These expressions will become informative for interpreting the effects of niche overlap and fitness differences on coexistence times.
Estimating intrinsic coexistence times Unlike in the deterministic model, in the stochastic model both species always go extinct in finite time. We define the length of time prior to the first species going extinct as the coexistence time. The distribution of coexistence times in the stochastic model will depend on the initial conditions of the system. In studies of single-species persistence, Grimm and Wissel [2004] defined a quantity, the intrinsic mean time to extinction, that provides a common approach for selecting the initial distribution of the population and, thereby, allows for cross-parameter and cross-model comparisons. Here, we extend this work to introduce an analogous concept, the intrinsic coexistence time. The intrinsic coexistence time assumes that competing species have coexisted for a sufficiently long period of time in the past to exhibit relatively stationary population dynamics i.e. quasi-stationarity. At quasi-stationarity, the distribution of community population sizes is given by the model's quasi-stationary distribution (QSD) [Méléard and Villemonais, 2012]. While in this quasi-stationary state, there is a constant probability, call it the persistence probability p, of losing one or both species each time step. We call the mean time 1/p to losing at least one of the species, the intrinsic coexistence time.
More precisely, the QSD for our stochastic model is a probability distribution π(N i , N j ) on positive population sizes N i > 0, N j > 0 such that if the species initially follow this distribution (i.e. Pr [N i , then they remain in this distribution provided neither goes extinct i.e.
By the law of total probability, the persistence probability p when following the QSD satisfies From a matrix point of view, π corresponds to the dominant left eigenvector of the transition matrix for the stochastic model and p is the corresponding eigenvalue [Méléard and Villemonais, 2012].
As the state space for the stochastic model is large even after truncating for rare events, directly solving for the dominant eigenvector is computationally intractable. Fortunately, there is an efficient simulation algorithm for approximating the QSD due to Aldous et al. [1988]. This algorithm corresponds to running a modified version of the stochastic model ( Fig. 1D). Whenever one of the species goes extinct in the simulation, the algorithm replaces this state with a randomly sampled state from the past. More precisely, if (Ñ i (1),Ñ j (1)), . . . , (Ñ i (t),Ñ j (t)) are the states of the modified process until year t, then compute (Ñ i (t + 1),Ñ j (t + 1)) according to (2). IfÑ i (t + 1) = 0 or N j (t + 1) = 0, then replace (Ñ i (t + 1),Ñ j (t + 1)) by randomly choosing with equal probability 1 t amongst the prior states approaches the QSD as t → ∞. This algorithm converges exponentially fast to the QSD [Benaïm and Cloez, 2015]. We use this algorithm to estimate coexistence times in our analyses. In Appendix S8, we also show how this algorithm also applies to a large class of models simultaneously accounting for demographic, environmental stochasticity, and spatial heterogeneity.
Having introduced our model and methods for determining coexistence times, we now describe how we apply these methods to address a series of questions about the relationship between deterministic coexistence and stochastic coexistence times for our empirically-parameterized models. We emphasize that while our models are empirically parameterized and so grounded in real biology, the focus of our analysis is a comparison between the predictions of the empirically-parameterized deterministic model and the predictions of the empirically-parameterized stochastic model. Understanding the relationship between our predicted coexistence times and observed coexistence times in nature is exceedingly difficult, likely requiring empirical studies over hundreds to thousands of generations in most systems. This is beyond the scope of our current analysis.
How do deterministic coexistence and competitive exclusion relate to coexistence times? We first explore how deterministic coexistence and deterministic competitive exclusion are related to coexistence times in the presence of demographic stochasticity. We do this through a mixture of analytical and numerical approaches. The analytical approach uses large deviation theory [Faure and Schreiber, 2014] to characterize how intrinsic coexistence times scale with community size S and how this scaling depends on whether the deterministic model predicts coexistence or exclusion. Based on our empirical parameterizations, competition between 82 species pairs is expected to result in deterministic competitive exclusion (i.e. parameter values result in min{I 1 , I 2 } < 1), while 8 species pairs are expected to stably coexist (i.e. parameter values result in min{I 1 , I 2 } > 1). The remaining species pairs were excluded from our analyses either because estimates for at least one of the model parameters were missing or because one of the species had an intrinsic fitness λ i of less than one. For the two groups of species pairs we focus on (i.e. either deterministically coexisting or resulting in deterministic exclusion), we used simulations of ten million years to compute intrinsic coexistence times C across a range of community sizes S. We describe the relationship between community size and coexistence times for the models of the species pairs in each group.
Do invasion growth rates predict coexistence times? We used linear regression to determine if deterministic invasion growth rates (4) influence intrinsic coexistence times. For this analysis, we calculated intrinsic coexistence times only for the species pairs that are expected to deterministically coexist. To estimate intrinsic coexistence times for these species, we set the community size S = 0.04. This size would be empirically justified for the smallest serpentine hummocks in our study landscape (those with a few m 2 of suitable habitat), which might contain as few as 20 individuals of the subdominant species (based on germinable densities projected from Gilbert and Levine [2013]). The species that compose the focal pairs in this analysis tend to be more common and are often found on larger hummocks, but this small habitat size allows us to evaluate the dynamics of systems where the effects of demographic stochasticity might be substantial. For each species pair we used simulations of 10 7 years to estimate the coexistence times predicted by the models. For seven of the species pairs predicted to coexist by the deterministic models, we were able to generate good estimates of intrinsic coexistence times predicted by the stochastic models because there were multiple extinction events in simulations of this length. For one species pair predicted to coexist by the deterministic model, there were no extinction events in the numerical simulations. Because we only have a lower bound of > 10 7 years for the mean intrinsic coexistence time for this pair, it was excluded from our analysis.
Coexistence time was the dependent variable in our regression, and the log of min{I 1 , I 2 } was the independent variable. We used this independent variable because the species with the lower invasion growth rate is, all else being equal, more likely to go extinct first. To examine the robustness of our conclusions, we repeated our analysis for 1, 000 randomly-drawn parameter values with community size S = 10. These random draws were performed on uniform distributions with [1.1, 2] for the λ i values, [0.2, 0.5] for the α ii values, and α ii × [0, 1] for the α ji values (to ensure deterministic coexistence).
Do equilibrium population sizes predict coexistence times? Invasion growth rates were a less-than-perfect predictor of coexistence times for the stochastic models (see Results). Therefore, we also used linear regression to explore the relationship between the equilibrium population sizes of the coexisting species and the intrinsic coexistence time. For this analysis, we used the same methods as described above, but with the log of the minimum of the two equilibrium population sizes i.e. min{Sn * 1 , Sn * 2 } as the independent variable in the linear regression. To examine the robustness of our conclusions, we repeated our analysis for 1, 000 randomly-drawn parameter values with community size S = 10 as described in the section Do invasion growth rates predict coexistence times?, above. Do greater niche overlap or greater fitness differences always reduce coexistence times? According to deterministic theory, greater niche overlap and greater fitness differences both have a negative impact on coexistence [May, 1975, Chesson, 2000a, Adler et al., 2007. It is important, therefore, to understand if these negative effects on deterministic coexistence also have predictably negative effects on the duration of species coexistence. We investigate these relationships by calculating coexistence times as a function of niche overlap and fitness differences.
Our goal is to assess the independent effects of these determinants of coexistence on coexistence times. However, because niche overlap and fitness differences are both functions of the same parameters (5), they are not quantitatively (or biologically) inherently independent [see also Song et al., 2019]. This precludes using the natural variation in niche overlap and fitness differences observed across our focal species pairs for our analyses. Therefore, to achieve our goal, we manipulated niche overlap and fitness differences separately for each species pair. Specifically, to assess the effects of niche overlap on model predicted coexistence times, we multiplied the interspecific competition coefficients α 12 , α 21 for each coexisting species pair by a common, fixed factor. This manipulation allows the niche overlap to vary while keeping the fitness ratio constant according to (5). Mechanistically, this may be interpreted as altering the degree to which the two species use the same resources. Similarly, to assess the effects of the fitness difference independent of any change in niche overlap, we multiplied the competition coefficients α 22 , α 21 within each coexisting pair by a fixed factor. This manipulation reduces the sensitivity of species 2 to competition, increasing its competitive ability according to (5), while having no effect on niche overlap. Mechanistically, this may be interpreted as increasing the efficiency with which species 2 uses shared resources. Based on these manipulations, for each coexisting species pair we explored the relationship between niche overlap and fitness differences, and coexistence times. We interpret these relationships via the effects of niche overlap and fitness differences on invasion growth rates and equilibrium population sizes as per (7).
Do stochastic competitive dynamics influence coexistence times? Competition may also influence coexistence times because the stochastic population dynamics of the two species are coupled. For example, a stochastic increase in the population size of one species might be expected to result in a concomitant decrease in the population size of its competitor, increasing its risk of extinction. The effect on coexistence times of this dynamic coupling of the two species is not accounted for by static metrics of coexistence such as invasion growth rates and equilibrium population sizes. To identify whether this dynamic coupling plays an important role in determining coexistence times, we built a simplified model, what we call the demographically uncoupled model, incorporating the effects of competition on invasion growth rates and equilibrium population sizes, but excluding the coupling of the stochastic fluctuations in the population sizes of the competitors. In our demographically uncoupled model, each species has a low density growth rate "λ i " and carrying capacity "1/α ii " that matches its invasion growth rate I i and equilibrium density n * i , respectively, from the deterministic two-species model. The update rule for this simplified model is Importantly, note that the change in population size of N i does not depend on N j . Therefore, in the demographically uncoupled model the effects of interspecific competition on invasion growth rates and equilibrium population sizes are retained, but the species are dynamically uncoupled. By quantifying the relationship between coexistence times calculated for the full (2) and demographically uncoupled (8) models, we can determine whether the combined effects of competition on invasion growth rates and equilibrium population sizes are sufficient to predict coexistence times, or if competitive dynamics arising from the coupled stochastic dynamics of the species are also important.
In Appendix S7, we show the coexistence times for the demographically uncoupled model reduce to calculating the intrinsic persistence times P i of each species in this simplified model. These persistence times P i are calculated independently for each species using the [Aldous et al., 1988] algorithm. Then, the coexistence time C simplified of the demographically uncoupled model equals: When the intrinsic persistence times P i for each species are sufficiently long (e.g. 10 years or more), the coexistence time predicted by this simplified model is approximately one half of the harmonic mean of the intrinsic persistence times for the two species: C simple ≈ 1 2 × . As the harmonic mean is dominated by the minimum of its arguments, the species with the shorter persistence time has the greater influence on the coexistence time.
We calculated intrinsic coexistence times for the full and demographically uncoupled models for the eight species pairs that are expected to deterministically coexist. We then used linear regression to determine the extent to which the simplified model explains coexistence times in the full model. To examine the robustness of our conclusions, we repeated our analysis for 1, 000 randomly-drawn parameter values with habitat size S = 10 as described in the section Do invasion growth rates predict coexistence times?, above.

Results
How do deterministic coexistence and competitive exclusion relate to coexistence times? The duration of species coexistence tended to be several orders of magnitude larger for species expected to deterministically coexist than for those pairs where deterministic exclusion is expected (Fig. 2). Indeed, our numerical results demonstrate that species expected to deterministically coexist will tend to do so for decades to millennia, even when community sizes are small. By contrast, our numerical simulations suggest that the vast majority of the 82 species pairs for which deterministic exclusion is expected will tend to coexist for less than five years ( Fig. 2A).
Our analytical and numerical results expose a fundamental dichotomy between deterministic coexistence and exclusion in terms of how coexistence times for the models scale with increasing community size. In particular, if deterministic exclusion is predicted (min{I 1 , I 2 } < 1), then our mathematical analysis (Appendix S3) implies that the intrinsic coexistence times are bounded above by 1/(1 − min{I 1 , I 2 }) and do not increase significantly with community size. This analytical result is consistent with our numerical results for the 82 models associated with the species pairs for which deterministic exclusion was predicted ( Fig. 2A). By contrast, if the two species are predicted to deterministically coexist (min{I 1 , I 2 } > 1), then our analysis implies that intrinsic coexistence times scale exponentially with the community size S. This means that small increases in community size lead to very large increases in the duration of species coexistence. This analytical result is also : Predicting coexistence times using the minimum, min{I 1 , I 2 }, of invasion growth rates (A) and the minimum, min{Sn * 1 , Sn * 2 }, of equilibrium population sizes (B). Both panels plotted on a log-log scale. Each of these species pair curves are colored by their min{I 1 , I 2 } value as in Figure  2.
consistent with our simulations of the 8 species pairs for which deterministic exclusion was predicted (Fig. 2B).
For species pairs predicted to deterministically coexist, the mathematical analysis only ensures there is a constant > 0 such that the intrinsic coexistence time is proportional to e S . There is no simple formula for as it will depend on many details of the individual-based model (see proofs in [Faure and Schreiber, 2014]). To understand some of these dependencies, we explored how coexistence times depend on key coexistence metrics from the deterministic model for a fixed community size S = 0.04, (a small hummock; see methods), the results of which we describe next.
Do invasion growth rates predict coexistence times? Across seven pairs of deterministically coexisting competitors, we found a positive correlation between the minimum invasion growth rate, min{I 1 , I 2 }, and intrinsic coexistence times, as determined from the quasi-stationary distributions of the empirically-parameterized model (adjusted R 2 = 0.6038, p = 0.0244; Fig. 3A). Intuitively, the species with the lower invasion growth rate tends to go extinct first, and so it is the lower of the two invasion growth rates that predicts coexistence times. However, there was substantial variation in coexistence times unexplained by invasion growth rates, such that similar invasion growth rates resulted in several orders of magnitude difference in coexistence times (Fig. 3A). Random sampling of parameter space produced similar results (Appendix S4: Fig. S1A).
Do equilibrium population sizes predict coexistence times? We also found a positive correlation between the minimum of the equilibrium population sizes, min{Sn * 1 , Sn * 2 }, and coexistence times (adjusted R 2 = 0.3386; p = 0.09964, Fig. 3B). As equilibrium population sizes did a worse job of explaining the variation in the coexistence times than the invasion growth rates, similar equilibrium population sizes also resulted in several orders of magnitude of difference in coexistence times. Random sampling of parameter space produced similar results (Appendix S4: Fig. S1B).

Do greater niche overlap or greater fitness differences always reduce coexistence times?
Greater niche overlap and greater fitness differences did not always reduce coexistence times, as deterministic theory would predict (Figs. 4,5). For 5 out of 7 species pairs, both greater niche overlap and greater fitness differences reduced the coexistence times predicted by the models (Figs. 4A,D; Appendix S6). For these pairs, the relationships were nonlinear, but negatively monotonic. However, for the remaining two species pairs, the effects of greater niche overlap and fitness differences on coexistence times were complex (Figs. 5A,D; Appendix S6). For these species pairs, increasing niche overlap consistently reduced coexistence times, but in a highly nonlinear fashion (Fig. 5A). By contrast, increasing fitness differences had both positive and negative effects on coexistence times for the models associated with these species pairs, depending on the magnitude of the change in the fitness difference (Fig. 5D). Importantly, for the species pairs where these complex effects emerged, the competitively inferior species (as detemined by the fitness ratio in equation (5)), had the higher (i.e not the minimum) equilibrium population size. When the inferior competitor has the higher equilibrium population size, complex effects on coexistence times emerge via the influence of niche overlap and fitness differences on the minimum of the equilibrium population sizes and the minimum of the invasion growth rates (see (7)). For example, as fitness differences increase, the superior competitor's equilibrium population size also increases (equation (7), Appendix S6, Fig. 5E,F). Consequently, when the superior competitor has the lower equilibrium population size, increasing the fitness difference increases the minimum of the equilibrium population sizes, which increases coexistence times (Fig. 5D). Increasing niche overlap can also increase the equilibrium population size of the superior competitor, but always decreases its invasion growth rate (equation (7), Appendix S6, Fig. 5B,C). When the superior competitor has the lower equilibrium abundance, these opposing trends result in coexistence times that decrease in a highly nonlinear manner with increasing niche overlap. A general result is that coexistence times decrease when both the minimum invasion growth rate and the minimum equilibrium population size decline with increasing niche overlap and increasing fitness differences (Figs. 4,5).
Do stochastic competitive dynamics influence coexistence times? To isolate the influence of coupled stochastic competitive dynamics on coexistence times we compared coexistence times calculated from the full model (2) including these coupled dynamics, to coexistence times calculated from the demographically uncoupled model (8) excluding these coupled dynamics. As described in the methods, the demographically uncoupled model is two uncoupled, individual-based single species models whose low-density growth rate and intraspecific competition coefficients equal I i and 1/n * i , respectively. For the 7 species pairs predicted to coexist in the deterministic model, the simplified model incorporating competition only through its effects on invasion growth rates and equilibrium population sizes did an exceptional job in predicting the actual coexistence time (log C = 1.031 log C simple with R 2 = 0.9977 and p < 10 −8 , Fig. 6). Random sampling of parameter space produced similar results (Appendix S4: Fig. S1C).  (E), the deterministic equilibrium population sizes Sn * i are plotted with the gray shading indicating the minimum of the two population sizes. In (C) and (F), the deterministic invasion growth rates of the deterministic, mean field model are plotted for the fitness inferior (dashed red) and fitness superior (solid black) with gray shading indicating the minimum of two invasion growth rates.

Discussion
We introduced a new metric, the intrinsic coexistence time, that characterizes the risk of species loss due to demographic stochasticity over any time horizon. This metric complements more traditional coexistence metrics based on invasion growth rates and can be computed for many existing databased models by following a two step procedure. First, extend the deterministic model to account for demographic stochasticity. Informing such a stochastic model with demographic data obtained in a field context is not much more difficult than parameterizing deterministic models. Aside from fecundity, most transition rates (i.e. survival, growth, dispersal probabilities) of the deterministic model directly transfer to the stochastic model. The fecundity distribution can be estimated from the raw data used to calculate mean fecundity in deterministic models, or one can assume that fecundity is Poisson distributed with the calculated mean. Second, estimate each intrinsic coexistence time by a single simulation using the algorithm of Aldous et al. [1988] (see Appendix S8 for extensions to a broad class of models). Using this two-step approach, we evaluated how well metrics  (E), the deterministic equilibrium population sizes Sn * i are plotted with the gray shading indicating the minimum of the two population sizes. In (C) and (F), the deterministic invasion growth rates of the deterministic, mean field model are plotted for the fitness inferior (dashed red) and fitness superior (solid black) with gray shading indicating the minimum of two invasion growth rates. from deterministic theory predict intrinsic coexistence times for 18 species of Californian annuals.
Consistent with the deterministic theory, we find invasion growth rates are, in general, good predictors of intrinsic coexistence times in models explicitly accounting for demographic stochasticity. Thus, reinforcing the general usefulness of these invasion growth rates in theoretical and empirical studies of species coexistence. Their usefulness stems from two conclusions of our study. First, when both species' invasion growth rates are positive (i.e. deterministic coexistence is predicted), intrinsic coexistence times increase exponentially with community size. Thus, for a given minimum invasion growth rate, doubling community sizes quadruples coexistence times. Strikingly, for the 8 species pairs of serpentine annuals predicted to deterministically coexist, community sizes corresponding to only tens of individuals of the less common species were still sufficient to ensure predicted coexistence times of greater than 1, 000 years (Fig. 2) -a time frame well beyond most empirically-relevant scenarios and much greater than considered in most conservation studies [Meine, 1999]. Second, for a given community size, we found that the minimum of the invasion growth rates (min{I 1 , I 2 }) explained 60% of the variation in coexistence times of the serpentine annuals predicted to deterministically coexist (Fig. 3), and 82% of the variation of coexistence times for a random sampling of parameter space (Appendix S4: Fig. S1A). The two conclusions are consistent with recent work on coexistence times for the lottery model which simultaneously account for environmental stochasticity as well as demographic stochasticity [Pande et al., 2020b, Ellner et al., 2020. The main differences being that invasion growth rates for the lottery model in the presence of temporal environmental stochasticity correspond to the geometric mean of fitness across time, with coexistence times increasing as a power law, rather than exponentially, with community size [Ellner et al., 2020]. Despite explaining a significant amount of variation, invasion growth rates are not perfect predictors of coexistence times in models accounting for demographic stochasticity. In particular, predicted coexistence times varied by an order of magnitude for species pairs of serpentine annuals with similar values of min{I 1 , I 2 } (Fig. 3). Thus, while higher invasion growth rates tend to be associated with longer coexistence times, and these coexistence times tend to be exceedingly long, the capacity of invasion growth rates to accurately predict coexist times should still be viewed with some caution. Pande et al. [2020b] raised a similar concern but for different reasons, when studying Chesson's lottery model. They demonstrated that increases in environmental variability can simul-taneously lead to an increase in invasion growth rates and shorter coexistence times. Our study finds that some disassociation between invasion growth rates and coexistence times can occur even without temporal fluctuations in the fitness of the rare species.
Beyond invasion growth rates, we found that population sizes away from the invasion boundary (i.e. equilibrium population sizes) are also important determinants of coexistence times. This finding is consistent with single-species models of extinction risk, where both low-density growth rate and equilibrium population size determine risk of extinction [e.g., Lande et al., 2003]. However, in and of themselves, equilibrium population sizes were poorer predictors of coexistence times than invasion growth rates. The minimum of the equilibrium population sizes explained 26% less of variation than min{I 1 , I 2 } in the empirical species pairs (Fig. 3), and 12% less of the variation of the coexistence times in the random sampling of parameter space (Appendix S4: Fig. S1B).
When combined, the invasion growth rates and the equilibrium population sizes are such good determinants of coexistence times that simply knowing these two quantities can explain over 99% of the variation of the predicted coexistence times (Fig. 6, Appendix S4:Fig. S1C). Notably, this predictive capacity emerges from a demographically uncoupled model parameterized only with the invasion growth rates and equilibrium population sizes, while excluding the effects of the coupled stochastic dynamics of interspecific competition. The success of the demographically uncoupled model emphasizes that coexistence times in the presence of demographic stochasticity can be predicted using minimal information gleaned directly from an appropriate deterministic model. Whether this tractable approach applies to stochastic models with more species or different nonlinearities in the per-capita growth rates is an important direction for future research.
As the predicted coexistence time for the demographically uncoupled model is approximately the harmonic mean of the persistence time for each of the (re-calibrated) single-species models, one can use it to gain additional insights into coexistence times for the competing species. For example, if both competitors have large invasion growth rates (i.e. I i 1), then each competitors' persistence time is approximately equal to its exponentiated equilibrium population abundance, exp(n * i S) (see Appendix S7). In which case, for a given community size n * 1 S + n * 2 S = N , the predicted coexistence time is maximized when both species have equal equilibrium abundances (see Appendix S7). For a more diverse community, this is equivalent to when Shannon's diversity index is maximized. However, if one species has a significantly lower invasion growth rate than its competitor, then maximizing the coexistence time requires that this species be over-represented in the community (i.e. the community would have a lower Shannon diversity index) (see Appendix S7).
Because invasion growth rates are imperfect predictors of coexistence times, coexistence metrics derived solely from invasion growth rates are also imperfect predictors of coexistence times. These coexistence metrics include modern quantitative definitions of niche overlap and fitness differences [Chesson, 2000a, 2013, Chesson, 2018. Because niche overlap and fitness differences can have opposing effects on the relative population sizes of the competitors, coexistence times do not necessarily decrease with increasing niche overlap or increasing fitness differences, as would be expected based only on deterministic theory (Fig. 5A,D). For models of two of the annual, serpentine species pairs, this unexpected outcome occurred when the inferior competitor had the higher equilibrium population size. Operationally, one can simply check whether this is the case before interpreting the effects of niche overlap and fitness differences as they are commonly derived and applied [Chesson, 2000, Adler et al., 2007, Gilbert and Levine, 2013, Chesson, 2013. According to equation (7), the competitive inferior has the larger equilibrium abundance when it is insensitive to intraspecific competition (i.e. has a large carrying capacity 1/α jj ) but highly sensitive to interspecific competition relative to the competitively superior species. This observation provides an indirect means of evaluating whether this inverted relationship is likely to occur in a given empirical system.
Our conclusions are based on an analysis of an annual plant competition model, which is a specific case of a large and more general class of models that can be re-expressed in Lotka-Volterra form. This class of models has been the focus of much of the theoretical and empirical work deriving and applying quantitative definitions of niche and fitness differences to studies of species coexistence [Macarthur and Levins, 1967, Chesson, 1990, 2000a, 2013. Given that niche and fitness differences in these models are derived from invasion criteria, we expect our conclusions about the effects of niche and fitness differences on coexistence times to apply across this class of models. More generally, equilibrium densities of competitors are likely to be important quantities in all models of competition regardless of their complexity, but are not often taken into account in recent assessments of coexistence. [Carroll et al., 2011, Chesson, 2013, Spaak and De Laender, 2020. For example, [Spaak and De Laender, 2020] provided a general method for defining niche and fitness differences using per-capita growth rates evaluated at densities where at least one species is absent. However, these metrics do not take into account, and are unlikely to correlate with, coexistence equilibrium densities. Therefore, coexistence times are unlikely to map intuitively to niche and fitness differences even for these more general definitions. In sum, our work suggests that any coexistence metric not explicitly taking into account densities at which the species coexist deterministically may not predict coexistence times correctly.
For models considered here, all conspecific individuals have the same invasion growth rate across time -there is no variation in the invasion growth rates among subpopulations of individuals. However, phenotypic, spatial, and temporal variation in vital rates can impact competitive outcomes [Levin, 1974, Warner and Chesson, 1985, Chesson, 1994, 2000b, Schreiber et al., 2011, Vasseur et al., 2011, Hart et al., 2016, Stump et al., 2021. Variation in vital rates often results in variation in invasion growth rates among subpopulations of individuals within generations (spatial or phenotypic variation) or between generations (temporal variation). In these situations, the distribution of these invasion growth rates, not only the mean, can play a role in coexistence times. Understanding their role and how coexistence times depend on the nature of the variation (phenotypic vs. spatial vs. temporal) is a major challenge for future work. As our methodology for computing intrinsic coexistence is applicable to models accounting for this variation (see details in Appendix S8), it may provide a unified approach to tackling this challenge.
Some progress has been made on understanding the role of temporal variation on coexistence times [Pande et al., 2020b, Ellner et al., 2020, Pande et al., 2020a. For example, using the methods presented here, [Ellner et al., 2020] found that the mean invasion growth rates for the lottery model [Warner and Chesson, 1985] are strongly correlated with intrinsic coexistence times. However, [Pande et al., 2020b,a] showed that ignoring the temporal variation in invasion growth rates can result in incorrect inferences about coexistence times. For example, greater environmental variation can increase mean invasion growth rates via the storage effect [Chesson, 1994], but can also lead to shorter coexistence times via the storage effect. Similar conclusions may also apply to within generational sources of variation. For example, increasing within-generational variation in fecundity decreases the persistence times of single species [Melbourne and Hastings, 2008]. We anticipate similar effects on coexistence times in our models.
The complexities we have identified in our study highlight the need to understand deterministic coexistence from both equilibrium-based and invasion-based approaches. Recent relevant equilibrium-based theories include Saavedra et al. [2017]'s structural stability of feasible equilibria and Barabás et al. [2014]'s sensitivity analysis of stable equilibria. Both of these approaches provide insights into the range of demographic parameter values for which the equilibrium abundances remain positive i.e. feasible. Moreover, Saavedra et al. [2017] have developed feasibility metrics analogous to stabilizing niche differences and fitness differences. These developments raise the promising possibility of developing more informative, integrative metrics for species coexistence times. This Appendix is for the manuscript "Does deterministic coexistence theory matter in a finite world?" submitted to Ecology by Sebastian J. Schreiber, Jonathan M. Levine, Oscar Godoy, Nathan J.B. Kraft, and Simon P. Hart Appendix S2 Deterministic approximation and extinction In this Appendix, we prove two results. First, we prove that for any given time frame [0, T ], the probability of stochastic model n t = (n 1,t , n 2,t ) deviating from the mean field modeln t = (n 1,t ,n 2,t ) is arbitrarily small for 0 ≤ t ≤ T provided the habitat size S is sufficiently large. Second, we prove that, despite this strong correspondence over finite time frames, the stochastic models goes to extinction in finite time with probability one. To facilitate the presentation of these conclusions, define f (n) = (f 1 (n), f 2 (n)) where f i (n) = n i λ i /(1 + (α ii n i + α ij n j )(λ i − 1)) where i, j ∈ {1, 2} and i = j. Let M = max{λ 1 /α 11 , λ 2 /α 22 }. Notice that f i (n) ≤ M for all n = (n 1 , n 2 ).
For the first assertion, let T ≥ 1 and ε > 0 be given. Assume that n 0 =n 0 . We will show that Continuity of f implies that there exists δ > 0 such that max 0≤t≤T n t −n t ≤ ε whenever f (n t−1 )− n t ≤ δ for 1 ≤ t ≤ T . As the mean and variance of a Poisson random variable are equal, Chebyshev's inequality implies from which (S-1) follows. For the second assertion, let Next we use the following standard result in Markov chain theory [Durrett, 1996, Theorem 2.3 in Chapter 5].
Proposition. Let X be a Markov chain and suppose that ; the last restriction ensures deterministic coexistence i.e. min{I 1 , I 2 } > 1. In panel A, a linear regression of log 10 C against min{I 1 , I 2 } explained approximately 82% of the variance (adjusted-R 2 of 0.824) with estimated intercept 1.00761 and slope 8.68233. In panel B, a linear regression of log 10 C against min{n * 1 , n * 2 } explained approximately 70% of the variance (adjusted-R 2 of 0.717) with estimated intercept 0.36795 and slope 1.85465. In panel C, a linear regression of log 10 C against log 10 C simple with intercept 0 explained nearly all of the variance (adjusted-R 2 of 0.9963) with estimated slope 1.006993.   This Appendix is for the manuscript "Does deterministic coexistence theory matter in a finite world?" submitted to Ecology by Sebastian J. Schreiber, Jonathan M. Levine, Oscar Godoy, Nathan J.B. Kraft, and Simon P. Hart Appendix S7 Coexistence times for the simplified model To determine the coexistence times for the simplified model, let P i be the mean persistence time for each species. As these times are geometrically distributed, the probability of extinction in each time step is 1/P i for species. As the two species dynamics are independent in the simplified model, the probability that neither goes extinct in a time step is (1 − 1/P 1 )(1 − 1/P 2 ), and the probability at least one goes extinct in a time step is 1 − (1 − 1/P 1 )(1 − 1/P 2 ). Hence, the mean time to losing one of the species is as claimed in the main text. For P 1 1 and P 2 1, we have P 1 P 2 P 1 + P 2 − 1 = 1 1/P 2 + 1/P 1 − 1/(P 1 P 2 ) ≈ 1 1/P 1 + 1/P 2 which equals 1/2 of the harmonic mean of the 1/P i as claimed in the main text.
Using these expressions, we can derive various insights about the two species model. First, we will show that P i ≈ exp(n * i S) whenever I i 1. Indeed, for N i,t ≥ 1 in (8) Thus, whenever N i,t ≥ 1, N i,t+1 ≈ Pois(Sn * i ) and the probability of extinction is approximately exp(−Sn * i ). As these approximations hold whenever N i,t ≥ 1, they describe both the quasi-stationary distribution and quasi-extinction probability. Hence, P i ≈ exp(−n i S)..

S-1
This Appendix is for the manuscript "Does deterministic coexistence theory matter in a finite world?" submitted to Ecology by Sebastian J. Schreiber, Jonathan M. Levine, Oscar Godoy, Nathan J.B. Kraft, and Simon P. Hart Appendix S8 Quasi-stationary distributions for spatially structured communities in fluctuating environments This Appendix describes how the methods used in the main text extend to ecological models accounting for any number of species, any types of species interactions, discrete spatial structure, and environmental fluctuations. After describing the model with and without demographic stochasticity, we present conditions that ensure the existence and uniqueness of a quasi-stationary distribution. Then we describe how the algorithm of Aldous et al. [1988] can be used to numerically estimate the quasi-stationary distributions. We conclude by discussing how the scaling of coexistence times with the community size S depend on the nature of the temporal fluctuations. While we focus on a model using Poisson updates for clarity, the same ideas apply to model with mixtures of multinomials and negative binomials.

Spatially structured community dynamics in fluctuating environments
The mean field model allows for k species and spatial patches. For simplicity, assume that all species have synchronized generations (e.g. competing annual plants, host-parasitoid interactions) and a time step correspond this generation time. Let n ij denote the density of species i in patch j and n i = (n i1 , . . . , n i ) the vector of species i's densities across the patches. Let n = ( n 1 , . . . , n k ) be the state of the community. To model the environmental fluctuations, let ξ 1 , ξ 2 , . . . be an ergodic, stationary sequence of random vectors that represent the fluctuations in the species parameters. For each species i, let A i ( n, ξ) be a non-negative, × projection matrix such that A( n, ξ) n i is species i's state in the next time step when the current community state is n and the current environmental state is ξ. The entries of A i ( n, ξ) correspond to reproduction and dispersal that may depend on the state of the community and the environment. The mean field model is n i,t+1 = A i ( n t , ξ t ) n i,t i = 1, . . . , k (S-1) To account for demographic stochasticity, we assume that births are Poisson distributed and dispersal occurs via multinomial sampling. As discussed in Faure and Schreiber [2014, Section 5.2], these assumptions imply that updates of each species in each patch are given by independent Poisson random variables whose means are given by mean field model. Namely, if S is the community size and N i = (N i,1 , . . . , N i, ) are the abundances of species i across the patches, then the stochastic model is N i,t+1 = Poisson A i ( N t /S, ξ t ) N i,t i = 1, . . . , k (S-2) A standard argument discussed in Schreiber [2017] implies that N t goes to zero (i.e. all species go extinct) in finite time with probability one. In particular, the process enters the extinction set in finite time with probability one where Z k + refers to the non-negative integer lattice of dimension k i.e. abundances of k species in patches.

Quasi-stationary distributions and their estimation
Let C = Z k + \E be the set of coexistence states. For a given community size S, a quasi-stationary distribution for (S-2) is a distribution π S : C → [0, 1] on the non-extinction states (i.e. N ∈C π S (N ) = 1) if there exists λ S ∈ (0, 1) such that N ∈C P [N t+1 = M |N t = N ]π S (N ) = λ S π S (M ). λ S corresponds to the probability of all species coexisting over one generation given the model is in the quasi-stationary state. Hence, 1/(1 − λ S ) is the intrinsic coexistence time. The argument in the proof of Faure and Schreiber [2014, Proposition 6.3] can be used to show that there exists a unique quasi-stationary distribution whenever the mean field model is a compact map (i.e. there exists a compact set K such that n t+1 ∈ K for all n t in (S-1)) and ξ takes values in a compact set.
To numerically estimate the quasi-stationary distribution using the Aldous et al. [1988] algorithm, assume that the stationary sequence ξ 1 , ξ 2 , . . . is given by a Markov process with transition operator T i.e. T (ξ, B) = P (ξ t+1 ∈ B|ξ t = ξ) for any measurable set B of environmental states. To estimate π S , we simulate a modified version of the process ξ t , N t . Namely, given ξ t , N t find ξ t+1 , N t+1 by two steps 1. Draw a multivariate Poisson N t+1 as given by (S-2) and draw ξ t+1 independently from the law of the probability measure T (ξ t , dξ).
For a sufficiently long run of this modified process, the empirical distribution of N t approximates the quasi-stationary distribution i.e. the fraction of time spent the modified process N t spends in community state N approximates π S (N ). Furthermore, the faction of time step 2 is taken approximates the extinction probability 1 − λ S .

Scaling of coexistence times with community size S
When the environmental dynamic is periodic (i.e. ξ 1 , ξ 2 , . . . is periodic) and the mean field model (S-1) has a positive attractor, the argument in the proof of Faure and Schreiber [2014, Theorem 2.7] implies that the coexistence time 1/(1 − λ S ) increases exponentially with the community size S. For stochastic environments (e.g. ξ 1 , ξ 2 , . . . are i.i.d.), the scaling of the coexistence time with community size S is more subtle [Guillin et al., 2019, Ellner et al., 2020. Theoretical and simulation results suggest that the coexistence time scale (at least) polynomially with S whenever the mean field model is stochastically persistent but can be exponential with the positive stationary distribution is bounded away from extinction. S-2