Alternative cliques of coexisting species in complex ecosystems

The possibility that some ecosystems can exist in alternative stable states has profound implications for ecosystem conservation and restoration. Current ecological theory on multistability mostly relies on few-species dynamical models, in which alternative states are intrinsically related to specific non-linear dynamics. Recent theoretical advances, however, have shown that multiple stable ‘cliques’—small subsets of coexisting species—can be present in species-rich models even under linear interactions. Yet, the mechanisms governing the appearence and characteristics of these cliques remain largely unexplored. In the present work, we investigate cliques in the generalized Lotka–Volterra model with mathematical and computational techniques. Our findings reveal that simple probabilistic and dynamical constraints can explain the appearence, properties and stability of cliques. Our work contributes to the understanding of alternative stable states in complex ecological communities.


Introduction
Multistability, the process by which a dynamical system can reach different stable states under the same environmental conditions, has key implications in ecology.The most dramatic one, albeit not completely understood, is that of 'catastrophic shifts': the phenomena by which ecosystems, from shallow lakes to microbial communities, can undergo unexpected and abrupt transitions between healthy and degraded states [1][2][3].Beyond examples of catastrophic shifts, the possibility of alternative stable states has long been known to play a role in ecological succession and community assembly.There, it provides a framework to understand historical community patterns and the contingency of species invasions [4,5].
Most mathematical theories explaining ecological multistability have been grounded on simple, few-species dynamical models [6,7].In these, alternative stable states can arise by imposing specific non-linear dynamics.A single dynamical equation can then stabilize at one of multiple stable abundances depending on the initial conditions [8,9].However, this oversimplified picture fails at explaining the mechanisms by which multistability can appear in complex communities where many different species interact [6,10].
Recent theoretical developments are starting to help us understand what multistability looks like in complex models [7,[11][12][13].More particularly, we now know that multiple stable states can appear in species-rich models without the need to impose specific non-linear dynamics [7,12].The most commonly used framework to investigate these systems is the generalized Lotka-Volterra (GLV) model, where the growth of each species is affected by self-regulation as well as positive and negative linear impacts by other species [14].Theoretical and empirical research on the species-rich GLV model has revealed a rich landscape of dynamical behaviors, from species coexistence to chaotic fluctuations and multistability (see e.g.[12,[14][15][16][17][18][19][20] and figure 1(A)).Depending on the strength and heterogeneity of interactions, the GLV model can exhibit four well-defined phases.These are (i) population divergence at high cooperation, (ii) a single state of species coexistence, (iii) a regime of competitive exclusion states, each with only one surviving species, and (iv) a much more complex and less understood phase, where chaotic fluctuations and multiple fixed points can arise (figure 1(B)).
In the present work we focus on the fourth phase.Recent developments have helped us understand the properties of chaotic dynamics in this phase [17,20].High heterogeneity in species interactions can drive the system towards a pinball-like scenario, where species approach a stable community state that is then invaded by other species from the pool [17].This rich behavior is potentially related to empirical evidence of chaotic ecological dynamics [21] and has been recently tested in microbial community experiments [22].
On the contrary, the properties of multiple stable states in this phase remain an uncharted territory [20].We know that these so-called clique states can be found across mathematical models and usually support small and excluding subgroups of coexisting species [7,23,24].Theoretical advancements have predicted the number and stability of cliques in the infinite species limit, which depend on the symmetry of the interaction matrix (see e.g.[16,18,24]).However, other key ecological properties of cliques remain to be explained.We do not totally understand the mechanisms by which these small cliques appear in the first place.Therefore, a fundamental question is uncovering the rules that determine which and how many species can coexist in a clique.
In the present work, we investigate cliques in the GLV model.Our goal is to uncover the mechanisms that dictate their presence and constrain their properties.Our results show that simple rules can predict how many species can coexist inside a clique, the patterns constraining their interactions and the impact of species migrations.Each of these processes depends predictably on the properties of the studied system, providing potentially testable signatures to detect cliques in experimental data.Our work provides a first step towards detecting the presence of alternative cliques of coexisting species in complex empirical systems.

Methods
In the present methods, we describe the GLV model (section 2.1), the rationale behind sampling interaction strengths at random (section 2.2) and how we parametrize them (section 2.3).The central aspect of the work is to systematically confront two mathematical approaches.First, we explain how we numerically simulate the dynamics of the model to find and characterize stable cliques (section 2.4).Second, we define a mathematical test to study if randomly selected subsets of species can be identified as a clique without simulating the dynamics (section 2.5).Finally, we describe how species migration and invasions are implemented (section 2.6).

The GLV model
The central point of our work is to study the GLV model of S interacting species [12].In this model, each species is labeled i = 1, 2, . . ., S and the abundance N i of species i follows In the absence of species interactions (a ij = 0 ∀i ̸ = j), the model does not show multistability.Our main goal is to analyze how interactions impact the presence of multiple stable states.A typical procedure is to set intra-species parameters r i and K i to have low heterogeneity, so that we can write r i = r and K i = K ∀i.We subsequently rescale time as t → rt, divide species abundances by carrying capacities (x i ≡ N i /K i ) and rescale interaction strengths (A ij ≡ a ij K j ) [12,20].This leaves us with a simpler set of differential equations to study: A fundamental question in theoretical community ecology is then to understand how the dynamical properties of this model depend on the properties of the interaction matrix A. In our notation, a positive interaction strength (A ij > 0) indicates cooperation and a negative strength (A ij < 0) indicates competition.While other researchers have used the opposite sign convention (see e.g.[12]), we believe that keeping the same sign for interaction strengths and their ecological meaning make the results easier to interpret.

Disordered systems and random interactions
The number of species (S) in natural ecosystems is often very large.Attempts at empirically characterizing the complete interaction matrix A become extremely complex.Inspired by the framework of disordered systems in statistical mechanics [25], another approach is to sample the elements of A from probability Its rich phase space can be visualized by measuring the fraction of initial conditions that end up in a stable attractor [12,20].In (A) we plot this fraction as a function of mean and standard deviation of interaction strengths A ij for systems of S = 50 species (see Methods).The GLV model has a rich phase space with four dynamical phases (divergence, coexistence, competitive exclusion and cliques+chaos).The circles mark two known intersections of the line separating the stable coexistence from the chaos and cliques regimes.In (B) we present an example of the complex dynamics of the cliques+chaos phase.We simulate the dynamics with µ = −0.5, σ = 0.5 and perturb population abundances every ∆t = 10 3 timesteps.Chaotic fluctuations, different stable community states and cycling behavior can be observed in a single system (see supplementary material section I.A for details on simulating panels (A) and (B)).
distributions.Since the influential work by Robert May [26], modeling species-rich communities with random interaction strengths has provided a unified lens to understand many ecological scenarios.When many species and interaction types become intertwined, the species-specific details cease to matter and community-level patterns emerge (see SI I.C and [26,27]).
In this context, parametrizing our model is simplified to tuning the statistical properties of A. Following the standard approach of the GLV model [12], we generate A ij from a normal distribution without imposing correlations between matrix elements [12,18].The interaction matrix A is built as an Erdős-Rényi graph, where A ij is normally distributed N (µ, σ) with probability p or A ij = 0 with probability 1 − p.Here we focus in the p = 1 scenario, so that σ is the only parameter that controls interaction heterogeneity, and discuss the impact of p < 1 in the supplementary material (SM) II.D.Each system is therefore characterized by S, a unique (µ, σ) pair and the interaction matrix A ∼ N (µ, σ).We investigate how the dynamics of (2) change with the mean µ and standard deviation σ of A (figure 1(A)).

Parametrizing (µ, σ) to find cliques
We first plot a complete phase space of (2) to visualize the parameter domain where cliques are found [12,20].We do so by setting S = 50, µ ∈ [−1.8, 0.5] and σ ∈ [0, 1.6], and generate a grid of 40 × 40 = 1600 systems.We integrate each system and count the fraction of simulations that reach a stable state.This allows us to have a first picture of the dynamical regimes of the model (figure 1(A)).The specific simulation procedure used to integrate each system is explained in the next section.
Our main task is to evaluate the properties of clique states.To explore the properties of cliques across a variety of systems, we set µ and σ values chosen from random uniform distributions falling inside the Cliques, cycles and chaos phase of figure 1(A).For simplicity, we choose the (µ, σ) area under study to be a rectangle that does not span the whole phase.This is done to simplify simulations and avoid contaminating results of clique properties with those of states in the neighbouring phases, as boundaries are not necessarily sharp for low S (figure 1(A), [12]).The choice of the shape might alter which specific cliques are found, but not their general properties (see e.g.figures S8 and S9 in [20]).
For sections 3.1 and 3.3, we sample µ and σ from uniform distributions (µ ∼ U(−1.3, 0.5), σ ∼ U(0.3, 0.6)) and study the properties of cliques appearing within this domain.We generate 1000 systems inside the rectangle, each with a unique (µ, σ) pair.Interaction strengths (µ) in this range are consistent with empirical estimates of GLV dynamics in bacterial communities with S ≈ 50 [22].Conversely, estimating empirical σ values for GLV-like dynamics remains an open task [28].
In section 3.2 we present the properties of cliques in a single system with µ = −0.5 and σ = 0.5.This is done to present our results in the simplest possible format for readability.The equivalent results for systems spanning the whole rectangle are presented in the SM (section II.C).
In section 3.2 we also explore the properties of cliques for systems with increasing S and close to the critical line separating stable coexistence from the cliques regime.This line begins at the point (µ, σ) = (−1, 0) and ends at (µ, σ) = (0, √ 2/S) [12,20,29] (figure 1(A)).The equation for the line is The steepness of this line then decreases with S. At higher S, the coexistence phase shrinks and the line lowers [12,20].We study the diversity of cliques as S grows by simulating systems with µ = −0.5 and displacing σ slightly above the critical line (σ = 0.5 √ 2/S + 10 −2 ).

Simulating the dynamical equations
Once a system is fully defined by µ, σ and S, we generate 100 different random initial conditions for the initial population abundances.For each initial condition, we integrate the system during ∆t 1 = 3 • 10 3 timesteps.Most clique states become stable much faster, but we choose a longer timeframe to allow for potentially long transients (see e.g.figure 1(B), where relaxation to clique states is of the order of ∆t ≈ 10 2 ).
To check if a state is stable, we further integrate the system during ∆t 2 = 10 2 additional timesteps and check if all species abundances have remained the same, as previously done in [7] (see SM I.A for additional details on the numerical implementation).
For those states that are stable, we define as x ′ * the vector containing only the abundances of surviving species, and A ′ is the reduced interaction matrix containing only the interaction strengths between these surviving species.A clique state is characterized by the properties of the system in which it appeared (S, µ, σ) and the diversity, abundances and interactions of the final surviving species (S ′ , x ′ * , (µ ′ , σ ′ )).This variables are the central object of study of the present work.
The interaction matrices of systems with S = 50 species already contain S(S − 1) = 2450 elements, so that the statistics of the finite-size sample A ∼ N (µ, σ) are likely to be very similar to those of the original distribution N (µ, σ).In systems with fewer species, smaller matrices might have statistics that deviate from those of the original distribution.Given that cliques support a small number of species S ′ < S, we try to find and explain the patterns underlying these deviations.First, we measure the difference in terms of mean and standard deviation between the initial system and each clique (∆µ = µ − µ ′ , ∆σ = σ − σ ′ ).Second, we study the presence of correlations between abundances x ′ * and interaction strengths A ′ (see [30] and Results section).Third, we measure how the statistics of interactions in a clique (µ ′ , σ ′ ) are related to the number of species S ′ .

An analytical framework to obtain cliques from A
We define a mathematical test to check if randomly selected subsets of species can be stable cliques of (2).Our framework is based on a two-step process: we first select a random set of S ′ < S species, and latter test if these selected species can coexist in a clique state.
Species that coexist in a clique are not expected to be a completely random subset of the original system.Yet, identifying the specific coexisting sets from (2) demands very advanced mathematical machinery [18].Instead, here we test four different methods to randomly select subsets of species from the original system.We either select species at random with probability γ (method 1) or proportional to their interactions (methods 2-4).For the latter method, we select a species with a probability proportional to the perceived impacts by other species (method 2), the exerted impacts towards other species (method 3) and the sum of both (method 4).For example, in the last method a species is chosen to participate in a potential clique with probability Here γ ∈ [0, 1] is constant across all species.By choosing larger or smaller values of γ we can control the average amount of species (S ′ ) that are selected.The max k term in the denominator indicates that, across all possible species k ∈ {1, . . ., S}, the filter selects the one for which the sum between parenthesis is largest.
For each ot the methods, we select S ′ species and obtain their interaction matrix A ′ from the original values of A. Now we set a test for the conditions by which it can be a feasible and stable state.First, the abundance of species can be calculated from x ′ * = (−A ′ ) −1 .Clique states will need to be feasible, so that only subsets with positive species abundances are considered (x ′ * i > 0 ∀i).Second, we compute the Jacobian matrix, with elements Clique states need to be stable, so that only subsets with negative eigenvalues of the Jacobian are considered.With this process, we filter random subsets of the original system that are feasible and stable.We next evaluate if the properties of these subsets are consistent with the properties of cliques found by numerically solving the system (see above).The details of the mathematical implementation can be found in the SM, section I.B.

Species migration
Albeit unstable, the species-extinction state x i = 0 is an absorbing state of the dynamics of (2).The relation between chaotic fluctuations and species extinctions has been long studied in theoretical ecology (see e.g.[31][32][33]) and has recently been analysed in the GLV model [17,20].Heterogeneous interactions in A can drive the GLV model to a state with chaotic abundance fluctuations that could drive species to extinction (figure 1(B)).In this section we study how the properties of cliques change if extinct species can reinvade.
In the previous sections, we have presented the conditions by which cliques are stable, assuming that extinct species cannot reinvade.However, if extinct species were able to reinvade, they could destabilize the stable community towards another fixed point or even a chaotic regime.This is commonly studied by introducing a minimal migration rate m ≳ 0: We implement this in the simulations (section 2.3) by admitting a very small population threshold below which the abundance of a species cannot keep decreasing (see [17] and SM I.A for details).In the mathematical tests (section 2.4), we study the effects of migration by asking if cliques are invadable.We test if at least one species not present in a chosen subset has a positive growth rate when invading at low abundance: Here, subscripts I and R indicate species belonging to the Invader (extinct) and Resident (clique) sets, respectively (see SM I.B.3).
In both the simulations and mathematical tests, we study the likelihood of finding uninvadable cliques as the size of the original system S grows and species can reinvade.In the simulations, by adding a small population threshold forbidding extinctions.In the mathematical tests, by testing if species not selected in a subset can invade it.The impact of migrations is analyzed because cliques are likely to become unstable in very large systems [18,20].The hypothesis we test is that the likelihood of finding a destabilizing invader increases rapidly as the size of the system S grows.

Results
The results of this paper are organized as follows.In section 3.1 and figure 2 we show that cliques display predictable patterns in species interactions.In section 3.2 and figure 3 we show how competition bounds the number of coexisting species in cliques.In section 3.3 and figure 4, we show how species migration and a larger system size can make cliques unstable.

Predictable interaction patterns in cliques
The GLV model harbors a phase with multiple stable states of coexisting species (figure 1(B)).We study how the interaction strengths in cliques (A ′ ) compare with the interactions in the original system (A).We start with strictly positive initial abundances (x i (t = 0) > 0 ∀i) and suppose that species which go extinct due to the dynamics are not able to reinvade (we explore the impact of invasions in section 3.3).
We first study if the interaction statistics of cliques (µ ′ , σ ′ ) deviate from those of their system of origin (µ, σ).We plot in figure 2(A) the difference in mean interaction (∆µ = µ − µ ′ ) and standard deviation (∆σ = σ − σ ′ ) between each original system and the stable cliques it generates.As a general rule of thumb, in the species-rich GLV model, we know that species will stably coexist if interactions are geared towards weak and homogeneous competition [12,26] (figure 1(A)).
We find that cliques systematically fall into a region where they are less competitive and less heterogeneous than in their system of origin A(µ, σ) (figure 2(A)).The 'less competitive, less heterogeneous' signature is found by numerically integrating (2) (see section 2.3).The same signature also appears by mathematically testing which randomly chosen species can coexist stably (see sections 2.4 and SM II.A).In the mathematical test, however, the condition is not as strict.This is because random species choices permit some competitive cliques with µ ′ ≲ µ that are hardly reached by the dynamics (see SM II.A for a detailed description and comparisons across methods).Cliques appear as subsets that harbor less competitive and less heterogeneous interactions than their system of origin.In (B) we test the theoretical feasibility condition from [30] see ( 6)).We find that the interaction strengths in cliques (A ′ ij ) correlate with species abundances (x ′ * i ), creating the predicted diffuse pattern.Roughly put, abundant species in a clique (larger circles in the community sketch) need to compete weakly (thinner interaction lines).Instead, less abundant species (smaller circles) can compete strongly (thicker lines) without impairing coexistence [30].Each dot represents an A ′ ij value and its mathematical expectation, while the dark line indicates the 1 : 1 relation expected from the theory.For both panels we plot in red the linear regression of the data to help visualize its qualitative trend.
We also study if there are deeper patterns at the level of species abundances and their impacts on other species.Barbier and colleagues described a species-level interaction pattern pervading coexistence in large communities [30].One of the mathematical results is that species that coexist at high abundance are found to compete weakly, while species with lower abundance are allowed to impact each other at higher strengths.This allows us to predict the expected value of an interaction between two species in the clique, A ′ ij , as a function of their abundances x ′ i , x ′ j (see [30] for the complete development): where x′ is a measure for the expected abundances that species would hold if all interactions had the same value A ′ ij = µ: Species-level interactions in cliques follow the same diffuse pattern (figure 2(B)).Interactions in a clique state (A ′ ) deviate from the randomness of A in a precise way, by which high-abundance species do not impact each other strongly.Species then coexist inside a clique following the same interaction patterns previously tested in large empirical communities [30].Cliques are then those interacting subsets (A ′ ) of the original system (A) that fulfill two constrains.First, they contain the most homogeneous and weakly competing interactions of the pool (figure 2(A)).Second, they ensure that abundant species in the clique are in turn weak competitors (figure 2(B)).Next we study what are the possible constraints behind these patterns.

The diversity of cliques is bounded by (µ, σ)
We have seen that cliques represent a subset of species of the original interaction matrix A that tend to hold weaker competition (µ ′ > µ) and lower heterogeneity (σ ′ < σ).We now study if there are explicit relations between the properties of cliques (S ′ , µ ′ , σ ′ ) and those of the system of origin (S, µ, σ).
We define a single system (S, µ, σ), and study the diversity of cliques that appear from it (S ′ ) and the lowest and highest interaction strength values these cliques support (min(µ ′ ), max(µ ′ )).We show that the relations between S ′ and µ ′ explain the maximum number of species that can survive in a clique given S, µ and σ.An equivalent analysis for a range of (µ, σ) values instead of a single one is shown in SM II.C. ).The maximum clique diversity scales as half of the total diversity, consistent with the diversity of the coexistence state before the critical transition [12].

Probabilistic constraint on highest µ ′
We first study the constraints of the random sampling process.Interactions in a clique, encoded in A ′ , are sampled from the original system A. As clique diversity S ′ increases and we sample more elements from A, the statistical properties of the clique (µ ′ , σ ′ ) will inevitably resemble those of the original system (µ, σ) [34].In figure 3(A) we plot in light gray the mean of randomly chosen subsets of A as their size increases.We see that a clique will have µ ′ → µ, σ ′ → σ as S ′ increases.The highest values of µ ′ that a clique can harbor follow the expected trend for the statistics of S ′ (S ′ − 1) elements sampled from a normal distribution [34].Because of the tendency to minimize competition (figure 2(A)), some cliques can fall as far as a 4σ distance from µ.Then, we can approximate that the clique with highest mean interaction strength is bounded by The most cooperative clique in the studied range can have a mean interaction as high as the original system permits (µ = −0.5),plus an added deviation proportional to the heterogeneity of the system and the number of species in the clique S ′ .The factor S ′ (S ′ − 1) is the number of elements sampled from A, as each of the S ′ species in the cliques interacts with S ′ − 1 other species.Equation ( 8) imposes a probabilistic constraint on µ ′ : cliques cannot be more cooperative than what the original system (µ, σ) and their size S ′ allow (dashed black line, figure 3(A)).The probability of finding cooperative cliques decreases very rapidly with their diversity S ′ .

Stability constraint on lowest µ ′
We next ask what is the smallest, most competitive µ ′ value that cliques can support.Observed µ ′ values (colored dots, figure 3(A)) do not appear as low as they could by sampling (light gray dots).The most competitive cliques observed (min(µ ′ )) are not as competitive as probability could allow (figure 3(A)).There is a stronger dynamical constraint limiting competition in stable cliques.
We know that there is a limit to how competitive and heterogeneous species can be to coexist stably [12,20,26,35].This (µ c , σ c )-limit induces a transition which separates stable community coexistence from the chaos and cliques regime (figure 1(A)).We hypothesize that species surviving in a stable clique fulfill the same coexistence condition: once other species have gone extinct, the remaining species must obey the GLV criteria for stable coexistence.
In the original parameter space, the condition is a line crossing from points [12,20,29]).The first point does not depend on system size: neutral coexistence (σ = 0) is always disrupted when competition overcomes self-regulation, no matter the number of species [15,36].The second point decreases with system size: a system with more species will have a smaller coexistence phase [12,26].For a given clique with mean competition µ ′ , heterogeneity σ ′ and number of species S ′ , the condition for stability is to ensure that average interaction strength is above A larger (S ′ ) or more heterogeneous (σ ′ ) clique requires more positive interactions to maintain stable coexistence [12,26].This imposes a lower bound to the competitive strengths that cliques can support given their size and interaction heterogeneity.
The mathematics behind (9) are valid for large systems: the critical transition is not sharp for S ′ ≲ 100, so that smaller systems could hold µ ′ ≲ min(µ ′ ) and still be stable.(see figure 5(b) in [12]).Because cliques are particularly small communities, it is likely that stable cliques of size S ′ can still exist beyond the threshold.
To estimate a correction of ( 9) for low S values, we perform a separate numerical test of destabilizing small GLV communities (see SM II.B for a complete description).In it, we generate systems with µ = 0 and find the value of σ 90% c for which 90% of states are unstable as S increases.As predicted [12], σ 90% c = √ 2/S beyond S = 50.However, at lower S values, we find that σ 90% c > σ: smaller communities can support higher interaction heterogeneity.We hypothesize that this is consistent with evidence of deviations in the most positive eigenvalue of a random finite-size matrix [37] (see SM II.B).Because cliques are particularly small, their stability threshold is better approximated by a line from 08 ), with γ = 14.118 (see SM II.B).Under this correction, the most competitive clique (min(µ ′ )) sampled from A that can be stable then sits at This value is lowest when the heterogeneity of clique interactions (σ ′ ) is also lowest.Given the original system with (µ, σ), the lowest standard deviation min(σ ′ ) that can be obtained when sampling S ′ (S ′ − 1) elements from A has an analytical expression depending only on σ and S ′ [34] (see SI II.B).Applying it, we show that the lowest interaction strength of cliques only depends on clique size (S ′ ) and the statistics of the original system: where σ is the standard deviation of the system of origin and n = S ′ (S ′ − 1).As S ′ increases, larger cliques need higher µ ′ values to maintain stable coexistence (figure 3(A), dashed red line).Conditions ( 8) and (11), corresponding to the black and red dashed lines in figure 3(A), explain the interaction constraints of cliques in the GLV model.Cliques cannot be more cooperative than what sampling from A allows, nor more competitive than what stable coexistence allows.

Competition bounds clique diversity
Our results explain the maximum diversity (S ′ M ) that clique states can attain given a system with fixed (µ, σ).This can be estimated by the crossing of the two constraints (black and red lines, figure 3(A)).Interestingly, a system with more species (S) but equal (µ, σ) will not allow for larger cliques: (S ′ M ) only depends on (µ, σ) (figure 3(B)).
Finally, we study the highest diversity that cliques can attain if (µ, σ) can vary.The highest diversity is found when µ is highest and σ is lowest, so that (µ, σ) is closest to the coexistence regime (figure 1(A)).Interestingly, the location of this transition depends on S. As S increases, the line decreases and the coexistence phase shrinks (figure 1(A)).Given a fixed µ, cliques appear at lower σ in larger systems (σ c = √ 2/S(1 + µ)).A clique generated from a system with lower σ is likely to hold more species (higher S ′ ).The prediction at the coexistence phase is that the number of surviving species decays as σ increases [12].At the critical line, the number of surviving species is about half of the system size.Right above the critical line, results indicate that cliques should maintain a similar diversity of S ′ = S/2 (see figure 3(a) in [12]).In figure 3(C) we plot the largest observed clique for a single system poised right above the critical line.We can see that the most diverse clique that each system harbors scales approximately as S/2.Smaller systems slightly overcome this threshold (see SM II.B), while in larger systems it becomes computationally harder to find the largest possible clique (figure 3(C)).

Species migration and clique instability
Theoretical results indicate that many clique states can exist in the GLV model (see e.g.[16,18] and [7] for a review of the literature).Yet, models estimate that most cliques are unstable in the many-species limit (S → ∞) [18], often leading to chaotic dynamics [20].In our analysis, however, we find many stable cliques.We plot the fraction of stable simulations against system size with (blue dots) and without (orange dots) species migration.We compare these with two theoretical predictions: the presence of successful invaders on cliques generated from the dynamics of the GLV model (dark line) and on cliques filtered by the mathematical procedure (gray line).Larger system sizes destabilize cliques due to adding invaders (blue dots, dark line).Randomly selected stable subsets do not accurately mimic the composition and resistance to invasions of cliques (gray line).
The main reason is that our studied system has a moderate number of species (S ∈ [10,100]) and no species migration (m = 0, ( 4)).If a species goes extinct, it remains so and cannot reinvade [33].Our hypothesis is that some of the observed cliques would become unstable if extinct species can reinvade.Moreover, if we increase system size, we would not observe larger cliques (see sections above).Instead, it would only increase the likelihood of finding potential invaders.In the many species limit, there would always be an invader to any clique state, reaching the predicted theroetical scenario of unstable cliques [16][17][18].
We test this hypothesis with two coherent methods (see SM I.B.3).First, we simulate (4) with m > 0 and count the fraction of initial conditions that reach a stable state.Second, we simulate (4) with m = 0 and test if there are species not present in each clique that fulfill the invasion criteria (5).Our results can be summarized into two main outcomes.
First, stable cliques become increasingly rare as system size increases and species do not get extinct (m > 0, figure 4, blue dots).The theoretical tests of successful invasions (m = 0, (5)) follow the same trend (figure 4, dark line).This is consistent with the hypothesis that larger systems harbor fewer stable cliques: system size S increases the likelihood of finding a destabilizing invader [16][17][18].On the contrary, a larger fraction of cliques is stable if species can go extinct or the system size is moderate (figure 4, yellow dots).
Second, the mathematical test provides a fast method to estimate aggregate clique properties (figures 2, 3 and SM II.A) but does not accurately capture the effects of invasions (figure 4).These errors become more likely as system size increases.This is because, in the mathematical tests, we select species at random and test if they can coexist stably without accounting for invasions by other non-selected species.The selected species set accurately recovers typical clique properties (see above and SM II.A), but it is likely that it is invadable when other species can migrate from the pool (figure 4, gray dashed line).To find uninvadable cliques, one would require a much slower test that checks each species subset first for stability and then against all potential invaders [38].

Discussion
Cliques have been identified in a variety of complex biological models [7], making their analysis relevant across life science disciplines.Previous research has uncovered the phase in which cliques emerge [12] or the number of cliques [18] in the species-rich GLV model and across comple biological systems [7].Here we have shown that simple rules can explain the basic ecological properties of these small cliques of coexisting species.

4.1.
How do cliques arise in random interaction models, and why are they so small?Multistability in community models appears when the interaction network (A) is too heterogeneous for species to coexist stably [12,26,35].Our results indicate that the specific subsets of stable coexisting species beyond this threshold fulfill three main constraints.
First, cliques are those subsets of species of the original pool that have the most positive and homogeneous interactions possible.This is consistent with the notion that strong competition [36] and interaction heterogeneity [26] destabilize community coexistence.At the species-level, coexistence also selects for a subtle non-random pattern, whereby pairs of species with high abundance are in turn weak competitors [30].
Second, the number of species that a clique can support is bounded by competition.Cliques with more species increase their average competition while at the same time reduce the competitive strength they can support.This two-fold constraint explains the small diversity of most clique states: too many coexisting species rapidly overcome the stability threshold.At most, clique states emerging very close to the stable coexistence regime can harbor about half of the species of the original system.
Third, species migration from the original pool can destabilize cliques.Larger systems cannot generate more diverse cliques, but in fact only provide destabilizing invaders.This explains mathematical predictions by which most cliques become unstable at very large community sizes [18].As S grows, the likelihood of finding a potential invader increases and stable cliques become rarer.

Empirical evidence of cliques
We outline potential links between the above results and field and laboratory experiments.In the field, community shifts have been observed across environmental and altitudinal gradients [39][40][41].Sharp shifts between different bacterial compositions have also been observed in experiments and clinical analyses of microbial communities [42][43][44].These two avenues hint towards the possibility that experiments with plant communities [45][46][47] or bacterial ecosystems [3,48,49] could hold key information about ecosystem multistability.
Confronting dynamical observations with theoretical predictions is a much more nuanced problem.Using the GLV model, recent work has observed the mathematical signatures of coexistence and chaotic dynamics in field and laboratory experiments [22,30].The GLV model has also been recently used to test multistability signatures in microbial communities at a much lower diversity of two to six bacterial strains [48,49].Connecting theoretical explanations of multistability and empirical species-rich observations is a promising avenue [50].
Towards this goal, our results provide a set of empirical fingerprints to assess the presence of cliques in complex empirical communities.These are (i) the selection for weaker and homogeneous competition (figure 2(A)), (ii) the correlation between abundances and competition strengths (figure 2(B)), (iii) how competition bounds species richness (figure 3(A)) and (iv) how community stability decays with species invasions (figure 4).

Limitations and next steps: the role of non-random interactions
Empirical parametrization of complex species-rich communities is particularly difficult [28].To overcome this barrier, we have followed the common convention of sampling interactions at random [26].This theoretical framework has recently showed its potential in capturing the behavior of experimental systems of many interacting species [22,30].
Ecological networks, however, typically hold non-random structures [51][52][53].Non-random interaction patterns have been suggested to play an essential role in the occurrence of multiple stable states [11,13], but the question of how the properties of the GLV model are be affected by non-random interactions remains open.We have shown that the properties of cliques are governed by the statistical properties of A. If the interaction matrices A were non-random, we hypothesize that the properties of cliques could change dramatically.

Conclusion
Multistability has been observed in natural species-rich communities, both in the form of large catastrophic transitions as well as sharp shifts between similar community states.Because of the importance of these shifts, an important body of literature has been devoted to understand the mechanisms underlying multistability.While simple dynamical models can provide an explanation for multistable behavior, the mechanisms explaining multistability in complex ecosystems remain mostly unknown.In species-rich models, a specific type of multiple stable states has been found to appear when heterogeneous competition is at play [7,12,16,18].While simple models typically require specific non-linearities for multistability to appear, these so-called cliques can arise even under linear interactions in species-rich models.
In the present work, we have investigated cliques in the GLV model.Our results provide both an explanation for their occurrence and a prediction for their ecological properties.Cliques are subsets of the original ecosystem that fulfill a set of simple dynamical and statistical constraints.These predictable signatures provide potential benchmark tests to confront experimental evidence of multistability.Our work, therefore, brings forward an explanation and a testable framework of multistability in empirical species-rich systems.

Figure 1 .
Figure 1.Multistability in the GLV model.The GLV model provides a benchmark model to understand species-rich ecosystems.Its rich phase space can be visualized by measuring the fraction of initial conditions that end up in a stable attractor[12,20].In (A) we plot this fraction as a function of mean and standard deviation of interaction strengths A ij for systems of S = 50 species (see Methods).The GLV model has a rich phase space with four dynamical phases (divergence, coexistence, competitive exclusion and cliques+chaos).The circles mark two known intersections of the line separating the stable coexistence from the chaos and cliques regimes.In (B) we present an example of the complex dynamics of the cliques+chaos phase.We simulate the dynamics with µ = −0.5, σ = 0.5 and perturb population abundances every ∆t = 10 3 timesteps.Chaotic fluctuations, different stable community states and cycling behavior can be observed in a single system (see supplementary material section I.A for details on simulating panels (A) and (B)).

Figure 2 .
Figure 2. Interaction patterns in cliques.In (A) we show that mean interaction strength in cliques is less competitive than in the original system (µ ′ > µ) and heterogeneity is lower (σ ′ < σ).Each dot represents the difference µ − µ ′ , σ − σ ′ between an original system A with S = 50 and a stable clique A ′ .The equivalent comparison with mathematical tests can be found in SM II.A. Cliques appear as subsets that harbor less competitive and less heterogeneous interactions than their system of origin.In (B) we test the theoretical feasibility condition from[30] see (6)).We find that the interaction strengths in cliques (A ′ ij ) correlate with species abundances (x ′ * i ), creating the predicted diffuse pattern.Roughly put, abundant species in a clique (larger circles in the community sketch) need to compete weakly (thinner interaction lines).Instead, less abundant species (smaller circles) can compete strongly (thicker lines) without impairing coexistence[30].Each dot represents an A ′ ij value and its mathematical expectation, while the dark line indicates the 1 : 1 relation expected from the theory.For both panels we plot in red the linear regression of the data to help visualize its qualitative trend.

Figure 3 .
Figure 3. Species diversity in cliques is bounded by µ and σ.In (A) we plot the observed mean interaction strength of cliques (µ ′ ) in terms of their diversity (S ′ ).Light gray dots are the statistics of randomly sampled elements of A. Colored dots indicate the observed µ ′ values of cliques appearing in a system with µ = −0.5 (dashed gray lines) and σ = 0.5.The maximum and minimum competition strengths in cliques are predictable and depend only on the parameters of the original distribution µ, σ and the clique size S ′ .Equating the two expressions allows us to compute the maximum diversity (S ′ M ) that a clique could ideally attain given the statistics of the original pool.In (B) we show that, even if we increase system size, the maximum observed diversity of cliques (max(S)) remains close to the predicted S ′ M (µ = −0.5, σ = 0.5) ≈ 12.In (C) we show the maximum size of a clique for a system close to the coexistence-to-cliques transition (fixed µ = −0.5, decreasing σ = 0.5 √ 2/S + 10 −2).The maximum clique diversity scales as half of the total diversity, consistent with the diversity of the coexistence state before the critical transition[12].

Figure 4 .
Figure 4. Species migration decreases clique stability.We plot the fraction of stable simulations against system size with (blue dots) and without (orange dots) species migration.We compare these with two theoretical predictions: the presence of successful invaders on cliques generated from the dynamics of the GLV model (dark line) and on cliques filtered by the mathematical procedure (gray line).Larger system sizes destabilize cliques due to adding invaders (blue dots, dark line).Randomly selected stable subsets do not accurately mimic the composition and resistance to invasions of cliques (gray line).