Search tips
Search criteria 


Logo of ncommsLink to Publisher's site
Nat Commun. 2017; 8: 0.
Published online 2017 February 24. doi:  10.1038/ncomms14389
PMCID: PMC5333123

Feasibility and coexistence of large ecological communities


The role of species interactions in controlling the interplay between the stability of ecosystems and their biodiversity is still not well understood. The ability of ecological communities to recover after small perturbations of the species abundances (local asymptotic stability) has been well studied, whereas the likelihood of a community to persist when the conditions change (structural stability) has received much less attention. Our goal is to understand the effects of diversity, interaction strengths and ecological network structure on the volume of parameter space leading to feasible equilibria. We develop a geometrical framework to study the range of conditions necessary for feasible coexistence. We show that feasibility is determined by few quantities describing the interactions, yielding a nontrivial complexity–feasibility relationship. Analysing more than 100 empirical networks, we show that the range of coexistence conditions in mutualistic systems can be analytically predicted. Finally, we characterize the geometric shape of the feasibility domain, thereby identifying the direction of perturbations that are more likely to cause extinctions.

Natural populations are faced with constantly varying environmental conditions. Environmental conditions affect physiological parameters (for example, metabolic rates1) as well as ecological ones (for example, the presence and strength of interactions between populations2,3,4,5). Therefore, in order to persist, ecological communities necessarily need, at the very least, to be able to cope with small environmental changes. Mathematically, this translates into an argument on the robustness of the qualitative behaviour of an ecological dynamical system: to guarantee robust coexistence, a model describing an ecological community needs at least to be (qualitatively) insensitive to small perturbations of the parameters6,7. This notion has been formalized in the measure of robustness8 or structural stability9, expressed as the volume of the parameter space resulting in the coexistence of all populations in a community.

While the local asymptotic stability (the ability to recover after a small change in the population abundances) of ecological communities has been studied in small10 and large11,12,13,14 systems, the study of structural stability (that is, the ability of a community to retain the same dynamical behaviour if conditions are slightly altered)—despite being proposed early on as a key feature in the context of the diversity–stability debate15,16,17,18—has historically been restricted to the case of small communities, with the first studies of larger communities appearing only recently9,19, and—because of mathematical limitations—dealing exclusively with the case of large mutualistic communities. Studies of structural stability have so far focused on the effect of ecological network structure (who interacts with whom) on the volume of parameter space leading to feasible equilibria, in which all populations have positive abundances.

Here we develop a geometrical framework for studying the feasibility of large ecological communities. We overcome the limitations that have hitherto prevented the study of consumer–resource networks, thereby providing a unified view of feasibility in ecological systems. Using a random matrix approach (which helped identify main drivers of local asymptotic stability), we pinpoint the key quantities controlling the volume of parameter space leading to feasible communities, as well as its sensitivity to changes in these parameters. We then contrast these expectations for randomly connected systems with simulations on structured empirical networks, quantifying the effects of network structure on feasibility.


Theoretical framework

For simplicity, we consider a community composed of S species whose dynamics is determined by a system of autonomous ordinary differential equations:

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m1.jpg

where ni is the density of population i, ri is its intrinsic growth rate and Aij measures the interaction strength between population i and j. In this paper we consider only the linear functional response (that is, A does not depend on n). In Supplementary Note 13 we discuss how and under which condition one could generalize our results to nonlinear functional responses. A fixed point n* (that is, a vector of densities making the right side of each equation zero) is feasible if An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m2.jpg*>0 for every population. A fixed point is locally asymptotically stable if, following any sufficiently small perturbation of the densities, the system returns to a small vicinity of the fixed point. The fixed point is globally asymptotically stable if the system eventually return to it, starting from any positive initial condition within a finite domain. A system with a fixed point is structurally stable if, following a sufficiently small change in the growth rates ri, the new fixed point is still feasible and stable.

To study the range of conditions leading to stable coexistence, we need to disentangle feasibility and local stability. This problem is well discussed in ref. 9, where it was solved for the case of one possible parameterization of mutualistic interactions. If A is diagonally stable or Volterra-dissipative (that is, there exists a positive diagonal matrix D such that DA+ATD is stable), then any feasible fixed point is globally stable20,21. Unfortunately, a general characterization of this class of matrices is unknown22. We proceeded therefore by considering only the matrices such that all the eigenvalues of A+AT are negative (that is, the matrix A is negative definite in a generalized sense23, corresponding to D being equal to the identity matrix; see Methods and Supplementary Notes 1 and 2). This choice reduces the number of parameterizations one can analyse, as not all the diagonally stable matrices are negative definite. However, as shown in Supplementary Fig. 1, only very few parameter combinations are excluded from this set. Moreover, the effects of negative definitness are well studied for random matrices24, and by using it we can extend the study of feasibility to any ecological network, including food webs.

Our goal is to measure the fraction of growth rate combinations, out of all possible ones, that lead to the coexistence of all S populations. Since we can separate stability and feasibility, we only need to find those ri leading to feasible fixed points, and the condition above ensures that these will be globally stable. As pointed out before9, the problem is not to find a particular set of ri leading to coexistence, but rather to measure how flexibly one may choose these rates. As shown in Fig. 1, this quantity—indicated by Ξ henceforth—can be thought of as a volume, or more precisely a solid angle, in the space of growth rates25 (see and Supplementary Note 3).

Figure 1
Geometrical properties of feasibility.

To calculate Ξ, one might naively wish to perform direct numerical computation of the fraction of growth rates, leading to a feasible equilibrium. While a direct calculation is viable when S is sufficiently small, this procedure becomes extremely inefficient for large S (ref. 9). We introduce a method that can be used to efficiently calculate Ξ with arbitrary precision, even for large S (see Supplementary Note 4). Using this method, we can accurately measure the size of the feasibility domain, with larger values of Ξ corresponding to larger proportions of conditions (intrinsic growth rates) compatible with stable coexistence. For reference, we normalize Ξ so that Ξ=1 when populations are self-regulated and not interacting (Methods), that is, when the interaction matrix A is a negative diagonal matrix, and thus equation (1) simplifies to S independent logistic equations.

Feasibility is universal for large random matrices

May's seminal work11 pioneered the use of random matrices as a reference, or null model, of ecological interactions. A particularly interesting feature of random matrices is that the distribution of their eigenvalues (determining local stability) is universal26. This means that local stability depends on just a few, coarse-grained properties of the matrix (that is, the number of species and the first few moments of the distribution of interaction strengths) and not on the finer details (for example, the particular distribution of interaction strengths; see Supplementary Note 5). In fact, these moments can be combined into just three parameters: E1, E2 and Ec (Methods). Together with S, they completely determine local asymptotic stability.

We tested whether universality also applies to feasibility. We considered different random matrix ensembles obtained for different connectance values and distributions from which the matrix entries were drawn, but with constant values of S and of E1, E2 and Ec. We then checked whether the size Ξ of the feasibility domain depended only on these four quantities or also on finer details. Surprisingly, we found that the feasibility of random matrices is also universal (Methods Supplementary Note 5 and Supplementary Fig. 2). Two very different (random) ecosystems, with completely different interaction types and distributions of interaction strengths, but having the same number of species S and the same E1, E2 and Ec, have the same Ξ in the large S limit. This result has important theoretical implications, as it indicates those moments as the drivers of feasibility, but also very practical consequences, namely that the parameter space one needs to explore is dramatically reduced.

An analytical complexity–feasibility relationship

The universality of Ξ suggests that it is amenable to analytical treatment. As explained in Supplementary Note 6 and shown in Fig. 2, when the mean and variance of interaction strengths are not too large and in the limit of large number of species, we are able to derive the following approximation for Ξ for large random interaction matrices A:

Figure 2
Feasibility domain in random and empirical webs.
An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m3.jpg

where S is is the number of species, d is the mean of A's diagonal entries and E1=, the product of the connectance C and the average interaction strength μ (see Methods). A more accurate formula is presented in Supplementary Note 6.

In analogy with the celebrated result of May11 connecting stability and complexity, equation (2) can be considered as a complexity–feasibility relationship. While in May's scenario and in its generalizations12 the effect of complexity and diversity on stability is always detrimental, it does depend on the interaction type in the case of feasibility. Given that d is negative by construction, having more species or connections can either increase (E1>0) or shrink (E1<0) the size of the feasibility domain, as a function of the sign of interaction strenghts (see Fig. 2). It is important to stress that we computed Ξ under the assumption of A being negative definite. When we consider how Ξ depends on S and other parameters, we need to take into account the conditions making the matrix negative definite (see Methods and Supplementary Notes 2 and 5). In the case of positive interaction strengths, this condition is d+SCμ<0, implying an upper bound for μ that depends on S.

Analytical prediction of the feasibility of empirical networks

Having explored the feasibility of random networks, we proceed to investigate the effects of incorporating empirical network structure. Ecological networks are in fact non-random27,28,29, and many studies have hypothesized that the structure of interactions could increase the likelihood of coexistence30,31,32. Having an analytical prediction for random matrices, we can study whether it predicts the size of the feasibility domain for empirical networks as well. Figure 2 shows the simulated values of Ξ for 89 mutualistic networks and 15 food webs (Supplementary Note 8 and Supplementary Table 1), parameterized multiple times and compared with our analytical approximation (see Methods). We find that Ξ of empirical mutualistic networks is well predicted by our formula, while it overestimates the feasibility domain of food webs, indicating that their non-random structure has a strong negative effect on feasibility.

We compared the effect of the empirical structure of mutualistic networks with randomizations, by controlling for the interaction strengths (see Supplementary Note 9 and Supplementary Figs 5–9). We show that, in the absence of variability in interaction strengths, the structure of empirical mutualistic networks has a positive effect on feasibility, which is strongly reduced when interaction strengths are allowed to vary. While this effect of empirical mutualistic networks is statistically significant, its effect on Ξ is negligible compared with the effect of the mean interaction strengths, and can only be detected by controlling very precisely for interaction strengths (see Supplementary Note 9). On a broader scale, as shown in Fig. 2, the size of the feasibility domain of empirical networks is well predicted by our analytical formula.

On the other hand, the negative effect of food web structure on Ξ is substantial. We compare each network with randomizations and also with predictions of the cascade model27, which has recently been shown to predict well the stability of empirical food webs14 (see Supplementary Note 9 and Supplementary Figs 10–12). By analysing different parameterizations we found that the feasibility domain of empirical structures is consistently and significantly smaller than that of both the randomizations and the cascade model. For most of the webs, the prediction obtained from the cascade model is better than that of randomizations, suggesting that the directionality of empirical webs plays a role in reducing feasibility, with other properties of the structure of empirical networks also contributing significantly to feasibility.

Shape of the feasibility domain

So far, we have focused on the volume of the parameter space resulting in feasiblity. However, two systems having the same Ξ can still have very different responses to parameter perturbations, just as two triangles having the same area need not to have sides of the same length (Fig. 1). The two extreme cases correspond to (a) an isotropic system in which if we start at the barycentre of the feasibility domain, moving in any direction yields roughly the same effect (equivalent to an equilateral triangle); (b) anisotropic systems, in which the feasibility domain is much narrower in certain directions than in others (as in a scalene triangle). For our problem, the domain of growth rates leading to coexistence is—once the growth rates are normalized—the (S−1)-dimensional generalization of a triangle on a hypersphere. For S=3, this domain is indeed a triangle lying on a sphere as shown in Fig. 1. If all the S(S−1)/2 sides of this (hyper-)triangle are about the same length, then different perturbations will have similar effects on the system. On the other hand, if some sides are much shorter than others, then there will be changes of conditions which will more likely have an impact on coexistence than others. We therefore consider a measure of the heterogeneity in the distribution of the side lengths (Fig. 1 and Supplementary Note 10). The larger the variance of this distribution, the more likely it is that certain perturbations can destroy coexistence, even when Ξ is large and the perturbation small. This way of measuring heterogeneity is particularly convenient because it is independent of the initial conditions. Moreover, the length of each side can be directly related to the similarity between the corresponding pair of species (Supplementary Note 10), drawing a strong connection between the parameter space allowing for coexistence and the phenotypic space. As in the case of Ξ, this measure is a function of the interaction matrix and corresponds to a geometrical property of the coexistence domain.

While Ξ is a universal quantity for random networks, the distribution of side lengths is not: it depends on the full distribution of interaction strengths (Supplementary Note 10). On the other hand, it is possible to compute it analytically in full generality, that is, for any distribution of interaction strengths and any interaction types. In particular, we are able to obtain an expression for its mean and variance, which depend only on S, E1, E2 and Ec (Supplementary Note 10). Figure 3 shows that the analytical formula, in the case of random A, matches the observed mean and variance of side lengths of random networks perfectly.

Figure 3
Distribution of side lengths in random and empirical networks.

Empirical feasibility domains have more heterogeneous shapes

As we have done for Ξ, we can now test how non-random empirical network topologies influence the distribution of side lengths. Figure 3 shows that empirical food webs and, in particular, empirical mutualistic networks display a much larger variation in side lengths than expected by chance. This result is particularly relevant, indicating that even if the feasibility domains of empirical mutualistic networks are equal or larger than those of random networks, their shapes are less regular than expected by chance, and thus we expect perturbations in certain directions to quickly lead out of the feasible domain of growth rates.


A classic problem in mathematical ecology is determining the response of systems to perturbations of model parameters. In the community context, one important application is getting at the range of parameters allowing for species coexistence33,34,35. Several methods exist for evaluating this range7,36,37,38, but they either rely on raw numerical techniques or else can only evaluate system response to small parameter perturbations. Here in the context of the general Lotka–Volterra model, we have given a method for the global assessment of all combinations of species' intrinsic growth rates compatible with coexistence—what we have called the domain of feasibility. Our geometrical approach can determine not only the total size of the feasibility domain, but also its shape: it is always a simply connected domain forming a convex polyhedral cone whose side lengths can be evaluated from the interaction matrix. Applying our method to empirical interaction networks, we were able to characterize the region of parameter space compatible with coexistence; the importance of this kind of information is underlined by a rapidly changing environment that is expected to cause substantial shifts in the parameters influencing these systems.

The geometrical framework we employed, pioneered by Svirezhev and Logofet25, allows for the formulation of a complexity–feasibility relationship. In analogy with the celebrated complexity–stability relationship, it relates the size of the feasibility domain with diversity, connectance and interaction strengths of a random interacting community. While communities are not random, this relationship sets a null expectation for the scaling of the proportion of feasible conditions. We obtain that the mean of interaction strengths sets the behaviour of feasibility with the number of species. If the mean is negative (for example, in case of competition or predation with limited efficiency), the larger the system is, the smaller is the set of conditions leading to coexistence, while for positive mean (for example, in the case of mutualism) the converse is true.

Here we have shown that the fraction of conditions compatible with coexistence is mainly determined by the number and the mean strength of interactions. In terms of network properties, the relevant quantity is the connectance, with other properties (for example, nestedness or degree distribution) having minimal effects. In particular, once the connectance and mean interaction strength are fixed, the matrices built using empirical mutualistic networks have feasibility domains very similar to that expected for the random case, as was also observed previously in a similar context39.

The empirical network structure of mutualistic networks has a statistically significant effect on the size of the feasibility domain. Whether this effect is ecologically relevant depends on the specific application at hand. For instance, the effect of structure could be neglected to quantify how the feasibility domain would change if a fraction of pollinators went extinct, and it could be evaluated using our analytical result. In contexts where the interaction strengths are strongly constrained, structure would play an important role. Our method provides, in this respect, a direct way of quantifying the importance of different factors, disentangling the way different interaction properties affect feasibility.

For mutualistic interaction networks, our results clearly show which properties determine the global health of the community, and therefore indicate which properties should be measured in the field. While not observing a link or measuring a wrong interaction coefficient could have strong effects on ecosystem dynamics, they have very little effect on the size of the feasibility domain and how the community copes with environmental perturbations and how likely extinctions are40. The major role is played by corse-grained statistical properties of the interactions, such as connectance or the mean and variance of the interaction strengths.

For food webs, on the other hand, empirical systems tend to have feasibility domains smaller than either their random counterparts or models conserving the directionality of interactions (cascade model). It is an open question which properties of real food webs are responsible for restricting the feasibility domain in this way. A possible candidate is the group structure observed in food webs41, corresponding to larger similarity of how certain species interact with the rest of the system than expected by chance, which in turn reduces the size of the feasibility domain.

These results parallel those for the distribution of the side lengths of the convex polyhedral cone delimiting the feasibility domain. The variance of side lengths for empirical structures is much higher than that in random networks. This implies that even if the total size of the feasibility domain is large, it will have a distorted shape that is very stretched along some directions and shortened along others (Fig. 1). Consequently, it will be possible to find growth rate perturbations of small magnitude that will drive the system outside its feasibility domain42.

We have shown that each side of the feasibility domain corresponds to a pair of species, with the length determined by how similarly the two species interact with the rest of the system. As two species interact more and more similarly (that is, have a larger niche overlap), the corresponding side becomes shorter and shorter, which in turns means greater sensitivity to parameter perturbations. Consistently with earlier results7,8, this fact establishes a relationship between niche overlap and the range of conditions that lead to coexistence: greater niche overlap means a more restricted parameter range allowing for coexistence, irrespective of the details of the interactions.

Several recent lines of work have studied the effect of network structure on coexistence in species-rich communities, with contrasting results9,30,31,32,43,44. For instance, on one hand nestedness was shown to increase biodiversity31, while, on the other hand, it is known to be associated with lower stability43,45. The differences between the size and shape of the feasibility domain shed light on these contrasting results. Most of these studies rely on numerical integration, and therefore strongly depend on initial conditions. Given the difference in the shape of the feasibility domains of random and empirical networks, different initial conditions and their perturbations could result in markedly different outcomes: the feasibility domain could appear to be large or small depending on the direction in which perturbations are made.

Our characterization of the geometrical properties of the feasibility domain contributes to the complete picture of the relation between feasibility and stability. It has been recently proposed that nestedness promotes larger feasibility domain sizes over stability19, suggesting the existence of a trade-off between feasibility and stability. As we showed, the (mild) increase of the feasibility domain size parallels with the increase of the variability of side lengths. The latter property is crucial to quantify the robustness to perturbations, and it might be interesting to explore more carefully the relation between stability and the shape of the feasibility domain.

Having established the general geometrical properties of the feasibility domain, we are in a much better position to critically evaluate the feasibility domains of real ecological communities. We consider this as a first step along the way of describing feasibility in more complex models and ecological scenarios.


Disentangling stability and feasibility

From equation (1), a feasible fixed point, if it exists, is given by the solution of

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m4.jpg

where the asterisk denotes equilibrium values. A fixed point is locally asymptotically stable if all eigenvalues of the community matrix

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m5.jpg

have negative real parts. As discussed in Supplementary Note 2, if A is diagonally stable or Volterra-dissipative (that is, there exists a positive diagonal matrix D such that DA+ATD is stable), then a feasible fixed point is globally stable in An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m6.jpg.

A general characterization of diagonally stable matrices is unknown for more than three species22. There exist algorithms46 that reduce the problem of determining whether a S × S matrix is diagonally stable into two simultaneous problems of (S−1) × (S−1) matrices. While this method can be efficiently used to determine the diagonal stability of 4 × 4 matrices, it becomes computationally intractable for large S.

A matrix A is negative definite if

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m7.jpg

for any non-zero vector x. A necessary and sufficient condition for a real matrix A to be negative definite is that all the eigenvalues of A+AT are negative23. A negative definite matrix is also diagonally stable, as the condition for diagonal stability holds with D being the identity matrix. Since it is extremely simple to verify this condition and it has been characterized for random matrices, we will study feasibility of negative definite matrices. In Supplementary Note 5 and Supplementary Fig. 2 we show that with this choice we are excluding only a small region of the parameter space.

Size of the feasibility domain

The quantity Ξ is the proportion of intrinsic growth rates leading to feasible equilibria. While a more rigourous definition is presented in Supplementary Note 4, with a slight abuse of notation, Ξ can be thought of as

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m8.jpg

The factor 2S is an arbitrary choice that does not affect the results. It has been introduced to have Ξ=1 in absence of interspecific interactions (Aij=0 if ij in equation (1)) and when all the species are self-regulated (Aii<0 if ij in equation (1)). Given the geometrical properties of the feasibility domain, the proportion of feasible growth rates can be calculated considering only growth rate vectors of length one (Fig. 1 and Supplementary Note 3), as this choice does not affect the value given by equation (6). In Supplementary Note 4 we provide an integral formula for Ξ (refs 47, 48), which makes both numerical and analytical calculations possible.

Our method is still valid if some of the species are not self-regulated (that is, Aii=0 for some i). In Supplementary Note 7 we explicitly discuss the properties of the feasibility domain of a community with consumer–resource interactions. In that case, Ξ=0 either when the diversity of consumers exceeds the diversity of resources or in the absence of interspecific interactions. Since consumers are regulated by their resources, they cannot survive in their absence and should therefore be characterized by negative intrinsic growth rates. We observe indeed that a necessary condition for an intrinsic growth rate vector to be contained in the feasibility domain is to have negative values for the components corresponding to consumers.

Random matrices and moments

E1, E2 and Ec are moments of the random distribution for the off-diagonal elements of the interaction matrix, and are simply and directly related to the interaction strengths. They can be calculated as

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m9.jpg

For random networks with connectance C, these expressions reduce to (ref. 26)

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m10.jpg

where μ is the mean of the interaction strengths, σ is their variance and ρ is the average pairwise correlation between the interaction coefficients of species pairs26.

Universality of the size of the feasibility domain

The size of the feasibility domain should, at least in principle, depend on all the entries of the interaction matrix. When these elements are drawn from a distribution, the size Ξ of the feasibility domain is then expected to depend on all the moments of that distribution. As S increases, the dependence of Ξ on some of those moments and parameters might become less and less important. Ξ is universal if, in the limit of large S, it depends only on a few properties of the interaction matrix (that is, on just the first few moments of the distribution).

Specifically, for each unique pair of species (i, j), we set Aij=0 with probability 1−C and assign a random pair of interaction strengths (Mij, Mji)=(x, y) with probability C. The pair (x, y) is drawn from a bivariate distribution with given mean μ, variance σ and correlation ρ between x and y (ref. 26). By considering different bivariate distributions, we can analyse the effect of different sign patterns (for example, only (+, −) or (+, +) interactions) and different marginal distributions (for example, drawing elements from a uniform or a lognormal distribution).

Non-universality of Ξ would mean that it depends on all the fine details of the parameterization:

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m11.jpg

where f(·) is an arbitrary function. The dependence on μ, σ and ρ can, without loss of generality, be expressed in terms of E1, E2 and Ec:

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m12.jpg

However, if Ξ is universal, then for large S, it is possible to express it as a function of E1, E2 and Ec only:

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m13.jpg

To verify this conjecture, we calculated Ξ for matrices with the same values of E1, E2 and Ec that differed for the values of the other parameters. As extensively shown in Supplementary Note 5, Ξ is uniquely determined by S, E1, E2 and Ec (equation (2)).

Parameterization of mutualistic networks

The 89 mutualistic networks (59 pollination networks and 30 seed-dispersal networks) were obtained from the Web of Life data set (, where references to the original works can be found. Empirical networks are encoded in terms of adjacency matrices L: Lij=1 if species j interact with species i and 0 otherwise. When the original network was not fully connected, we considered the largest connected component.

In the case of mutualistic networks, the adjacency matrix L is bipartite, that is, it has the structure

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m14.jpg

where Lb is a SA × SP matrix (SA and SP being the number of animals and plants, respectively). The adjacency matrix contains information only about the interactions between animals and plants, but not about competition within plants or animals.

We parameterized the interaction matrix in the following way:

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m15.jpg

where the symbol o indicates the Hadamard or entrywise product (that is, (A o B)ij=AijBij), while WA, WAP, WPA and WP are all random matrices. WA and WP are square matrices of dimension SA × SA and SP × SP, while WAP and WPA are rectangular matrices of size SA × SP and SP × SA. The diagonal elements An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m16.jpg and An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m17.jpg are set to −1, while the pairs (An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m18.jpg, An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m19.jpg) and (An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m20.jpg, An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m21.jpg) are drawn from a bivariate normal distribution with mean μ, variance An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m22.jpg=cAn external file that holds a picture, illustration, etc.
Object name is ncomms14389-m23.jpg and correlation ρAn external file that holds a picture, illustration, etc.
Object name is ncomms14389-m24.jpg. Since these two matrices represent competitive interactions, μ <0. The pairs (An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m25.jpg, An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m26.jpg) were extracted from a bivariate normal distribution with mean μ+, variance An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m27.jpg=cAn external file that holds a picture, illustration, etc.
Object name is ncomms14389-m28.jpg, and correlation ρAn external file that holds a picture, illustration, etc.
Object name is ncomms14389-m29.jpg, where μ+>0. For each network and parametrization we computed the size of the feasibility domain Ξ.

We considered different values of μ, μ+, c, and ρ. Their values cannot be chosen arbitrarily, since A must be negative definite. For a choice of c, ρ, and a ratio μ/μ+, the largest eigenvalue of (A+AT)/2 is linear in μ+(as an arbitrary μ+ can be obtained by multiplying A by μ+ and then shifting the diagonal). Given the values of μ/μ+, c and ρ, one can therefore determine μmax, the maximum value of μ+ still leading to a negative definite A (that is, the value of μ+ such that the largest eigenvalue of (A+AT)/2 is equal to 0). Figure 2 was obtained by considering more than 1,000 parameterizations. Both the ratio μ/μ+ and the coefficient of variation c could assume the values 0.5 or 2, while the correlation ρ assumed values from the set {−0.9, 0.5, 0, 0.5, 0.9}. The value of μ+ was set equal to 0.25μmax and 0.75μmax. In addition, we considered the case μ =0.

Parameterization of food webs

In the case of food webs the adjacency matrix L is not symmetric, Lij=1 indicating that species j consumes species i. We removed all cannibalsistic loops. Since Lij and Lji are never simultaneously equal to one (there are no loops of length two), we parameterized the off-diagonal entries of A as

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m30.jpg

while the diagonal was fixed at −1. Both W+ and W are random matrices, where the pairs (An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m31.jpg, An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m32.jpg) are drawn from a bivariate normal distribution with marginal means (μ+, μ) and correlation matrix

An external file that holds a picture, illustration, etc.
Object name is ncomms14389-m33.jpg

We considered considering different values of μ, μ+, c and ρ. As explained above, given the values of μ/μ+, c and ρ, one can determine μmax, the maximum value of μ+ still corresponding to a negative definite A. Figure 2 was obtained by considering more that 350 parameterizations. Both the ratio μ/μ+ and the coefficient of variation c could assume the values 0.5 or 2, while the correlation ρ assumed either the value −0.5 or 0.5. The value of μ+ was set either to 0.25μmax or 0.75μmax.

Data availability

The code needed to replicate the results presented here can be found at

Additional information

How to cite this article: Grilli, J. et al. Feasibility and coexistence of large ecological communities. Nat. Commun. 8, 14389 doi: 10.1038/ncomms14389 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Material

Supplementary Information:

Supplementary Figures, Supplementary Tables, Supplementary Notes and Supplementary References


J.G. was funded by the Human Frontier Science Program, S.A. and G.B. were supported by NSF-1148867. We thank J. Hidalgo, D. Logofet, M.J. Michalska-Smith, R. Rohr, S. Saavedra and M. Wilmes for comments.


The authors declare no competing financial interests.

Author contributions J.G., S.S., J.R.B., S.A. and A.M. devised the study. J.G. performed the analysis and drafted the main text. J.G., M.A., S.S., G.B., J.R.B. and A.M. discussed the project, contributed to the derivations and edited the manuscript.


  • Gillooly J. F., Brown J. H., West G. B., Savage V. M. & Charnov E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001). [PubMed]
  • Walther G.-R. et al. . Ecological responses to recent climate change. Nature 416, 389–395 (2002). [PubMed]
  • Tylianakis J. M., Didham R. K., Bascompte J. & Wardle D. A. Global change and species interactions in terrestrial ecosystems. Ecol. Lett. 11, 1351–1363 (2008). [PubMed]
  • Harmon J. P., Moran N. A. & Ives A. R. Species response to environmental change: impacts of food web interactions and evolution. Science 323, 1347–1350 (2009). [PubMed]
  • Tylianakis J. M., Laliberté E., Nielsen A. & Bascompte J. Conservation of species interaction networks. Biol. Conserv. 143, 2270–2279 (2010).
  • Meszéna G., Gyllenberg M., Pásztor L. & Metz J. A. J. Competitive exclusion and limiting similarity: a unified theory. Theor. Popul. Biol. 69, 68–87 (2006). [PubMed]
  • Barabás G., Pásztor L., Meszéna G. & Ostling A. Sensitivity analysis of coexistence in ecological communities: theory and application. Ecol. Lett. 17, 1479–1494 (2014). [PubMed]
  • Barabás G., Pigolotti S., Gyllenberg M., Dieckmann U. & Meszéna G. Continuous coexistence or discrete species? A new review of an old question. Evol. Ecol. Res. 14, 523–554 (2012).
  • Rohr R. P., Saavedra S. & Bascompte J. On the structural stability of mutualistic systems. Science 345, 1253497–1253497 (2014). [PubMed]
  • May R. M. & Leonard W. J. Nonlinear aspects of competition between three species. SIAM J. Appl. Math. 29, 243–253 (1975).
  • May R. M. Will a large complex system be stable? Nature 238, 413–414 (1972). [PubMed]
  • Allesina S. & Tang S. Stability criteria for complex ecosystems. Nature 483, 205–208 (2012). [PubMed]
  • Suweis S., Grilli J. & Maritan A. Disentangling the effect of hybrid interactions and of the constant effort hypothesis on ecological community stability. Oikos 123, 525–532 (2014).
  • Allesina S. et al. . Predicting the stability of large structured food webs. Nat. Commun. 6, 7842 (2015). [PMC free article] [PubMed]
  • MacArthur R. H. Fluctuations of animal populations and a measure of community stability. Ecology 36, 533–536 (1955).
  • Gilpin M. E. Stability of feasible predator-prey systems. Nature 254, 137–139 (1975).
  • Roberts A. The stability of a feasible random ecosystem. Nature 251, 607–608 (1974).
  • Goh B. & Jennings L. Feasibility and stability in randomly assembled Lotka-Volterra models. Ecol. Model. 3, 63–71 (1977).
  • Saavedra S., Rohr R. P., Olesen J. M. & Bascompte J. Nested species interactions promote feasibility over stability during the assembly of a pollinator community. Ecol. Evol. 6, 997–1007 (2016). [PMC free article] [PubMed]
  • Volterra V. Leçons sur la Théorie Mathématique de la Lutte Pour la vie. Bull. Amer. Math. Soc 42, 304–305 (1936).
  • Goh B. S. Global stability in many-species systems. (1977).
  • Logofet D. O. Stronger-than-Lyapunov notions of matrix stability, or how flowers help solve problems in mathematical ecology. Linear Algebra Appl. 398, 75–100 (2005).
  • Johnson C. R. Positive definite matrices. Am. Math. Monthly 77, 259 (1970).
  • Tang S. & Allesina S. Reactivity and stability of large ecosystems. Front. Ecol. Evol. 2, 2–21 (2014).
  • Svirezhev Y. M. & Logofet D. O. The Stability of Biological Communities Nauka (1978).
  • Allesina S. & Tang S. The stability-complexity relationship at age 40: a random matrix perspective. Popul. Ecol. 57, 63–75 (2015).
  • Cohen J. E., Briand F. & Newman C. M. Community Food Webs: Data and Theory Springer (1990).
  • Williams R. J. & Martinez N. D. Simple rules yield complex food webs. Nature 404, 180–183 (2000). [PubMed]
  • Bascompte J., Jordano P., Melián C. J. & Olesen J. M. The nested assembly of plant-animal mutualistic networks. Proc. Natl Acad. Sci. USA 100, 9383–9387 (2003). [PubMed]
  • Bascompte J., Jordano P. & Olesen J. M. Asymmetric coevolutionary networks facilitate biodiversity maintenance. Science 312, 431–433 (2006). [PubMed]
  • Bastolla U. et al. . The architecture of mutualistic networks minimizes competition and increases biodiversity. Nature 458, 1018–1020 (2009). [PubMed]
  • Thébault E. & Fontaine C. Stability of ecological communities and the architecture of mutualistic and trophic networks. Science 329, 853–856 (2010). [PubMed]
  • Armstrong R. & McGehee R. Coexistence of species competing for shared resources. Theor. Popul. Biol. 9, 317–328 (1976). [PubMed]
  • Abrams P. A. The theory of limiting similarity. Annu. Rev. Ecol. Syst. 14, 359–376 (1983).
  • Armstrong R. & McGehee R. Competitive exclusion. Am. Natural. 15, 151–170 (1980).
  • Yodzis P. The indeterminacy of ecological interactions as perceived through perturbation experiments. Ecology 69, 508–515 (1988).
  • Dambacher J. M., Li H. W. & Rossignol P. A. Relevance of community structure in assessing indeterminacy of ecological predictions. Ecology 83, 1372–1385 (2002).
  • Aufderheide H., Rudolf L., Gross T. & Lafferty K. D. How to predict community responses to perturbations in the face of imperfect knowledge and network complexity. Proc. R. Soc. Lond. B 280, 20132355 (2013). [PMC free article] [PubMed]
  • James A., Pitchford J. W. & Plank M. J. Disentangling nestedness from models of ecological complexity. Nature 487, 227–230 (2012). [PubMed]
  • Barabás G. & Allesina S. Predicting global community properties from uncertain estimates of interaction strengths. J. R. Soc. Interface 12, 20150218 (2015). [PMC free article] [PubMed]
  • Allesina S. & Pascual M. Food web models: a plea for groups. Ecol. Lett. 12, 652–662 (2009). [PubMed]
  • Suweis S., Grilli J., Banavar J. R., Allesina S. & Maritan A. Effect of localization on the stability of mutualistic ecological networks. Nat. Commun. 6, 10179 (2015). [PMC free article] [PubMed]
  • Suweis S., Simini F., Banavar J. R. & Maritan A. Emergence of structural and dynamical properties of ecological mutualistic networks. Nature 500, 449–452 (2013). [PubMed]
  • Lever J. J., van Nes E. H., Scheffer M. & Bascompte J. The sudden collapse of pollinator communities. Ecol. Lett. 17, 350–359 (2014). [PubMed]
  • Staniczenko P. P. A., Kopp J. C. & Allesina S. The ghost of nested-ness in ecological networks. Nat. Commun. 4, 1391 (2013). [PubMed]
  • Redheffer R. Volterra multipliers II. SIAM J. Algebr. Discrete Methods 6, 612–623 (1985).
  • Ribando J. M. Measuring solid angles beyond dimension three. Discrete Comput. Geom. 36, 479–487 (2006).
  • Gourion D. & Seeger A. Deterministic and stochastic methods for computing volumetric moduli of convex cones. Comput. Appl. Math. 29, 215–246 (2010).

Articles from Nature Communications are provided here courtesy of Nature Publishing Group