Quantifying invasibility

Abstract Invasibility, the chance of a population to grow from rarity and become established, plays a fundamental role in population genetics, ecology, epidemiology and evolution. For many decades, the mean growth rate of a species when it is rare has been employed as an invasion criterion. Recent studies show that the mean growth rate fails as a quantitative metric for invasibility, with its magnitude sometimes even increasing while the invasibility decreases. Here we provide two novel formulae, based on the diffusion approximation and a large‐deviations (Wentzel–Kramers–Brillouin) approach, for the chance of invasion given the mean growth and its variance. The first formula has the virtue of simplicity, while the second one holds over a wider parameter range. The efficacy of the formulae, including their accompanying data analysis technique, is demonstrated using synthetic time series generated from canonical models and parameterised with empirical data.


I N T RODUC T ION
The capacity of natural systems to maintain diversity (e.g. genetic polymorphism and coexistence of multiple competing species) poses a major challenge to the theory of population genetics and community ecology (Chesson, 2000;Dempster, 1955). Interspecific interactions in diverse communities reflect many complex factors, and their effects are intertwined with strong fluctuations caused by environmental and demographic stochasticity (Kalyuzhny et al., 2014;Lande et al., 2003). A generic method that provides a simple metric to quantify persistence properties is therefore highly desirable. Turelli's (1978Turelli's ( , 1981 approach to this problem focuses on mutual invasibility, that is, on the 'conditions under which a rare invading species will tend to increase when faced with an array of resident competitors in a fluctuating environment'. This approach greatly simplifies the analysis of the complicated dynamics of a community by reducing it to a series of single-species invasion problems. All complex interactions are encapsulated in a few effective parameters that reflect the overall influence of the resident species and the environmental fluctuations on a given rare population. The rarity of the focal species allows one to neglect nonlinear (density dependent) effects, facilitating the analysis even further. An analogous approach is taken in permanence (uniform persistence) theories (Hofbauer, 1981;Hutson & Schmitt, 1992;Schreiber et al., 2011;Schuster et al., 1979).
Chesson's modern coexistence theory (Barabás et al., 2018;Chesson, 1982Chesson, , 1994Chesson, , 2000Chesson & Warner, 1981) bases its invasibility analysis on a single parameter, the mean growth rate of the invading species [r] (Chesson, 1994;Lewontin & Cohen, 1969;Turelli, 1978). The mean is taken over all the instantaneous growth rates r(t) where t denotes time. The same parameter, [r] , is employed in the adaptive dynamics theory of evolution (Brännström et al., 2013;Metz et al., 1995), in the study of epidemics and in many other fields [see Table 1 of Grainger et al. (2019)]. Consequently, ecologists have

Quantifying invasibility Jayant Pande
| Yehonatan Tsubery | Nadav M. Shnerb developed a collection of techniques to infer [r] from empirical time series or from models parameterised using empirical data (Grainger et al., 2019). Many contemporary studies of species coexistence and the maintenance of biodiversity analyse this 'invasion criterion' [r] (Chesson, 2008;Grainger et al., 2019) and partition it between underlying mechanisms such as niche differentiation, fitness differences and so on (Ellner et al., 2016(Ellner et al., , 2019Letten et al., 2018). Some authors employ [r]based criteria for mutual invasibility as a metric for stability in empirical studies (Usinowicz et al., 2012(Usinowicz et al., , 2017. However, recent work has revealed that [r] is not a reliable quantitative indicator of invasibility (Pande et al., 2020a). Systems with different underlying parameters may yield different probabilities of invasion for a species despite having the same numerical value of [r]. Even worse, invasibility may decrease even as [r] increases. Although [r] provides a fair binary classification-when its value is positive (negative), the extinction state is a repeller (attractor), from which important asymptotic predictions follow (Chesson, 1982(Chesson, , 1994Ellner et al., 2020;Hofbauer, 1981;Schreiber et al., 2011)-it does not reliably measure invasibility.
The failure of [r] as a metric reflects two problems (Ellner et al., 2020;Pande et al., 2020a;Pande et al., 2020b). First, since the growth rate r measures the logarithm of abundance ratios, zero-abundance (extinction) states lead to infinite negative contributions. To avoid this problem, [r]based analyses must neglect demographic stochasticity, the intrinsic stochasticity arising from the birth and death of discrete individuals, despite its crucial importance when the invading population is small. Second (Szilágyi & Meszéna, 2010), [r] does not fully include the effect of environmentally induced variation in abundance over time-that is, environmental stochasticity-which affects entire populations collectively.
Thus, given the limitations of [r], a better metric is required. Following former studies (Lande et al., 2003), Dean and Shnerb (2020) suggested such a metric for systems with environmental stochasticity. However, they did not take demographic stochasticity into account, and assumed that environmental fluctuations are weak such that the diffusion approximation (Crow & Kimura, 1970;Karlin & Taylor, 1981) is valid.
To provide a more general metric, here we present two new formulae that predict, quantitatively, the chance of invasion. The first of these is based on the diffusion approximation, differing from the work of Dean and Shnerb (2020) in its inclusion of demographic stochasticity. As is usual with the diffusion approximation, this formula is valid only for small levels of stochasticity. Our second formula is based on the large-deviations WKB (Wentzel-Kramers-Brillouin) approximation, and works for weak as well as strong levels of both demographic and environmental stochasticity. The WKB approach leads to more accurate results over a wide parameter range, and we recommend its use except in cases of very weak stochasticity, where the diffusion approximation formula may work better.
Moreover, our WKB formula converges to the known classical results for establishment probability (Gillespie, 2004;Haldane, 1927;Kimura, 1962) when the environment is fixed through time, and to the expression suggested in Dean and Shnerb (2020) when the diffusion approximation holds. Additionally, our analysis clarifies how the effect of demographic stochasticity may be taken into account indirectly by introducing an 'extinction threshold' at the right density. This allows one to extract quantitative predictions from infinite population models, like those used in permanence (uniform persistence) theories (Hofbauer, 1981;Hutson & Schmitt, 1992;Schreiber et al., 2011;Schuster et al., 1979).
Below we derive our two formulae, discuss their accuracy, explain how to extract the parameters required for their use from abundance time series, and finally compare their predictions to simulations of a few canonical models.

DE R I VAT ION OF I N VA SION FOR M U L A E The backward Kolmogorov equation and the diffusion approximation
In a community of N individuals, we define the chance of a species of n individuals to grow in abundance and reach a target abundance n f ≪ N as Π n→n f (or, in shorter notation, Π n ). Π n satisfies the following backward Kolmogorov equation, where W n→n+Δn are the transition probabilities from n to n + Δn individuals during a given period of time . This equation reflects the fact that the chance of a population of abundance n to ultimately reach a target abundance n f is equal to the chance to jump from n to any other destination, W n→n+Δn , multiplied by the chance to reach n f from this destination. Normalisation requires ∑ Δn W n→n+Δn = 1 , which means Π n = constant is always a solution of Equation (1). The lower bound Δn = − n corresponds to extinction and the upper bound Δn = N − n to ultimate fixation. The boundary conditions are Π 0 = 0 (reflecting no chance of invasion from the extinction state) and Π n f = 1 (because n ≥ n f implies successful invasion).
In Equation (1), W n→n+Δn are time-averaged transition probabilities, with the average taken over all the states of the environment. This averaging assumes the time period , unspecified so far, to be larger than or equal to the dwell time of the environment (the typical duration for which the environment stays constant before changing), otherwise the transition probabilities are correlated. From here on we assume to be the dwell time. We also assume that the typical time to invasion is much larger than ; in the opposite case the environment varies only slightly over the time required for the population to reach n f (Cvijović et al., 2015;Mustonen & Lässig, 2008) and the chance of invasion may be found by averaging over fixed environments, as explained in the discussion section.
The diffusion approximation (Karlin & Taylor, 1981) replaces the complicated difference Equation (1) by the following simpler differential equation that regards Π n as a continuous function of n (writing it as Π(n); throughout this paper, discrete arguments of functions are written as subscripts and continuous ones are put in parentheses), Here primes represent derivatives with respect to n. In the invasion regime, the invading population is rare, so we can neglect density-dependent effects: different individuals of the invading species do not interact with each other, and the reproductive success of an invader, its chance to produce m offspring, is dictated solely by the environment (including abiotic factors and the effect of the resident species). The total growth in the entire invading population is thus n times [m], the expected net growth associated with a single individual, that is, [Δn] = [m]n. The variance of this quantity is made up of two terms: demographic stochasticity, which is proportional to n and appears even if the environment is fixed, and environmental stochasticity, which is proportional to n 2 and appears only when the environment fluctuates (Lande et al., 2003). We employ ℰ,  d and  e as the relevant constants of proportionality, defined by Putting the forms of [Δn] and Var [Δn] from Equation (3) in Equation (2), the two independent solutions of the differential equation are found to be.
where The boundary conditions Π(0) = 0 and Π n f = 1 must be satisfied by a linear combination of Π I (n) and Π II (n), and this yields (see Supplement S1) the invasion formula under the diffusion approximation (DA),

Beyond the diffusion approximation: A large-deviations (WKB) formula
The diffusion approximation is based on two assumptions. It requires Δn to be small with respect to n (so Π(n + Δn) is approximated by terms up to (Δn) 2 in a Taylor series), and it assumes that the mean displacement in each step is negligible with respect to its standard deviation (technically, nℰ 2 ≪  d + n e ) (Karlin & Taylor, 1981). When abundance fluctuations and/or ℰ are large, as occurs in many biological scenarios, the diffusion approximation may fail. To overcome this difficulty, it is necessary to take recourse to a largedeviations technique. In recent years, such techniques, based often on the WKB (Wentzel-Kramers-Brillouin) approximation, have been used in several studies (Assaf & Meerson, 2017;Kamenev et al., 2008;Kessler & Shnerb, 2007;Ovaskainen & Meerson, 2010). Here we employ the WKB technique with the two-destination approximation (Pande & Shnerb, 2020) scheme that facilitates the analysis and allows us to derive the required invasion formula.
In this analysis, it is more convenient to use the logarithmic parameters z ≡ ln n and Δz ≡ ln(n + Δn) − ln n. On the logarithmic axis, Equation (1) becomes where the bounds of the sum are simply those of Equation (1) transformed to the zspace. Below we switch between the nand znotations or even mix them as convenient, bearing in mind the reciprocal relations z = ln n and n = exp(z).
To use the WKB technique, it is helpful to simplify Equation (7) while retaining its character as a difference equation. Our strategy for doing this employs the two-destination approximation (Pande & Shnerb, 2020), wherein we replace the actual random process in Equation (7) by a simpler process that allows, for each initial location z, only two values for the destination z + Δz. Despite this drastic reduction in the process complexity, no significant errors arise in the dynamical quantities of interest, as long as the first two moments of each jump are preserved (Pande & Shnerb, 2020). The first moment is the mean jump size [Δz] during the dwell time , and the second moment is (Δz) 2 . Note that [Δz] is related to the invasion growth rate by the relation [r] = [Δz]∕ (Chesson, 1982).
Here we stipulate that the two possible destinations lie on equal distances to the right and left of the initial position z, though the probability of jumping to the right and to the left may differ (see Supplement S2 for alternative formulations). Under this scheme, which treats Π as a continuous function of z, Equation (7) yields the simpler equation (4) where is the jump size in zspace and is the forward jump probability. and are in general functions of z, and must satisfy the following relations to leave the first two moments unchanged, Here the mean jump in the two-destination process is found by summing the product of the forward jump ( ) and the forward jump probability ( ) with the product of the backward jump (− ) and the backward jump probability (1 − ), and similarly for the second moment of the jump. Equations (9) imply To solve for Π(z) in this two-destination process, we apply our WKB-based approach. The basic idea is simple: the diffusion approximation requires Π(z) to be smooth over z, so one may expand Π(z + Δz) in a Taylor series and keep only the lowest orders. The large-deviations technique is based on a less restrictive assumption, that the logarithm of Π, that is, S(z) ≡ ln Π(z), is smooth over z. Therefore, we express Π(z) as e S(z) , and expand S(z + Δz) as S(z + Δz) ≈ S(z) + ΔzS � (z), where a prime denotes a derivative with respect to z. For brevity of notation, we define q ≡ S � (z). Equation (8) then gives where sinh is the hyperbolic sine function. This equation has two solutions. The trivial one, q = 0, yields the trivial constant solution Π I . The nontrivial solution for q, given by determines the nontrivial solution Π II (z) = exp[S(z)] = exp ∫ q(z) dz . A linear combination of these two solutions must satisfy the appropriate boundary conditions.
Since and are both z dependent, it is difficult to integrate q(z) over z. To obtain an analytic expression, we therefore interpolate between two limits: the large-n (outer) region where demographic stochasticity is negligible and q is constant, and the small-n (inner) region where environmental stochasticity is negligible.
The outer regime: We have everywhere assumed the invading population to be small (n ≪ N) such that there are no density-dependent effects and the only influence of n on the displacement Δz is through demographic stochasticity. Therefore, if n is large enough such that demographic stochasticity is negligible-this condition defining the outer regime-Δz is z-independent and we can define the Here q is obtained by substituting the expressions for and from Equation (10) in Equation (12). The inner regime: If the invading population is very small, birth-death events (among the invaders) are rare, and the environment changes a few times between two successive such events. Since the birth of different individuals in a species takes place in different environments, there are no species-wide, collective responses to the environment. Accordingly, although environmental variations can still affect the mean growth in the abundance (such as through the storage mechanism; Chesson & Warner, 1981;Chesson, 2000), their effect on random fluctuations in the abundance is inherently demographic in nature. The net effect is that in the inner regime the dynamics are those of a population in a constant environment, characterised by the mean growth rate E 0 and the strength of demographic stochasticity V d . For such a system, with an effectively constant environment, the diffusion approximation is known to be highly accurate (Parsons et al., 2010). In the nspace, the demographic stochasticity enters Equation (2)  where E 0 and V d are zindependent constants. Again, there are two solutions to Equation (15), the trivial solution Π I (z) = constant and the nontrivial solution Because z = ln n, integration over dz equals integration over dn ∕ n, therefore Interpolation: Equations (13) and (16) give Π II (z) when n is large and small, respectively. To interpolate between these two solutions, we interpolate between their corresponding values of q, using Clearly, q inter (n) approaches q in at small n and q at large n values, as desired. Therefore, where Combining this solution with the trivial (constant) solution so that the boundary conditions (Π n=0 = 0 and Π n=n f = 1) are satisfied, we obtain the final WKB-based invasion formula, with R and q defined in Equations (19) and (13).
Importantly, when nR ≫ 1 (which also means n f R ≫ 1 , as n ≤ n f ), the chance of invasion tends to 1 if q < 0 and to n∕n f q if q > 0. Since sign(q) = − sign E 0 , when the invading population is large (n ≫ 1∕R) and has a fitness advantage (E 0 > 0) it is certain to establish. If n f ≫ n, then the opposite statement also holds: an invading population with negative E 0 cannot establish. Therefore, only when 1 ∕R ≪ n ≪ n f does invasibility become a binary property for which the mean invasion growth rate [r] is a proper metric (Ellner et al., 2020;Pande et al., 2020b).

I N F ER R I NG T H E PA R A M ET ER S F ROM A BU N DA NC E T I M E SER I E S
Our invasion formulae, Equations (6) and (20), each require three parameters, reflecting the mean and the variance of abundance variations during the dwell time of the environment. Here we explain how to infer these parameters and the value of the dwell time from abundance time series. We start by considering cases in which artificial time series can be generated using simulations of a community dynamics model parameterised using empirical data. The length and the accuracy of such time series are, in principle, unlimited. Later, we discuss the case of empirical time series, which may be short and plagued by observation errors and thus difficult to accurately estimate the relevant parameters from.
The starting point of our analysis is a time series of the abundance n t (Figure 1a). To ascertain the parameters required for the diffusion approximation formula, the mean and the variance of Δn (the change in n during the duration ) are measured, as in panels (b) and (c), and from these ℰ,  d and  e are determined as illustrated in the insets in these panels. Likewise, to find the parameters for the WKB-based formula, the abundance series n t is first transformed to a series for z t = ln n t [panel (d)], from which the mean and the variance of Δz are found [panels (e) and (f)], and the parameters E 0 , V d and V e determined [insets in panels (e) and (f)].
Importantly, demographic stochasticity has negligible effect on our ability to measure ℰ. Demographically induced jumps are symmetric in the nspace, so they leave the value of [Δn] unchanged. In contrast, demographic stochasticity does diminish the value of [Δz] (because a symmetric, zero-mean change in n implies a net negative change in ln n). This effect becomes more pronounced as n decreases, so [Δz] becomes negative at small values of n [see Equation (15)]. Moreover, because the logarithm of zero diverges, transitions to the state of zero abundance must be excluded from the time series, which generates an artificial bias towards invasion when n is very small. For these two reasons, if the available abundance time series has only a small number of invaders such that demographic stochasticity plays an important role in the dynamics, then the inference of the parameters ℰ,  d and  e from the nspace analysis is more accurate than the inference of E 0 , V d and V e from the zspace. In such cases, the diffusion approximation formula Equation (6) may work better than the WKB-based Equation (20).
On the other hand, once the relevant parameters are known, the WKB-based formula performs better than the diffusion approximation-based one. When the constant selection terms (either ℰ or E 0 ) and the environmental stochasticity terms ( e or V e ) are small, both formulae estimate the chance of invasion well, but in cases of strong selection or large fluctuations the diffusion approximation becomes inaccurate. The superiority of our WKB formula is demonstrated in Figure 2.
The preceding discussion assumes the existence of a long and accurate time series, which is rare when the relevant datasets are empirical. We note nevertheless that the inference of the main parameter used in modern coexistence theory, [r] = E 0 ∕ , encounters the same challenges as the inference of the three parameters required for our formulae (see Supplement S3). Moreover, it is often possible to estimate the typical number of propagules per individual, say, or of hatchlings that reach maturity in a nest during a good or a bad year. From such estimates, approximate values for demographic and The raw data is an abundance time series (as obtained from direct observations or a numerical simulation of a properly parameterised model). Here this dataset of n t , n t+ , n t+2 , … is illustrated in panel (a). Each pair of adjacent points yields a single value Δn = n t+ − n t . The mean value [Δn] is plotted against n in panel (b , we filtered from the dataset all the cases where n t+ = 0; this yields an artificial positive bias when n is very small. Finally, the variance of Δz is shown (in panel (f)) to increase as n decreases, and the inset shows that it satisfies Var[Δz] = V e + V d ∕ n. All the datasets were obtained for the individual-based version of the lottery model (Chesson, 1982;Chesson & Warner, 1981) with = 0.2, s 0 = 0 and e = 0.25 (see Supplement S4 for model details). environmental stochasticity may be determined (as in Saether et al. (2004), for example), and then our formula may be used to obtain an order of magnitude estimate for the chance of invasion. Finally, a note on , the dwell time of the environment. In our model, the environmental conditions are fixed on average for periods of duration and then take an arbitrary and uncorrelated new value. This means that the correlation time of the environment is the mean time before the next shift, that is, ∕2. Therefore, equals twice the value of the autocorrelation time of the environment.

T H E I N VA SION FOR M U L A I N AC T ION
The results shown in Figure 2 (and, more expansively, in Supplement S5) indicate that our WKB-based formula estimates Π n well, and does so for a much wider range of parameters than our diffusion approximation formula. Henceforth, we treat the WKB-based formula as the recommended one. In this section, we examine its accuracy using a few model systems, and contrast it with quantitative predictions based solely on the invasion growth rate.
We use simulated data from three individual-based versions of canonical models: the discrete-time lottery model of Chesson and Warner (1981), its continuous-time (Moran) analogue, and the Leslie-Gower competition model for trees and saplings employed by Usinowicz et al. (2012) and Usinowicz et al. (2017). Detailed descriptions of these models and of our simulations are provided in Supplement S4. In all these models, there is a negative covariance between environment and competition, and thus environmental stochasticity may promote invasion via the storage mechanism (Chesson, 2000;Chesson & Warner, 1981). This makes the relationship between the underlying process parameters and the invasibility properties intricate and presents a demanding test to our formula. collapse on to a line also means that in this case one cannot choose a trajectory along which the predicted and observed values of the invasibility have different trends.  .1 and 1). In the left panel, the true chance of invasion is plotted against [r], demonstrating that [r] is a poor metric: For a given value of [r], the true invasibility may have many different values, and when the model parameters and change (e.g. along the path indicated by the grey arrow), the invasibility may even decrease while [r] increases. In the right panel, the true chance of invasion is plotted against the predictions of Equation (20), showing not only a data collapse but also close agreement with the numerical values. All the results were obtained for the case where there is no mean fitness advantage of the invading species, that is, s 0 = 0 ( [r] is positive even in this case due to the storage effect), and with N = 10 5 . Similar graphs for n f = 500 and for n f = 20 are presented and discussed in Supplement S6, with the latter case providing a good example where Equation (6), based on the diffusion approximation, works better than Equation (20) based on WKB, due to higher accuracy of parameter extraction. to establish, with n varying. We allow the mean fitness of the invading species to differ from the mean fitness of the resident species, a property quantified by the parameter s 0 . Again, [r] performs poorly as a quantitative predictor for the chance of invasion, whereas Π WKB n→n f predicts the true invasibility well. The limitations in [r] are more prominent here, since [r] is by definition oblivious to n , so all points for a given value of s 0 have the same value of [r] irrespective of n, despite being widely separated in their values of the observed invasibility.
In Figure 5, we show results for the individual-based version of the Leslie-Gower model for the dynamics of trees and saplings. We have used the recruitment rates for the species Spondias mombin and Spondias radlkoferi (scaled by a constant factor), as provided in Usinowicz et al. (2012). While in the models used for Figures 3 and 4 the dwell time is a model parameter, here it needs to be identified as well from the abundance time series (just like the parameters E 0 , V d and V e ). As a run-through of our whole method, we have provided a full step-by-step analysis for this model, including the generation of the abundance time series and the identification of all the required parameters, in Supplement S8.
For this model, the two panels in Figure 5

DI SC US SION
The chance of an invading population to grow in abundance to a large size governs many ecological and evolutionary processes, and is one of the main metrics used to assess the stability of an ecological community. In this paper, we have provided formulae that quantify this F I G U R E 5 Chance of invasion in the Leslie-Gower forest dynamics model of Usinowicz et al. (2012) and Usinowicz et al. (2017), as described in Supplement S4. For a community size of N = 10000, we calculated the chance to reach n f = 1000, starting from varying n, and plotted the results against In our simulations, we employed the empirical recruitment rates for Spondias mombin and Spondias radlkoferi, as provided in Usinowicz et al. (2012), scaled by a constant factor. A step-by-step description of the analysis is presented in Supplement S8. chance in fluctuating environments in terms of the dwell time of the abundance fluctuations, the mean and the variance of these fluctuations, and the initial and target densities.
In the existing literature, the sign or the numerical value of [r] is used as a proxy for the chance of invasion (Ellner et al., 2016(Ellner et al., , 2019Letten et al., 2018;Usinowicz et al., 2012Usinowicz et al., , 2017. However, as noted by Grainger et al. (2019), 'demographic stochasticity can have important consequences for low-density invasion growth rates, but this is frequently overlooked … the impact of demographic stochasticity on outcomes predicted by ( [r]) remains largely an open question'. As shown above, these are two independent and important problems. First, the estimation of [r] from low-density time series yields reduced values for [r] and may even invert its sign (see Figure 1e). Second, even if [r] is known, it cannot be used as a quantitative metric unless the invading population is large enough (1 ∕R ≪ n ≪ n f ).
The formulae provided here, Equations (6) and (20), facilitate a more informed analysis. First, the invasibility may be calculated using these formulae for any size of the invading population, as our treatment takes demographic stochasticity into account. Second, the diffusion approximation formula Equation (6) depends on the parameters in the nspace, in particular, on ℰ = [Δn] ∕n, which is not affected by demographic stochasticity (and does not diverge upon transitions to the zero-abundance state, unlike [r]). This formula may thus be used even when the relevant datasets are limited to small abundances. The WKB-based formula Equation (20) may be employed when [Δz] is calculated in the regime where it is independent of the abundance.
The price one has to pay for using our formulae is the need to calculate two more quantities, V d and V e (or  d and  e ). As illustrated in Figure 1, the calculation of these additional parameters is of a similar level of complexity as that of [r] (see Supplement S3). In some cases, even when the time series are not long enough, it may be possible to estimate these parameters to within acceptable error through heuristic means, in which case ) is again a more informative metric of the chance of invasion than [r].
Given the simplicity and accuracy of our formulae for Π n , we suggest that they replace [r] as an invasion metric in ecological and evolutionary analysis. A larger Π n implies a larger chance for a given population to invade or to recover from disturbances. Therefore, all other things (rate of disturbances, rate of migration from a regional pool, etc.) being equal, an increase in the value of Π n implies more successful colonisation and lower extinction rates, unlike the case for [r] (Pande et al., 2020a).
As an invasion parameter Π n may be used for the same purposes for which [r] is currently employed. In particular, one may partition its numerical value to study the contribution of specific mechanisms, like covariance between the competition and the environment (Ellner et al., 2016). Moreover, our formulae allow a quantitative estimate of the invasibility even when a precise numerical assessment through simulations is very difficult or computationally expensive, for instance, when n f is huge or when the chance of invasion is very close to 0 or to 1. This latter case is the situation in Usinowicz et al. (2017), where a comparison between the stability of different tree communities requires one to quantify the invasibility even though it is very close to 1.
Beyond their potential use in coexistence theory and related problems, our formulae provide some important insights.
First, they help to define a value of n f , the target abundance at which a species may be said to have established successfully. Two simple alternatives are choosing n f to be some fixed large number (independent of the overall community size) or a fixed frequency (a given fraction of the community), but these choices are somewhat arbitrary. Our WKB formula suggests that a natural definition for n f (when [r] is positive, so q is negative) is the abundance above which the chance of extinction is small, that is, n f ∼ exp(1∕|q|)∕R.
Second, because the incorporation of both demographic and environmental stochasticity in a model is a complicated task, many population models, like those used in permanence (uniform persistence) theories (Hofbauer, 1981;Hutson & Schmitt, 1992;Schreiber et al., 2011;Schuster et al., 1979), deal only with environmental variations, whose effect is generally stronger than that of demographic variations (Lande et al., 2003). However, without demographic stochasticity, the system never reaches the zero-abundance state, so in these theories one must introduce an arbitrary threshold below which a population is considered extinct. In Supplement S9, we show that our formulae converge to the expressions obtained in the literature when demographic stochasticity is neglected, if the extinction threshold is chosen at that is, at the abundance level at which the strengths of demographic and environmental stochasticity become equal. This result clarifies the range of parameters for which infinite population models are applicable: the total population has to be much larger than n th . Furthermore, by introducing an absorbing boundary at this threshold value n th , one may obtain quantitative approximations (for quantities like the mean time to extinction) from models that disregard demographic stochasticity.
Our invasion formulae have a few limitations. First, to derive them we assumed that the dwell time of the environmental variations is smaller than the time required for a population to successfully invade. When is larger than the typical invasion time (the "quenched" scenario of Mustonen & Lässig, 2008), the actual value of becomes insignificant, since the environment typically remains fixed throughout the invasion process (see (21) n th = V d ∕V e , Supplement S6). Our formulae do not apply in this regime, but they are not required either, since in this case one may use the well-known formula for the probability of invasion in a fixed environment (Haldane, 1927;Kimura, 1962) and determine the overall invasibility by taking an appropriately weighted expectation value over all possible environmental states. This fixed-environment invasion formula is, by construction, the limit V e → 0 of Equation (20) [and of  e → 0 of Equation (6)].
Second, as mentioned earlier, Equation (20) holds only if the parameters [Δz], V d and V e are constant in the invasion regime, which requires the ratio n f ∕N to be small, otherwise density-dependent effects appear.
Third, to use abundance time series, one must contend with sampling errors, which are common in ecological surveys (Clark & Bjørnstad, 2004). If these errors are unbiased, they do not affect [r] (or [Δz]), but they can overestimate V e . In some cases, these errors may be filtered out using standard techniques (Clark & Bjørnstad, 2004;Kalyuzhny et al., 2014). Otherwise, an informed guess about their effect on V e is required.
Species richness and genetic polymorphism reflect the balance between the rate of extinction and the rates (of colonisation, mutation, speciation, etc.) at which new types get established in the community. Invasibility, the chance of a population to invade or to recover from low densities, plays a crucial role in both kinds of processes. Consequently, our formulae should facilitate a more accurate and reliable assessment of these important aspects of life science systems.

AU T HOR S' CON T R I BU T ION S
J.P. and N.M.S. developed the formulae and analysed the data. Y.T. simulated the Leslie-Gower model for treesapling dynamics.