Permanence via invasion graphs: incorporating community assembly into modern coexistence theory

To understand the mechanisms underlying species coexistence, ecologists often study invasion growth rates of theoretical and data-driven models. These growth rates correspond to average per-capita growth rates of one species with respect to an ergodic measure supporting other species. In the ecological literature, coexistence often is equated with the invasion growth rates being positive. Intuitively, positive invasion growth rates ensure that species recover from being rare. To provide a mathematically rigorous framework for this approach, we prove theorems that answer two questions: (i) When do the signs of the invasion growth rates determine coexistence? (ii) When signs are sufficient, which invasion growth rates need to be positive? We focus on deterministic models and equate coexistence with permanence, i.e., a global attractor bounded away from extinction. For models satisfying certain technical assumptions, we introduce invasion graphs where vertices correspond to proper subsets of species (communities) supporting an ergodic measure and directed edges correspond to potential transitions between communities due to invasions by missing species. These directed edges are determined by the signs of invasion growth rates. When the invasion graph is acyclic (i.e. there is no sequence of invasions starting and ending at the same community), we show that permanence is determined by the signs of the invasion growth rates. In this case, permanence is characterized by the invasibility of all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-i$$\end{document}-i communities, i.e., communities without species i where all other missing species have negative invasion growth rates. To illustrate the applicability of the results, we show that dissipative Lotka-Volterra models generically satisfy our technical assumptions and computing their invasion graphs reduces to solving systems of linear equations. We also apply our results to models of competing species with pulsed resources or sharing a predator that exhibits switching behavior. Open problems for both deterministic and stochastic models are discussed. Our results highlight the importance of using concepts about community assembly to study coexistence.


Introduction
Understanding the mechanisms allowing interacting populations to co-occur underlies many questions in ecology, evolution, and epidemiology: When are species limited by a common predator able to coexist? What maintains genetic diversity within a species? Why do multiple pathogen strains persist in host populations? One widely used metric for understanding coexistence is invasion growth rates: the average per-capita growth rates of populations when rare. This approach has a long history going back to the work of MacArthur and Levins (1967), Roughgarden (1974), Chesson (1978) and Turelli (1978). This earlier work focused on the case of two competing species and lead to the mutual invasibility condition for coexistence. Namely, if each competitor has a positive invasion growth rate when the other species is at stationarity, then the competitors coexist. Mathematically, this form of coexistence corresponds to all species densities tending away from extinction. This corresponds to permanence or uniform persistence for deterministic models (Schuster et al. 1979;Sigmund and Schuster 1984;Butler et al. 1986;Garay 1989;Hofbauer and So 1989;Hutson and Schmitt 1992) and stochastic persistence for stochastic models (Chesson 1982;Chesson and Ellner 1989;Schreiber et al. 2011).
A key feature of the mutual invasibility criterion is that coexistence is determined by the signs of invasion growth rates. For many classes of multispecies models, positive invasion growth rates of at least one missing species from each subcommunity is a necessary condition for permanence to persist under small structural perturbations, i.e., robust permanence (Hutson and Schmitt 1992;Schreiber 2000). However, it need not be sufficient as in the case of three competing species exhibiting a rock-paperscissor dynamic (May and Leonard 1975). In this case, all single species equilibria can be invaded by a missing species but coexistence depends on quantitative information about the invasion growth rates at these equilibria (Hofbauer 1981;Hofbauer and Sigmund 1998;Schreiber 2000;Hofbauer and Schreiber 2010). This raises the question, when is it sufficient to know the sign structure of the invasion growth rates? Are rock-paper-scissor type dynamics the main barrier to qualitative conditions for permanence?
Invasion growth rates are the basis of what has become known as modern coexistence theory (MCT) or Chesson's coexistence theory (Chesson 1994;Letten et al. 2017; Barabás et al. 2018;Chesson 2018;Ellner et al. 2018;Grainger et al. 2019b, a;Godwin et al. 2020;Chesson 2020). To understand the mechanisms underlying coexistence, positive invasion growth rates are decomposed into biologically meaningful components and compared to the corresponding components for the resident species (Chesson 1994;Ellner et al. 2018). This decomposition allows ecologists to identify which coexistence mechanisms may or may not be operating in their system (Chesson 1994;Adler et al. 2010;Ellner et al. 2016Ellner et al. , 2018. For many applications of MCT, simple variants of the mutual invasibility condition determine which invasion growth rates need to be positive for coexistence (Chesson and Kuang 2008;Chesson 2018). For example, for two prey species sharing a common resource and a common predator, coexistence is determined by the invasion growth rates of each prey species into the three species community determined by its absence (Chesson and Kuang 2008). However, for more complex models, it is less clear how invasion growth rates determine coexistence (Barabás et al. 2018;Chesson 2018).
Here, we address these issues by introducing a mathematically precise notion of the invasion graph that describes all potential transitions between subcommunities via invasions (Sect. 2). For our models, we focus on a specific set of generalized ecological equations (Patel and Schreiber 2018). However, the proofs for our main results should hold for more general classes of models (e.g. reaction-diffusion equations (Zhao 2003), discrete-time models (Garay and Hofbauer 2003;Roth et al. 2017)) where conditions for permanence are determined by Morse decompositions (Conley 1978) and invasion growth rates. Under suitable assumptions about the ecological dynamics described in Sect. 2.2, we prove that if the invasion graph is acyclic and each subcommunity is invadable, then the species coexist in the sense of robust permanence (Sect. 3). In fact, we show that invasibility only needs to be checked at −i communities, i.e., communities without species i that are uninvadable by the remaining missing species. This sufficient condition always is a necessary condition for robust permanence. We show that our assumptions in Sect. 2.2 hold generically for Lotka-Volterra systems and provide a simple algorithm for computing invasion graphs (Sect. 4). We also apply our results to models of two prey sharing a switching predator and a periodically-forced chemostat with three competing species (Sect. 4). We conclude with a discussion about open problems and future challenges (Sect. 5).

Models, assumptions, and permanence
To cover ecological models accounting for species interactions, population structure (e.g. spatial, age, or genotypic), and auxiliary (e.g. seasonal forcing or abiotic variables), we consider a class of ordinary differential equations introduced by Patel and Schreiber (2018). In these equations, there are n interacting species with densities x = (x 1 , x 2 , . . . , x n ) taking values in the non-negative cone of R n . In addition, there are auxiliary variables, y = (y 1 , y 2 , . . . , y m ), taking values in a compact subset Y of R m . These auxiliary variables may describe internal feedbacks within species (e.g. genetic or spatial structure) or external feedbacks (e.g. environmental forcing or abiotic feedback variables). Let z = (x, y) denote the state of the system. In this framework, the dynamic of species i is determined by its per-capita growth rate f i (z), while the dynamics of the auxiliary variables are determined by some multivariate function g(z) = (g 1 (z), . . . , g m (z)). Thus, the equations of motion are (1) The state-space for these dynamics is the non-negative orthant K = [0, ∞) n × Y . The boundary of this orthant, K 0 = {(x, y) ∈ K : i x i = 0}, corresponds to the extinction of one or more species. The interior of this orthant, K + = K \ K 0 , corresponds to all species being present in the system. For any initial condition z ∈ K at time t = 0, we let z.t denote the solution to (1) with initial condition z.0 = z. Our first two standing assumptions for the equations (1) are: x, y) and (x, y) → g(x, y) are locally Lipschitz and, consequently, there exist unique solutions z.t to (1) for any initial condition z = (x, y) ∈ K. A2: The system is dissipative: There exists a compact attractor We are interested in when the equations (1) are robustly permanent, i.e., species persist following large perturbations of their initial conditions z and small perturbations of the equations governing their dynamics (Hutson and Schmitt 1992;Schreiber 2000;Garay and Hofbauer 2003;Patel and Schreiber 2018) for t sufficiently large. Namely, for all initial conditions supporting all species, the species densities are eventually uniformly bounded away from the extinction set K 0 .
(1) is robustly permanent if it remains permanent under perturbations of f i and g that satisfy assumptions A1-A2. More precisely, given any compact neighborhood V of , there exists δ > 0 such that dx i dt = x ifi (x, y), dy dt =g(x, y) is permanent whenever ( f (z), g(z)) − (f (z),g(z)) ≤ δ for all z = (x, y) ∈ V and (f ,g) satisfy assumptions A1-A2 with a global attractor˜ contained in V .

Invasion growth rates, schemes, and graphs
To understand whether species coexist in the sense of permanence, we have to consider species per-capita growth rates when rare, i.e., invasion growth rates. These are best described using ergodic probability measures that correspond to indecomposible dynamical behaviors of the model. Recall, a Borel probability measure μ on K is invari- for any continuous function h : K → R and any time t. Namely, the expected value of an "observable" h does not change in time when the initial condition is chosen randomly with respect to μ. An invariant probability measure μ is ergodic if it can not be written as a non-trivial convex combination of two invariant probability measures, i.e., if μ = αμ 1 + (1 − α)μ 2 for two distinct invariant measures μ 1 , μ 2 , then α = 1 or α = 0. The simplest example of an ergodic probability measure is a Dirac measure μ = δ z * associated with an equilibrium z * of (1). This Dirac measure is characterized by h(z)μ(dz) = h(z * ) for every continuous function h : K → R. Alternatively, if z * .t is a periodic solution with period T , then the measure μ defined by averaging along this periodic orbit is an ergodic measure (Mañé 1983;Schreiber 2000). Specifically, h(z) for all continuous h : K → R. More generally, the ergodic theorem implies that for every ergodic measure μ there exists an initial condition z * such that μ is determined by averaging along the orbit of z * , i.e., h(z)μ(dz) = lim T →∞ For any subset of species S ⊂ [n] = {1, 2, . . . , n}, we define to be the open face of K supporting the species in S. For an ergodic measure μ, we define the species support S(μ) ⊂ [n] of μ to be the smallest subset of [n] such that μ(F(S(μ))) = 1.
To understand whether a missing species i / ∈ S(μ) not supported by an ergodic measure μ can increase or not, we introduce the non-autonomous, linear differential equation to approximate the dynamics of species i's density x i when introduced at small densities. The solution of this linear differential equation satisfies Birkhoff's Ergodic Theorem implies that for μ almost every initial condition z = (x, y).
Consequently, we define the invasion growth rate of species i at μ as is also defined for the resident species i ∈ S(μ) supported by μ. In this case, we don't interpret r i (μ) as an invasion growth rate. Indeed, the following lemma shows that r i (μ) = 0 in this case, i.e., resident species have a zero invasion growth rate.
The proof of this lemma follows from the argument given for models without auxiliary variables y found in (Schreiber 2000, Lemma 5.1).
Proof Let i ∈ S(μ) be given. Let π i : K → R be the projection onto the i-th component of the x coordinate, i.e., π i (z) = x i when z = (x, y) ∈ K. Since μ(F(S(μ)) = 1, Birkhoff's Ergodic Theorem implies that there exists an invariant Borel set U ⊆ F(S(μ)) such that μ(U ) = 1 and Choose an open set V such that its closure V is contained in F(S(μ)), V is compact, and μ(V ∩ U ) > 0. By the Poincaré recurrence theorem, there exists z ∈ V ∩ U and an increasing sequence of real numbers t k ↑ ∞ such that z.t k ∈ V for all k ≥ 1. Since is V is compact, there exists a δ > 0 such that (2) and (3) imply that We make the following additional standing assumption: A3a: For each ergodic invariant Borel probability measure μ supported by 0 , r j (μ) = 0 for all j / ∈ S(μ), and A3b: sgn r j (μ) = sgn r j (ν) for any two ergodic measures μ, ν with S(μ) = S(ν), and all j.
Assumption A3a requires the invasion growth rates r i (μ) are non-zero for species not supported by μ. This assumption holds typically for dissipative Lotka-Volterra systems or systems with a finite number of ergodic measures. Due to their time averaging property, assumption A3b holds for all Lotka-Volterra systems and replicator equations (Hofbauer and Sigmund 1998) and certain types of periodically forced versions of these equations (Patel and Schreiber 2018). This assumption automatically holds when each face supports at most one invariant probability measure (e.g. there is a unique equilibrium, periodic orbit, or quasi-periodic motion in a given face). Sometimes this sign parity also can be verified when the per-capita growth functions f i exhibit the right convexity properties (e.g. Kon 2004;Schreiber 2004). For non-Lotka Volterra systems, it is possible for the per-capita growth rates of a missing species to have opposite signs at different ergodic measures. In this case, assumption A3b fails. For example, this failure arises in models of two predator species competing for a single prey species (McGehee and Armstrong 1977). If one predator has a type II functional response, then the predator-prey subsystem may simultaneously have an unstable equilibrium (defining one ergodic measure) and a stable limit cycle (defining another ergodic measure). McGehee and Armstrong (1977) showed that the invasion growth rates of the other predator species may be positive at the stable limit cycle but negative at the unstable equilibrium. A similar phenomenon arises in models of two prey species sharing a common predator (Schreiber 2004).
In light of assumption A3, we can uniquely define for each subset S ⊂ [n] of species. Let S be the set of all subcommunities: all proper subsets S of [n] such that S = S(μ) for some ergodic measure μ. For S ∈ S, there are at most n − 1 species and, consequently, μ(K 0 ) = 1 for any ergodic measure μ supported by S. Furthermore, S isn't empty as it always contains the empty community ∅ ⊂ [n]. Let |S| be the number of elements in S, i.e., the number of subcommunities. Now, we introduce our two main definitions. We define the invasion scheme IS to be the table of the signs of invasion growth rates {(r i (S)) i∈ [n] : S ∈ S}. When viewed as a |S| × n matrix, the invasion scheme is the signed version of the characteristic matrix introduced in (Hofbauer 1994). The rows of this matrix correspond to the different subcommunities while the columns correspond to the different species. We define the invasion graph IG as the directed graph with vertex set S and a directed edge from S ∈ S to T ∈ S if The first condition implies that there are no self-loops in the invasion graph. The second condition implies that all the species in T missing from S can invade S. The third condition allows for the loss of species from S that are not in T and ensures that these lost species can not invade T . One can view the invasion graph as describing all potential transitions from one subcommunity S ∈ S to another subcommunity T ∈ S due to invasions of missing species.
Remark 1 If z ∈ K is such that its α-limit set lies in F(S) for a proper subset S ⊂ [n] and its ω-limit lies in F(T ) for another proper set T = S, then there is a directed edge from S to T . The proof follows from the arguments presented in Appendix A.
Remark 2 As we focus on determining whether or not [n] is permanent, the invasion graph IG doesn't include [n]. However, for visualization purposes, we include [n] in our plots whenever [n] is permanent (see, e.g., Fig. 1). Schreiber 2022 provides R code for computing and plotting these invasion graphs.

Main results
The following theorem partially answers our main question,"when is knowing only qualitative information of the invasion growth rates r i (μ) (namely, their sign) sufficient for determining robust permanence?" Recall, that a directed graph is acyclic if there is no path of directed edges starting and ending at the same vertex, i.e., there are no cycles.
Theorem 1 Assume that A1-A3 hold and IG is acyclic. Then (1) is robustly permanent if for each S ∈ S there is i such that r i (S) > 0 i.e., each subcommunity is invadable.  Fig. 1 Invasion graphs, −i communities, and sample simulations for two 5 species competitive Lotka-Volterra models. Vertex labels correspond to the species in the community. −i communities for which species i has a positive invasion growth rate are colored olive green, and otherwise gold. Lighter shaded vertices correspond to non-permanent communities, all others correspond to permanent communities. Thicker directed edges correspond to single species invasions. Green directed edges indicate transitions due to species i invading a −i community. Both models have acyclic invasion graphs, but only the model in the top panels allows for robust permanence. Sample simulations of the models are shown in the right hand panels. Parameter values in Appendix C Remark 3 If ( f , g) are twice continuously differentiable and there exists S ∈ S such that r i (S) < 0 for all i ∈ [n]\S, then Pesin's Stable Manifold Theorem (see, e.g., Pugh and Shub (1989)) implies there exists z ∈ K + such that ω(z) ⊂ K 0 and, consequently, (1) is not permanent.
The proof of Theorem 1 is given in Appendix A. The idea of the proof is as follows. The invasion graph being acyclic allows us to construct a Morse decomposition (Conley 1978) of the flow on the extinction set K 0 . Each component of the Morse decomposition corresponds to a subcommunity in the invasion graph. The invasibility conditions ensure that the stable set of each component of the Morse decomposition doesn't intersect the non-extinction set K\K 0 (Schreiber 2000;Garay and Hofbauer 2003;Patel and Schreiber 2018). Then one can apply classic results about permanence (Butler et al. 1986;Garay 1989;Hofbauer and So 1989).
We derive two useful corollaries from this theorem. First, we show that if the invasion graph is acyclic, then checking for robust permanence only requires checking the invasibility conditions for a special subset of subcommunities. Specifically, given In words, a −i community is a community S that doesn't include species i and that can not be invaded by any of the missing species j / ∈ S except possibly species i.
S is not a −i community). As, by assumption, r i (S) > 0 for the −i communities, applying Theorem 1 completes the proof.

Remark 4
In general, −i communities correspond to subcommunities S for which there is an initial condition z = (x, y) such that (i) x i = 0 and x j > 0 for j = i, and (ii) the ω-limit set of z intersects F(S). Indeed, if there is a community S ⊂ [n]\{i} and initial condition z satisfying (i) and (ii), then the proof of Lemma 2 in Appendix A implies S is a −i community. Conversely, if S is a −i community and the functions f i , g are twice continuously differentiable, then Pesin's Stable Manifold Theorem (see, e.g., Pugh and Shub (1989)) implies there exists an initial condition z satisfying (i) and (ii).
Our second corollary concerns average Lyapunov function condition for permanence due to Hofbauer (1981). This sufficient condition is H: There exist positive constants p 1 , . . . , p n such that i p i r i (μ) > 0 for all ergodic measures μ with S(μ) ∈ S.
One can ask "when does knowing only the signs of the r i (μ) ensure that condition H holds?" To answer this question, we say the invasion scheme IS is sequentially permanent if there is an ordering of the n species, say 1 , 2 , . . . , n , such that column i of IS i has only non-negative entries where IS 1 = IS and IS i for i ≥ 2 is defined by removing the rows of IS i−1 where species i−1 has a positive per-capita growth rate. As sequential permanence implies that the invasion graph is acyclic and max i r i (S) > 0 for all S ∈ S, we get the following result: Corollary 2 Assume that A1-A3 hold. If IS is sequentially permanent, then (1) is robustly permanent.
If IG is sequentially permanent, then the invasion schemes of (1) restricted to each of the communities { 1 }, { 1 , 2 }, . . . , [n] are also sequentially permanent. Hence, each of the communities in this sequence is also permanent. Such a sequential way to prove permanence has been used by Hofbauer et al. (2008, p. 877 ff). A simple example of a system which is not sequentially permanent but has an acyclic invasion graph is given in Sect. 4.2.
To see how sequential permanence relates to condition H, consider the special case where r i (μ) = r i (ν) for ergodic measures satisfying S(μ) = S(ν) e.g. a Lotka Volterra model or models where each face supports at most one ergodic measure.

Applications
To illustrate the use of our results, we first describe how to verify them for Lotka-Volterra systems and also illustrate their application to two non Lotka-Volterra systems: a model of competing prey sharing a switching predator, and a periodically-forced model of three competing species. For the first two applications, the conditions of Theorem 1 are evaluated analytically, while in third application, we verify the conditions numerically.

Lotka-Volterra systems
Consider the Lotka-Volterra equations where x is the vector of species densities and there are no y variables (see, however, below for several extensions involving auxiliary variables). Let A be the n × n matrix corresponding to the species interaction coefficients and b the n × 1 vector of intrinsic rates of growth. Then f (x) = Ax + b. Assume A and b are such that the system is dissipative (i.e. A2 holds). Hofbauer and Sigmund (1998, ch. 15.2) provide various algebraic conditions that ensure dissipativeness. Furthermore, assume that each face of the non-negative orthant has at most one internal equilibrium. Under these assumptions, the Lotka-Volterra system exhibits the time averaging property. Namely, if z = x is an initial condition such that the ω-limit set of x.t is contained in F(S) for some S ⊂ {1, . . . , n}, then Therefore, r i (μ) = j A i j x * j + b i for any ergodic measure μ supported by F(S). In particular, assumption A3b is satisfied.
These observations imply that computing the invasion scheme and graph involves three steps.
Step 1: Find the set E of all feasible equilibria with at least one missing species: for each proper subset S ⊂ [n], solve for x ∈ R n such that (Ax) i = −b i and x i > 0 for i ∈ S, and x j = 0 for j / ∈ S. By assumption, E is a finite set. The vertices S of the invasion graph are given by S such that F(S) ∩ E = ∅.
Step 3: Compute the invasion graph by checking the directed edge condition for each pair of subcommunities in S, i.e., there is a directed edge from S to T iff r j (S) > 0 for all j ∈ T \S and r j (S) < 0 for all j ∈ S\T .
Two examples of using this algorithm for different 5 species competitive communities are shown in Fig. 1. In the case of the community in the top panel, we also plot the vertex [n] and transitions to this community. For both examples, the invasion graph is acyclic. For the community in the top panel, there is a unique −i community for each species and species i has positive invasion growth rates at this community. Hence, Corollary 1 implies that this system is robustly permanent, see sample simulation in the upper right panel of Fig. 1. Three −i communities (i = 2, 3, 4) are co-dimension one and, consequently, species i invading these communities (green directed edges) results in all species coexisting. The other two −i communities (i = 1, 5) have more missing species. For example, the −1 community {2, 3, 4} also misses species 5. When species 1 invades this community, species 3 and 4 are displaced leading to the −5 community {1, 2}. Successive single species invasions by species 5, 3 (or 4), and then 4 (or 3) assemble the full community. For the community in the lower panel of Fig. 1, there are nine −i communities. For three of these −i communities, species i has negative invasion growth rates. Hence, the system isn't permanent. Two of these −i communities ({1, 2, 5} and {1, 3, 4}) correspond to permanent subsystems where all the missing species have negative invasion growth rates. Hence, these −i communities correspond to attractors for the full model dynamics. Moreover, each is a −i community for each of the missing species e.g. {1, 3, 4} is a −2 and −5 community. The third of these uninvadable −i communities ({1, 3, 4, 5}) is a non-permanent system due to the attractor on the boundary for the {1, 3, 4} community. This explains the directed edge from {1, 3, 4, 5} to {1, 3, 4}. The lower, right hand panel of Fig. 1 demonstrates that the dynamics approach a three species attractor corresponding to one of the −i communities.
Certain modifications of the classical Lotka-Volterra equations also satisfy the timeaveraging property. Hence, for these modifications, computing the invasion scheme and the invasion graph also reduces to solving systems of linear equations. For example, Patel and Schreiber (2018) showed that if the intrinsic rates of growth are driven by a uniquely ergodic process (e.g. periodic, quasi-periodic), then this reduction is possible. In this case, one uses auxiliary variables dy dt = f (y) that are uniquely ergodic and replace b with vector valued functions b(y).

Competitors sharing a switching predator
Theoretical and empirical studies have shown that predators can mediate coexistence between competing prey species (Paine 1966;Hutson and Vickers 1983;Kirlinger 1986;Schreiber 1997). For example, Hutson and Vickers (1983) and Schreiber (1997) showed that a generalist predator with a type I or II functional response can mediate coexistence when one prey excludes the other, but can not mediate coexistence when the prey are bistable, i.e., the single prey equilibria are stable in the absence of the predator. Here, we re-examine these conclusions by considering a modified Lotka-Volterra model accounting for predator switching (Kondoh 2003).
Let x 1 , x 2 be the densities of two prey competitors. Let x 3 be the density of a predator whose prey preference is determined by the relative densities of the two prey species. Specifically, if y is the fraction of predators actively searching for prey 1 and 1 − y is the fraction actively searching for prey 2, then we assume predators switch between prey at a rate proportional to the prey densities, i.e., dy dt = x 1 (1 − y) − x 2 y. For simplicity, we assume the two prey species have a common intrinsic rate of growth r , normalized intraspecific competition coefficients and a common interspecific competition coefficient α. Under these assumptions, the predator-prey dynamics are where a is the attack rate of the predator and d is the per-capita death rate of the predator. The state space is K = [0, ∞) 3 × [0, 1]. As we are interested in predator mediated coexistence, we assume that a > d to ensure the predator always persists. Under this assumption, the single prey-predator subsystem x i − x 3 − y with i = 1, 2 has a unique globally stable equilibrium given by x i = d/a, x 3 = r (1−d/a)/a and y = 1, 0 for i = 1, 2, respectively. If α ∈ [0, 1), then the x 1 − x 2 − y prey subsystem has a globally stable equilibrium x 1 = x 2 = 1/(1 + α) and y = 1/2. Alternatively, if α > 1, then the x 1 − x 2 − y prey subsystem is bistable with a saddle at x 1 = x 2 = 1/(1 + α) and y = 1/2. As the invasion graphs for α = 0 are acyclic, Theorem 1 implies that robust permanence occurs if and only if the three equilibria associated with the −i communities {1, 2}, {1, 3}, {2, 3} are invadable. Invasibility of {1, 2} and {1, 3} requires that r (1 − αd/a) > 0. This occurs whenever a/d > α, i.e., predation is sufficiently strong relative to interspecific competition. Invasibility of {2, 3} requires a/d > 1+α. In particular, unlike the case of non-switching predators (Hutson and Vickers 1983), predator-mediated coexistence is possible in the case of bistable prey. In this case, the invasion scheme is not sequentially permanent.

Three competing species in a periodically forced chemostat
Chemostat are used in laboratories to study the dynamics of interacting microbial populations. As they provide a highly controlled environment, they are the basis of many mathematical modeling studies (Smith and Waltman 1995). For example, when species compete for a single limiting resource with a constant inflow, Hsu et al. (1977) proved that (generically) the species with the lowest break-even point excludes all others. When the resource inflow, however, fluctuates, coexistence is possible. For example, using a mixture of numerics and analysis, Lenas and Pavlou (1995) and Wolkowicz and Zhao (1998) showed that three competing species can coexist when the inflow rate varies periodically. Specifically, if R denotes the density of the resource in the chemostat and x i the density of competitor i, then they considered a model of the form: where R 0 is the incoming resource concentration, D(t) = D 0 + a cos(ωt) is a periodically fluctuating dilution rate, and f i (R) = α i R β i +R corresponds to a type II functional response. We can put this model into our coordinate system by defining y 1 = R and (y 2 , y 3 ) to be points on the unit circle: dy 2 dt = −ωy 3 dy 3 dt = ωy 2 with y 2 2 + y 2 3 = 1 (6) where D(y 2 ) = D 0 + a y 2 . The state space for (6) is K = R 4 + × S 1 where S 1 ⊂ R 2 denotes the unit circle.
To better understand these effects of the amplitude of fluctuations on species coexistence, we numerically calculated the Lyapunov exponents for all subsystems and created the invasion graphs for different amplitudes of the dilution rate (Fig. 3). At the a=0.325 Fig. 3 Invasion graphs for three competing species in a periodically-forced chemostat for increasing values of the amplitude a of the dilution rate. Shaded nodes correspond to −i subcommunities; olive green shading corresponds to a positive invasion growth rate of species i and yellow shading a negative invasion growth rate. Parameters as in Fig. 2 amplitude value a = 0.3 used by Wolkowicz and Zhao (1998), we recover the invasion graph suggested by their analysis. Specifically there are only two -i communities, {1, 2}, {1, 3}, and both of these communities can be invaded by the missing species. Hence, Theorem 1 implies robust permanence. At a higher amplitude of a = 0.325, the community {1, 2} can no longer be invaded by species 3 and permanence no longer occurs, consistent with the loss of species 3 in Fig. 2B at a = 0.325. At a lower value of the amplitude, a = 0.275, the −2 community is not invadable, a prediction consistent with species 2 being excluded in Fig. 2B at a = 0.275. At an even lower value of the amplitude, a = 0.2, the invasion graph in Fig. 3 dramatically changes with the community determined by species 3 resisting invasion from the other two species, consistent with only species 3 persisting in Fig. 2 at a = 0.2.

Discussion
Modern coexistence theory (MCT) decomposes and compares invasion growth rates to identify mechanisms of coexistence (Chesson 1994(Chesson , 2000Letten et al. 2017;Chesson 2018;Barabás et al. 2018;Ellner et al. 2018;Grainger et al. 2019b, a;Godwin et al. 2020;Chesson 2020). Our work addresses two key question for this theory: When are signs of invasion growth rates sufficient to determine coexistence? When signs are sufficient, which positive invasion growth rates are critical for coexistence? To answer these questions, we introduced invasion schemes and graphs. The invasion scheme catalogs all invasion growth rates associated with every community missing at least one species. The invasion graph describes potential transitions between communities using the invasion growth rates. Potential transition from a community S to a community T occurs if (i) all the species in T but not in S have positive invasion growth rates when S is the resident community and (ii) all the species in S but not in T have negative invasion growth rates when T is the resident community. Our definition of invasion graphs is related to what is often called an assembly graph in the community assembly literature (Post and Pimm 1983;Law and Morton 1996;Morton et al. 1996b;Serván and Allesina 2021). For example, Serván andAllesina (2021, page 1030) define assembly graphs for Lotka-Volterra systems. Like our definition applied to Lotka-Volterra systems (see Sect. 4.1), vertices correspond to feasible equilibria of the model. Unlike our definition, Serván and Allesina (2021) only consider transitions between communities due to single species invasions. This more restrictive definition, however, may not exclude heteroclinic cycles between equilibria due to multiple species invasion attempts, i.e., the "1066 effect" of Lockwood et al. (1997). These heteroclinic cycles may exclude the possibility of determining permanence only based on the signs of the invasion growth rates (Hofbauer 1994).
We show that the signs of the invasion growth rates determine coexistence whenever the invasion graph is acyclic, i.e., there is no sequence of invasions starting and ending at the same community. For acyclic graphs, we identify a precise notion of what Chesson (1994) has called "−i communities", i.e., the communities determined in the absence of species i. Specifically, these are communities where (i) species i is missing, and (ii) all other missing species have a negative invasion growth rate. −i communities can be found, approximately (see Remark 4), by simulating initial conditions supporting all species but species i for a sufficiently long time, removing "atto-foxes" (Sari and Lobry 2015;Fowler 2021), and seeing what species are left. This characterization ensures that each species i has at least one −i community associated with it.
When the invasion graph is acyclic, we show that robust permanence occurs if, and only if, at each −i community, species i has a positive invasion growth rate. Thus, this result helps define the domain of modern coexistence theory which relies on the signs of invasion growth rates determining coexistence (MacArthur and Levins 1967;Chesson 1994Chesson , 2000Letten et al. 2017;Barabás et al. 2018;Ellner et al. 2018;Grainger et al. 2019b, a;Godwin et al. 2020;Chesson 2020).
Our work also highlights the importance of going beyond average Lyapunov functions when only using qualitative information about invasion growth rates. The average Lyapunov function condition for permanence requires the existence of positive weights p i such that i p i r i (μ) > 0 all ergodic measures μ supporting a strict subset of species (Hofbauer 1981). This sufficient condition for permanence has received more attention in the theoretical ecology literature (Law and Blackford 1992;Morton 1993, 1996;Chesson 2018Chesson , 2020 than sufficient topological conditions using Morse decompositions (Garay 1989;Hofbauer and So 1989), or conditions using invasion growth rates with Morse decompositions (Schreiber 2000;Garay and Hofbauer 2003;Schreiber 2004, 2010;Roth et al. 2017;Patel and Schreiber 2018). This is likely due to the more technical nature of these latter papers. However, when one only knows the signs of the invasion growth rates, the average Lyapunov condition only works for specific types of acyclic graphs. Specifically, those graphs corresponding to a nested sequence of permanent communities, {1}, {1, 2}, {1, 2, 3}, . . . , {1, 2, 3, . . . , n}, where species i + 1 has non-negative invasion growth rates for all communities including species 1, 2, . . . , i. While these special graphs arise in some situations (e.g. diffusive competition (Chesson 2018;Mierczyński and Schreiber 2002) or certain generalizations of mutual invasibility (Chesson and Kuang 2008)), many communities do not exhibit this special structure.
There remain many mathematical challenges for an invasion-based approach to coexistence. First, while our main assumption A3b naturally holds for Lotka-Volterra and replicator systems, it is too strong for many other systems. Notably, the assumption does not allow for the invariant sets supporting multiple ergodic measures at which the invasion growth rates for a species have opposite sign. What can be done in these cases is not clear as they can cause complex dynamical phenomena such as riddled basins of attraction (Alexander et al. 1992;) and open sets of models where permanence and attractors of extinction are intricately intermingled (Hofbauer and Schreiber 2004). More optimistically, for stochastic models accounting for environmental stochasticity, the story may be simpler. For these models, permanence corresponds to stochastic persistence -a statistical tendency of all species staying away from low densities (Chesson 1982;Benaïm et al. 2008;Schreiber et al. 2011;Benaïm 2018;Hening and Nguyen 2018;Benaïm and Schreiber 2019). Under certain natural irreducibility assumptions (Schreiber et al. 2011;Hening and Nguyen 2018;Hening et al. 2020), each face F(S) supports at most one ergodic measure; A3b naturally holds for these models. Using the stochastic analog of invasion growth rates, one can define invasion schemes and invasion graphs as we do here. For these models, it is natural to conjecture: if the invasion graph is acyclic and all −i communities are invadable, then the model is stochastically persistent.
Dealing with cyclic invasion graphs is another major mathematical challenge. When these cycles are sufficiently simple, their stability properties can be understood using either average Lyapunov functions or Poincaré return maps (Hofbauer 1994;Krupa and Melbourne 1995;Krupa 1997). For more complex heteroclinic cycles (even between equilibria), the path forward for characterizing coexistence via invasion growth rates is less clear (Hofbauer 1994;Brannath 1994). Even for cyclic graphs where invasion growth rates characterize coexistence, it remains unclear how to carry out the second step of modern coexistence theory, i.e., how best to decompose and compare invasion growth rates to identify the relative contributions of different coexistence mechanisms. We hope that future mathematical advances on these issues will be incorporated into a next version of the modern coexistence theory (MCT v2.1).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Proof We prove this lemma by induction on k.
Suppose the Lemma holds for k − 1, and all the M I for |I | < k are given. Then we define M I for |I | = k as the maximal compact invariant subset in F(I ). It exists (i.e. there are no invariant sets arbitrarily close to the boundary of F(I )) as the family of M J with J I forms a Morse decomposition of the boundary ∂F(I ) of F(I ), each M J is isolated, and hence ∂F(I ) is isolated. Hence properties 1. and 2. hold.
Finally, the assumption that IG is acyclic implies property 4 by choosing a suitable order on S k .
Taking k = n − 1 we get Lemma 4 If IG is acyclic then the family of invariant sets {M I : I ∈ S} is a Morse decomposition of 0 .
We also need the following lemma which follows from Assumption A3 and Schreiber (1998, Corollary 1) (see, also, Mañé (1983, Exercise I.8.5). Assume that for every I , there is j such that r j (I ) > 0. Then Lemmas 4 and 5 and (Patel and Schreiber 2018, Theorem 2) imply (1) is robustly permanent.

Appendix B. Proof of Proposition 1
Assume IS is sequentially permanent. Without loss of generality, we can assume the sequence is {1, 2, . . . , n}. Let C be a matrix with sgn(C) = IG. Let c + > 0 and −c − < 0, respectively, be the minimum of the positive and negative elements of C. By the definition of sequential permanence, C i ([n − 1]) = r i ([n − 1]) = 0 for i ≤ n − 1 and sgn(C n ([n − 1])) = r n ([n − 1]) = +1. Hence, for any choice of p 0, i p i C i ([n − 1]) > 0. Choose p n = 1. Next, for any S ∈ S such that {1, . . . , n − 2} ⊂ S, C i (S) = 0 for i ≤ n − 2, C n−1 (S) ≥ c + , and C n (S) ≥ −c − . Hence, for any choice of p = ( p 1 , . . . , p n−1 , 1) 0, i p i C i (S) ≥ p n−1 c + − c − . Choose p n−1 = 2c − /c + . Then, i p i r i (S) > 0 for any S ∈ S such that {1, . . . , n − 2} ⊂ S. Proceeding inductively (i.e. choosing p k ≥ 2( p k+1 + · · · + p n−1 + 1)c − /c + ), we get p = ( p 1 , . . . , p k ) 0 such that C p 0. Now suppose IG is not sequentially permanent. We begin by assuming this failure of sequential permanence occurs at the first step i.e. there is no species such that r i (S) ≥ 0 for all S ∈ S. Then for each species i, there is S i ∈ S such that r i (S i ) < 0. Let C be such that its positive entries equal 1/n, its negative entries equal −1, and sgn(C) = IS. Suppose, to the contrary, there exists p = ( p 1 , . . . , p n ) 0 such that C p 0. Then 0 < i p i C i (S j ) ≤ −p j + i = j p i /n for any j ∈ [n]. Adding these n inequalities leads to a contradiction. This completes the proof when sequential permanence fails at the first step. Now suppose that the definition fails at the k + 1 step with k < n i.e. (i) there exist distinct species 1 , 2 , . . . , k such that column i of IS i has only non-negative entries, and (ii) every column of IS k+1 has at least one negative entry. Then we can use the same argument restricted to IS k+1 as r i (S) = 0 for all i ≤ k and S ∈ S satisfying [k] ⊂ S. Fig. 1 Lotka-Volterra parameters for the coexistence panels in Fig