article: The consensus that complexity begets stability in ecosystems was challenged in the seventies, a result recently extended to ecologically-inspired networks. The approaches assume the existence of a feasible equilibrium, i. e. with positive abundances. However, this key assumption has not been tested. We provide analytical results complemented by simulations which show that equilibrium feasibility vanishes in species rich systems. This result leaves us in the uncomfortable situation in which the existence of a feasible equilibrium assumed in local stability criteria is far from granted. We extend our analyses by changing interaction structure and intensity, and find that feasibility and stability is warranted irrespective of species richness with weak interactions. Interestingly, we find that the dynamical behaviour of ecologically inspired architectures is very different and richer than that of unstructured systems. Our results suggest that a general understanding of ecosystem dynamics requires focusing on the interplay between interaction strength and network architecture.
  A central question in ecology is to understand the factors and conditions that ensure ecological systems to persist, a requisite for the sustained provisioning of vital ecosystem services. This question of a “balance of nature” has a long history in science [1,2], and the consensus that “complexity begets stability” emerged among ecologists in the fifties. MacArthur [3] had a radical view on this question, arguing that stability will increase with the two fundamental ingredients of complexity, the number of species and of interactions. The argument for this claim was borrowed from Odum [4]: stability increases with the number of paths through which energy can flow up in a food web. Later, Elton [5] provided a suite of arguments for this positive relationship. The first one is the following: mathematical systems composed of one predator and one prey exhibit conspicuous fluctuations. Implicit in this argument is that more complex systems should be more stable, which remained untested at that time. In the seventies, Levins [6,7], Ashby and Gardner [8], and May [9] showed numerically that large random systems may be expected to be stable up to a certain connectance threshold, contradicting the earlier ideas that complex natural systems are more likely to be stable. In his impactful work, May [10,11] showed mathematically using random matrix theory [12,13] that large and random ecosystems are inherently unstable. His approach was based on a mathematical study of community matrices, which represented unstructured random networks of interacting species. He used a local-stability analysis assuming these systems were at equilibrium. He derived a simple and elegant criterion for system stability, which is a milestone in the stability-complexity debate [14]. May [11] concluded that there was no comfortable theorem assuring that increased complexity will lead to stable systems, and that the task was therefore to “elucidate the devious strategies which make for the stability in enduring natural systems”. Recently, the work of May was revisited by Allesina and collaborators [15,16]. They established stability criteria for systems where species interact specifically via either competition, mutualism, or predation. As such, their contribution is a refinement of May’s approach that considered mixtures of interaction types. The general conclusion is that, in species-rich communities, interactions should be moderate to ensure the stability of the equilibrium. They also performed simulations to study stability in randomly generated systems whose architecture mimics empirically observed food webs [17,18]. They arrived at the counter-intuitive conclusion that such ecologically-inspired structured graphs are less stable than unstructured ones. They also focused on interaction weights, and interestingly found that weak interactions should increase the stability for mutualistic and competitive webs, but decrease the stability of food webs (see also [19], where more precise statements are obtained for interactions weights of different intensity and symmetry). It was also found that more realistic structures seem to be detrimental for stability, and that the structure alone plays a minor role for stability compared to interaction weights [20,21]. All the above approaches are local-stability analyses, where systems are linearised at the equilibrium point and stability is evaluated only in the close vicinity of this point. In ecological systems, this point is meaningful only if the equilibrium-abundances are all strictly positive; in other words, if the equilibrium is feasible. Fig 1 illustrates the fact that a stable equilibrium of a simple two-level food web can be feasible when the absolute value of the interspecific competition coefficient is small, but infeasible when it becomes too large. Hence, a stable and feasible equilibrium can be transformed into an infeasible stable equilibrium by changing competition coefficients. In these previous works, equilibria were simply assumed feasible, without further analysis; basically, the true Jacobian, evaluated at the equilibrium, was replaced by a random matrix. Interestingly, soon after the work of May [11], Roberts [22] noted that May’s approach indeed remained silent concerning the feasibility of the equilibria. In a simulation study, he found that, in feasible systems, stability increased with the number of species. However, Roberts did not explore the very question of the relationship between system size and feasibility probability. This is a key issue since feasibility is a prerequisite to local-stability analysis. In the particular case of competitive systems, this question had to await the work of Logofet [23] who found that “equilibriumness” was vanishing with species richness. Note that already in 1970, Vandermeer [24] studied a question very related to feasibility, the expected number of species that can coexist in competitive communities. The same question has been treated by Rossberg [25], where the coefficient of variation of species abundance is studied in a mean-field approach. Since then, few contributions have explored the question of feasibility, and if so in particular systems and with regard to the characteristics of the interactions rather than to species number (e. g. , for asymmetric competition [26,27] or for mutualism [28]). Very recently the question of feasibility reappeared with the aim of providing general criteria leading to feasible equilibria [29], following an analogue approach as in ref. [28]. The idea of this important work is to compute the fraction of demographic parameters leading to feasibility under given conditions. Very interestingly, from the perspective of system complexity whose effect on stability has been thoroughly investigated in previous works [4,6, 8,9, 15], keeping this fraction well-defined (different from zero and infinity) when the number of species increases requires that the absolute mean interaction strength is small enough [29, Fig 2]. Here, we present an extensive study of the fundamental question of equilibrium feasibility, and expose its underlying mathematical principles in order to investigate how complexity influences feasibility probability. Two classes of networks are studied: 1) networks which do not possess any particular topological structure, for which we consider random, mutualistic, competitive, and predator-prey interactions, and 2) predator-prey networks designed to capture food web architecture, following the cascade [17], niche [18] and nested-hierarchy [30] models. We show that, in situations compatible with the above works [10,15,16], the existence of a feasible equilibrium is not guaranteed, since the probability to observe such a point decreases exponentially towards zero as the size of the network grows. Interestingly, we find that with weak interactions, non-trivial feasible equilibria are found only in ecologically-inspired systems, in a way that is ecologically sensible.
  We will first focus on webs whose topologies are based on Erdös-Renyi random graphs (unstructured webs), where edges occur with probability C, the so-called connectance of the graph. It is given by C = L/ (S ⋅ (S − 1) ) with L the total number of links in the network. Next, we will consider graphs whose topologies are drawn randomly according to the cascade [17], the niche [18] and the nested-hierarchy food web models [30] (structured webs), that capture more accurately trophic interactions of real food webs. Fig 2 provides illustrations of these various topologies. If species i and j interact; that is if the web contains an edge between these two species, the nature of the interaction is encoded in the sign of the coefficients aij: aij < 0 and aji < 0 for competition; aij > 0 and aji > 0 for mutualism; aij > 0 and aji < 0 for predation of i on j. The interactions aij form the S × S interaction matrix A. In order to exhibit the effects of network structure on the dynamics, the aijs are randomised with fixed mean and standard deviation. Their magnitude plays a key role for the local stability of networks [32,33]. We introduce a parameter 0 ≤ δ ≤ 1 ruling the average interaction strength. The aijs are divided by the linkage density (a measure of complexity given by CS [10]) raised to the power δ (see ref. [28,34]). This normalisation is mathematically sensible since Wigner’s theory [12], on which May’s and Allesina’s results are based, is built for the case δ = 1/2. Analytically, we will show that three regimes emerge when S becomes large: strong (0 ≤ δ < 1/2), moderate (δ = 1/2), and weak interactions (1/2 < δ ≤ 1), see Fig 3a. Intraspecific competition is included as customarily [10,11,15] with a common coefficient θ < 0, separated from the interaction matrix and unaffected by the normalisation constant. The first type of networks (unstructured) are based on Erdös-Renyi graphs. For May’s random model [10], an entry of the adjacency matrix α is set to one with probability C, zero otherwise. For mutualistic, competitive and predator-prey interactions, a pair (i, j), i < j, is linked such that αji = αij = 1 with probability C. The adjacency matrices for the second type of networks (structured networks modeling predation) are built according to the cascade [17], niche [18] and nested-hierarchy [30] models. Their adjacency matrix α indicates that i preys upon j when αij = 1 and their connectance is set to C/2. Their detailed construction is reported in S1 Text. See also Fig 2 for an illustration. In May’s random model [10], interactions strengths are sampled from i. i. d. random variables, leading to mixture of interactions types, including also commensalism and amensalism. Concerning models for mutualism (competition), aij and aji are independently randomly chosen from a positive (negative) random variable. For predation models, any pair (i, j) such that i < j is linked (i. e. αji = αij = 1 in the adjacency matrix) with probability C and the interactions strengths are also independently sampled, with the restriction that sign (aij) = −sign (aji). For the structured networks, the entries of α are multiplied with positive i. i. d. values and the entries of αt, the transpose of α, with negative i. i. d. values. The interaction matrix A is the addition of the two. Gaussian or folded-Gaussian random variables with standard-deviation σ are chosen in our simulations, but our analytical results are not limited to the particular choice of these distributions, see S1 Text. Our simulations are performed by randomly sampling A and then computing the equilibrium x* according to Eq (1). This Monte-Carlo procedure is used for the estimation of PS under the different scenarios.
  A system is feasible if all abundances at equilibrium are positive. If this equilibrium exists, the vector x* corresponds to a point for which the dynamics of the Lotka-Volterra system stops varying, and is given by x * = - (θ I + (C S) - δ A) - 1 r. (1) If the solution to this equation involves negative abundances, then the system is not feasible. Since the interactions aij are random variables, it is natural to study PS, the probability to observe a feasible equilibrium, which is given by P S = P (x i * > 0, for all species i). (2) This probability depends on θ, C, S, δ, A and r. For a given δ, C, and type of network, we choose the growth rates vector r = (ri) most favourable for feasibility, or, in case of trivially feasible equilibria, ri to be independent and bounded random variables. The most favourable growth rates correspond to the deterministic vector making the equilibrium to be, on average, as far as possible from the boundaries of the feasibility domain. This deterministic vector is referred to as the mean structural vector. By definition, the feasibility domain is the positive orthant of the phase space. We show in the S1 Text that this vector is related to a particular way of choosing intrinsic growth rates, which have been used in ref [20] to avoid negative abundances. It is now possible to compute PS with regard to the randomness of the aijs. Extending results on large systems of random equations [35], we develop analytical formulas for PS for our first class of unstructured models when S tends towards infinity, and use Monte Carlo simulations to estimate PS for the second class of structured models (see S1 Text). Before exposing our results, we briefly discuss the link between stability and feasibility. It is only when A is such that x* is feasible that the local-stability analysis of x* is sensible. This analysis involves the Jacobian matrix evaluated at x* (the so-called community matrix), J (x *) = diag (x *) (θ I + (C S) - δ A), (3) which of course depends on the random vector x* and on the random interactions matrix A. The equilibrium x* is linearly stable when the real parts of all eigenvalues of J (x*) are negative. Under May’s approach, and similarly in ref. [15,16,19], J (x*) is replaced by the random matrix J ˜ = θ I + (C S) - δ A. All the information on x* and on the inherent relations between x* and A (see Eq (1) ) are overlooked. The problematic point is that one obtains stability results by focusing on the eigenvalues of J ˜ even when a feasible equilibrium x* does not exist. For example, May [10] considers that the elements of A are such that all diagonal entries are set to 0, and that all off-diagonal random elements are independent and set to 0 with probability 1 − C. A non-zero entry is sampled from any centred distribution with standard deviation σ. Applying results from random matrix theory [13], the eigenvalues of J ˜ will then have negative real parts for large S when σ (C S) δ C S < | θ |, (4) which is May’s stability condition. We show in S1 Text, Proposition S. 4. 1, that May’s criterion still holds for the matrix J (x*) under the additional assumption that x* is feasible (see Fig 3b). Thus it appears critical to study the feasibility of such systems. May’s paradox that complexity decreases stability in mathematical models led to the exploration of the effects of topological structure on systems’ dynamics. Several works have shown that particular architectures can stabilise ecological networks [37–40]. Nevertheless, the feasibility of different network architectures for different interaction regimes has not been explored. Consider for example Cohen’s cascade model of trophic networks [17], where species are hierarchically ordered so that they feed only on lower indexed species, generating sub-triangular matrices with random structure. For the three regimes of δ, simulations yield similar behaviour of PS as for unstructured models (Table 1). We now test two other food web models that more accurately capture the structure of trophic interactions in real systems, the niche [18] and the nested-hierarchy models [30]. The former generates purely interval food webs (predators consume all species in a niche interval) while the latter is based on evolutionary processes and relaxes this constraint. The results are qualitatively similar to the ones of unstructured random models for all types of interactions. However, two interesting features emerge. First, the mean structural vector now contains negative growth rates, which is the essence of predator-prey systems where predators will die in the absence of prey (this feature also appears in the cascade model, but not in the unstructured food web model, see Table 2 in S1 Text). Second, with weak interactions, each x i * no longer converges to a deterministic limiting vector, but instead to a non-trivial random variable having a distribution of finite support and positive standard-deviation. Thus, contrary to the unstructured and cascade models, a particular realisation of the interaction matrix A will influence the equilibrium values. Indeed, these highly structured networks and their related mean structural growth rates induce correlations among the abundances at equilibrium, which prevent a convergence to a deterministic value. We find that the size of the support depends on the hierarchical position of the species (Fig 6). We performed additional simulations based on empirical food webs, and find similar results (Figs G and H in S1 Text). Since the mean structural vector is constructed for the purpose of maximizing the probability of feasibility, we have also considered random Gaussian growth rates of positive mean to complete our analysis. For the niche and nested-hierarchy models, the standard deviations σ i * of the abundances at equilibrium are ordered according to the trophic position of the species within the food web (from basal autotrophs, to intermediate species that are predator and prey, and to top predators). We find σ basal * < σ intermediate * < σ top *, a result similar to that obtained with the structural vector, see Figs 3e and 7. This effect, also visible in empirical food webs (Figs G and H in S1 Text), disappears in two particular cases: 1) for zero mean Gaussian random growth rates; 2) when the growth rate standard-deviation becomes large. In both cases, the equilibria then become independent of species trophic position (Figs G and I in S1 Text).
  Our results can be summarised as follows: for the cases compatible with May’s [10] and Allesina’s [15,16] framework where δ = 1/2 (moderate interactions), there exist almost surely no feasible equilibria in species-rich systems for any of the models considered. For weak interactions in randomly structured systems, feasibility is granted; also, the stability criteria are trivially asymptotically satisfied. Therefore, in such situations, the criteria are basically void of information. However, including realistic structure in trophic systems, we find that the results are no longer trivial. As a general message, future work tackling dynamical stability must be preceded by a feasibility analysis. In our Lotka-Volterra framework, this involves knowledge of the structure and intensity of interactions, and of the growth rates. Feasibility can then be easily evaluated. Our results indicate that sensible information on the dynamics of natural systems will be obtained by concentrating on the interplay between ecological structure, strengths of interaction weights, and how interaction strengths are distributed in their architecture. There are still many other aspects that should be taken into account to reach general conclusions about feasibility and stability that better fit natural systems. They are very different from unstructured systems, and many constraints affect their architecture and dynamics. A first question simply concerns system size in terms of species richness S. Most described systems are middle sized, but large webs do exist [41,42], so that mathematical study focusing on large S limits can be justified. An intriguing aspect here is how smaller and well-definable webs are dynamically embedded in larger species-rich systems. For example, it has been found that top predators have a stabilising effect on food webs by coupling fast and slow energy channels in natural aquatic systems [43]. This question could be framed in a general theory for the “inverse pyramid of habitat”, which explores the consequences of the ubiquitous (but often overlooked) observation that species located at higher trophic levels tend to have larger home ranges [44]. In this respect, the spatial structure of the system is relevant and has non-trivial consequences for its dynamics [45,46], a research framework that is currently expanding (e. g. , [47–49]). Relevant to system size lies also the question of how network nodes are defined: taxonomic resolution is often heterogeneous in observed food webs, with basal species being more often pooled (e. g. , phytoplankton being considered as a single species). Secondly, species richness S and connectance C are usually treated as independent parameters, while it has been described that C decreases with S in different systems [50], thus restricting the parameter space to be explored. Finally, even if topology per se has been suggested to play a minor role in system stability [20], we showed that network architecture does play a non-trivial role on equilibrium abundances. Layer architectures as those produced by the niche and nested-hierarchy models yield ecologically sensible results; exploring other structural models [21,51–53] and especially natural architectures is necessary. Apart from these structural considerations, ecologists are confronted by a more elusive issue, the estimation of demographic parameters. To make sense, the investigation of feasibility and stability should be based on meaningful parameter values. The growth rate vector r that appears in Lotka-Volterra dynamics is mostly unknown or very difficult to obtain experimentally. One can either keep this vector as an intrinsic free parameter or assume that such growth rates change with interaction weights, as in [20] or with our mean structural vector. The estimation of interaction parameters is also far from trivial [54–57]. Experimental studies (see, e. g. , [20,58,59]) choose interaction weights indirectly, according to particular methods like predator-prey mass ratio models, biomass flux, or the Ecopath method (see, e. g. , [60–62]). Additionally, one unresolved question is how intraspecific competition scales with average interaction strength, which is perhaps the most important ingredient for stability (e. g. , [10,15,58,63]), and also plays a role in feasibility (see Proposition S. 3. 2 in S1 Text). Also, if the magnitude of intraspecific competition can be assumed to be independent of S, this is likely not the case with interspecific interactions (captured by our parameter δ), which has been shown to play a key role in mutualistic networks [28]. This question is related to the study of weak interactions, which have been shown to promote stability [19,32,58,59] and here feasibility. Finally, one must consider that these parameters do scale with body size [64], which defines a constraint between parameter values and network structure, as large-bodied species typically populate higher trophic levels. Interestingly, this allometric relationship may underlie the result that stability depends on the link between trophic position and interaction strength [20], and our finding that abundance variability scales with trophic position. For mathematical tractability, many studies on the dynamics of food webs, including ours, rely on Lotka-Volterra models with mass-action type interactions (the so-called Holling type I functional response). However, investigations on the number of prey eaten per predator and unit time in simple systems do not support this modelling assumption. There exists a vast literature on the subject and the choice of a sensible functional response is still debated [65,66]. In Section S. 5 of the S1 Text, we studied a predation unstructured model with a Holling type II functional response. We found, independently of the choice of the regime δ, that the system behaves similarly to the Lotka-Volterra model with weak interactions. Here, the saturation inherent to Holling type II functional response makes interactions weak, irrespectively of the parameter δ. It shows that interaction dynamics can strongly affect system feasibility in predator models. Other aspects related to interspecific interactions affect community dynamics, e. g. , evolving trait-mediated direct and indirect interaction weights [67–69], ontogenetic shifts in interactions [70], or high-order interactions whereby the interaction between two species is modulated by one or more other species [71]. Also, it has been shown that “foraging adaptation” by predators affects the dynamics of food webs [72]. Such non-linearities in species interactions have been explored from the stability point of view of systems adjusted to be feasible [73–75]. However, since the precise way interactions are modelled has non-trivial consequences on system dynamics, the relationship between feasibility and complexity should also be considered in such situations. May’s approach considered systems with any kind of interaction and not only trophic ones. An articulate answer on the relationship between stability and complexity should obviously incorporate all types of interactions relevant for system dynamics (including the effects of ecosystem engineers [76,77]). This question has been the focus of recent theoretical developments [78,79], and empirical studies specifically addressing this question start to emerge [80]. Science faces increasingly complex situations where high-dimensional parameters occur. A good example is statistical mechanics, for which relevant results have been obtained for highly complex systems without having a precise knowledge of the microscopic details of a model. Our results like many cited here follow the same line. However, due to the complexity of natural systems, the gap between theoretical and empirical investigations is likely to remain open. To attain a consensus, the feasibility of the model systems should not be forgotten [23,24,26–29], and more effort must be devoted to obtaining empirical and experimental time series as necessary benchmarks.
summary: Finding conditions leading to stable equilibria, i. e. equilibria robust to small perturbations, is a key feature in understanding how persistent are ecological systems. The number of species involved in the ecosystem as well as the number of relations between them has historically been pointed out to be the prominent ingredients determining stability. The latter is however always computed for a particular equilibrium, which needs and is usually assumed to contain only positive abundances. This assumption has until now only barely been tested in relation to system size. We provide here analytical results complemented by numerical simulations which show that it is almost impossible to find equilibria containing only positive abundances in species rich systems or, if so, the parameters are constrained to lead to ecologically-nonsensical abundances in unstructured systems. Interestingly, only ecologically-inspired architectures permit one to obtain feasible equilibria with a rich and sensible dynamical behavior.
keywords: ecology and environmental sciences
  ecological niches
  particle physics
  systems science
  mathematics
  theoretical ecology
  ecological metrics
  computer and information sciences
  system stability
  food web structure
  fundamental interactions
  species diversity
  elementary particle interactions
  community ecology
  physics
  ecology
  predation
  trophic interactions
  weak interaction
  biology and life sciences
  physical sciences
  species interactions
title: The feasibility of equilibria in large ecosystems: A primary but neglected concept in the complexity-stability debate