The feasibility of equilibria in large ecosystems: A primary but neglected concept in the complexity-stability debate

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.


Model
Moderate interactions Weak interactions

S.2.1 Feasibility and stability in dynamical systems
Consider dynamics driven by a general system of first order autonomous ordinary differential equations dx i (t) dt = f i (x(t)) , i ∈ S = {1, . . . , S} , t 0, with f = (f 1 , . . . , f S ) : R S → R S a differentiable function and x(t) = (x 1 (t), . . . , x S (t)) ∈ R S . Solutions x(t) of this system evolve in time and their trajectories are rarely describable as such, especially in high dimensions. However, there are some important features that can be exhibited without solving the system. One of the most instructive is the existence of equilibria, which are points x * in the phase space R S where the solution of the system does not vary in time, that is f i (x * ) = 0 for all i ∈ S. A trajectory reaching an equilibrium remains indefinitely at this point. If x * belongs to the admissibility domain of the model (i.e., all x * i > 0 for ecological systems), then it is called feasible. Unfeasible equilibria are well-defined in the mathematical sense, but have to be rejected in the perspective of the application.
If small perturbations of the trajectory in the neighbourhood of the equilibrium x * fade over time, so that the system tends to restore the equilibrium, then x * is said to be locally stable.
This property is related to the derivatives of the function f . More specifically, an equilibrium x * is linearly stable when the Jacobian matrix J(x * ) = ∂f i ∂x j (x * ) evaluated at the equilibrium has only eigenvalues with negative real parts (see [11], [67] or [68]).

S.2.2 Lotka-Volterra model
In the context of ecological networks, the dynamics of interacting species is commonly described by the Lotka-Volterra equations (see e.g. [69,70,29]). The function f can be written as where x i denotes the abundance of species i, θ < 0 is a friction coefficient (intraspecific competition, assumed to be the same for all species), r i describes the intrinsic growth rate of i and the interaction coefficients a ij stand for the per capita effect of species j on species i. The connectance C denotes the proportion of links present in the network with respect to the number of all possible links S(S − 1) or S 2 , depending on the model. The product CS is a measure of the complexity of the system, and is the average number of links between two species. The term (CS) δ is introduced as a normalisation of the interactions strength, with δ 0 a parameter controlling this renormalisation (see 2.3.3 below). In matrix form, the system writes where I denotes the S × S identity matrix, A = (a ij ) 1 i,j S is the interaction matrix, r = (r 1 , . . . , r S ) is the growth rates vector and • denotes the Hadamard product, that is x • y = (x 1 y 1 , . . . , x S y S ).
The admissibility domain of this model is the positive orthant, that is x ∈ R S with x i > 0 for all i. A feasible equilibrium must then satisfy x * i > 0 for all species i ∈ S and so that the right hand side of (S2) equals zero, assuming this inverse matrix exists. Such an equilibrium is locally stable if the Jacobian matrix (also named community matrix in this context) has all eigenvalues with negative real part. As one can see in the previous equation, this matrix depends explicitly on the equilibrium abundances x * (diag (x * ) denotes the matrix with x * in the diagonal and 0 everywhere else) and it should not be confused with the interaction matrix A.
Several ecological models have been developed with the goal of understanding network structure and its effects on system dynamics [71,10,17,18,28]. Those models generate essentially different types of network structures, and consequently particular interaction matrices A. In order to exhibit the particular effects of a structure on the dynamics of systems, the coefficients a ij can be randomised. This allows the exploration of some behaviour of the models and the detection of their key features. A direct consequence in considering random interactions is that feasibility and stability have to be considered from a probabilistic point of view. We define therefore the probability of feasibility of an equilibrium as The purpose of this work is to study this probability, which expresses the likelihood of a model to provide feasible equilibria, for which the stability can eventually be examined.

S.2.3 Interactions
The interaction matrix A describes firstly who interacts with whom in the network (the structure), and secondly what is the type of these interactions: if a ij < 0 and a ji < 0, species i and j compete with each other; if a ij > 0 and a ji > 0, their interaction is mutualistic and they both benefit from the presence of each other; if a ij > 0 and a ji < 0, then species i preys upon species j; if a ij = 0, species j has no direct effect on i.
The models of complex ecological networks that we consider can be divided into two main categories: unstructured and structured models. The former consist in webs for which the topological structure is completely free and random. This equates to considering networks based on Erdős-Rényi graphs, which are constructed by adding randomly edges between nodes with the same probability C. Here, we consider four of those models depending on the type of interactions: random (as in May's formalism [10]), mutualistic, competitive and predatorprey. Structured models define stochastic rules for the construction of the underlying graphs so that not all graphs are equally likely. This is the case in the three models studied here: the cascade [17], the niche [18] and the nested-hierarchy model [28]. Since these three models were built to represent food webs, only predator-prey interactions are considered. The section Method in the Main Text describes their construction in detail. Suppl. Fig. L provides an example of network samples for the unstructured predator-prey model and for the structured cascade and niche model.

S.2.3.1 Interaction models on unstructured networks
Random network. This model introduced by May in 1972 [10] considers both the structure and the type of interactions as fully random. The interactions are independent and identically distributed (i.i.d.) centred random variables with common standard deviation σ and can thus be either positive or negative, corresponding to a mixture of competition, mutualism and predation. The coefficients equal 0 independently of their position, with a fixed probability 1 − C. Note that a ij = 0 does not necessary imply a ji = 0, so that amensalism (0, −) and commensalism (0, +) interactions are also possible. In this model, the connectance is given by C = L/S 2 , where L is the number of links in the network.
Mutualistic network. This model considers only mutualistic interactions (+, +). The diagonal coefficients are a ii = 0 and any pair of species (i, j), i < j is linked with probability C. If the pair (i, j) is not linked, then a ij = a ji = 0, otherwise the interactions strength a ij and a ji are positive i.i.d. random variables. The elements of A are almost all i.i.d., except for the pairs a ij and a ji , which are only independent conditionally on the fact that they are different from 0. In this model, the connectance is given by C = L/S(S − 1).
Competitive network. This model is built in the exact same manner as the previous one, except that the interactions are nonpositive random variables in such a way that there are only competitive interactions (−, −).
Predator-prey network. This is an example of an unstructured model for predation. The diagonal coefficients are a ii = 0 and any pair of species (i, j), i < j is linked with probability C. If the pair (i, j) is not linked, then a ij = a ji = 0, otherwise the interactions strength are independently sampled, with the restriction that sign(a ij ) = − sign(a ji ). This results in a sign antisymmetric interaction matrix whose elements are identically distributed and almost all independent, except for the pairs (a ij , a ji ), which are correlated by their sign. The connectance in this model is given by C = L/S(S − 1).

S.2.3.2 Interaction models on structured networks
Cascade model. This model introduced by Cohen et al. [17] is an example of structured food webs. The interactions are of predator-prey type. The species are ordered on a line and they can feed only on species with a strictly lower rank, excluding any loop in the network. The resulting interaction matrix A has an upper diagonal with nonnegative i.i.d. entries and a lower diagonal with nonpositive i.i.d. entries. An upper diagonal entry a ij , i < j can be zero with probability 1 − C and in this case a ji = a ij = 0. The diagonal of A consists in zeros and the connectance is C = L/S(S − 1).
Niche model. In 2000, Williams and Martinez [18] proposed a new model that permits loops (including cannibalistic loops), and that generates purely interval food webs. Each species is randomly assigned three numbers: a niche value, a range radius proportional to the niche value and a range centre. Species feed on all species whose niche value falls into their range. The species with the smallest niche value has a range 0 to ensure the presence of, at least, one basal species. This defines the so-called adjacency matrix of the network: α with α ij = 1 if i preys upon j and α ij = 0 otherwise. The parameters of the probability distributions of the niche and the range are chosen in order to obtain an average connectance C/2 in the matrix α. For the construction of the corresponding interaction matrix, one multiplies the entries of α with positive i.i.d. values and the entries of α t with negative i.i.d. values, and finally one adds up the two corresponding matrices. This results in a interaction matrix A with upper diagonal with mostly (but not exclusively) nonnegative entries, and lower diagonal entries mostly nonpositive. The average connectance of A is C = L/S 2 . Notice that this construction is slightly different from the one proposed by Allesina and Tang in [15], since if species i preys on j and vice-versa, those interactions do not cancel out but are added in A. The combination of both interactions can either be mutualistic or competitive, or can remain a predator-prey interaction.
Nested-hierarchy model. Cattin et al. [28] proposed a model that tries to implicitly take evolution into account and relaxes the intervality of the diets of the niche model. The model creates a nested hierarchy between species driven by phylogenetic constraints. Two numbers are randomly assigned to each species: a niche value and a number of prey that depends on the niche value. The species with the smallest niche value has no prey. After having reordered the species according to their niche value, prey are assigned to species starting with the species that has the smallest niche value, according to the following rule: for a consumer i, one chooses randomly a prey j with smaller niche value; then, one considers a pool of prey consisting of all prey consumed by other consumers of j; consumer i will then be assigned prey among this pool; if the pool is too small, choose another pool, and if this is still not possible, choose randomly a new prey. This algorithm builds the adjacency matrix of the network α, whose connectance is set on average to C/2 by tuning the parameters of the probability distribution of the niche values and the number of links of each species. The interaction matrix A is then constructed in the exact same way as for the niche model and the connectance is C = L/S 2 .

S.2.3.3 Interaction strength
Intensity of interactions plays a key role for the local stability of networks [30,31]. As already suggested by May's criterion [10], in the framework of the random model, interactions strength should decrease at least at rate 1 √ CS in order to preserve stability when complexity CS increases. This consideration motivates the introduction of the normalising parameter δ 0 in (S2) for the study of ecological feasibility of equilibria. Three regimes naturally emerge: strong interactions (δ = 0), moderate interactions (δ = 1/2) and weak interactions (δ = 1). A useful approach for the interpretation of these different regimes is the weight of a node in the network. Consider for example the random model. As the network is unstructured, all species have on average the same role, and the expected total weight of their interactions is since a ij are i.i.d. by assumption. One sees that for δ = 1, W i is constant and does not depend on the complexity. For δ = 1 2 , the total weight is proportional to the square root of the complexity, and for δ = 0, W i depends linearly on the complexity. From the mathematical point of view, these three regimes correspond to three normalisation of the sum of i.i.d. random variables. Indeed, let X 1 , X 2 , . . . be a sequence of i.i.d. centred random variables. The rescaled sum 1 n δ n i=1 X i converges almost surely to zero if 1 2 < δ 1 (law of large numbers). If δ = 1 2 , the central limit theorem states that 1 n δ n i=1 X i converges in distribution to a normally distributed random variable with mean zero and finite standard deviation. Proposition S.3.2 and Theorem S.3.5 below show the link between these probabilistic regimes and the nature of the equilibrium as the size of the network goes to infinity, depending on the parameter δ. While x * converges to a deterministic value when interactions are weak (δ = 1), it obeys a central limit theorem and converges to a well-defined random variable when the interactions are moderate (δ = 1 2 ). Finally, if 0 δ < 1 2 , the rescaled sum does not converge and its standard deviation explodes. Consequently, the case of strong interactions (δ = 0) is not studied, as the mathematical problem is asymptotically ill-posed.

S.2.4 Growth rates
The parameter r is of particular interest when studying feasibility as illustrated in [26]. For any fixed realisation of the random interaction matrix A, it is always possible to tune the parameter r in order to get any desired equilibrium. Choose any vector u ∈ R S and set as a structural growth rates vector. A glance at (S3) shows immediately that the resulting equilibrium is x * = u. Consider now the all-ones vector 1 in R S . This direction 1 is in a sense the most feasible, since it is as far as possible from the boundaries of the admissibility domain (i.e. the positive orthant). Setting the growth rates vector to the so-called structural vector r(1) leads the equilibrium deterministically to x * = 1. Interestingly, the authors of [20] studied stability using non-random interaction weights. They used r(u) with u = (−d/θ)1, for some positive constant d > 0, to avoid transcritical bifurcations where the abundances of the steady state x * i become negative. They argue that with this particular choice, the only way in which a species can become extinct is via a degenerate bifurcation that makes a positive equilibrium unstable, causing the system to move suddenly to a different equilibrium point.
However, choosing the parameter r in such a way, i.e., contingently to A, is clearly uninformative as it erases the whole structure of the model, and nips the purpose of the randomisation in the bud. The parameter r should thus not be set according to a particular realisation of A. This can be achieved with, for example, the mean structural vector v defined by This vector is deterministic, therefore clearly independent of the randomness of A, and it provides equilibria that are, on average, the most feasible. Choosing this v as growth rates vector allows then to study feasibility in the most favourable conditions. The mean structural vectors of the seven models of interest are given in Table 2.
The analytical results of Sections S.3.1.1 and S.3.2.1 hold for more general growth rates vector than v. They allow to consider random vectors r = (r i ) with i.i.d. entries. In the case of unstructured models, v is indeed a particular case of this type of vectors, as the entries are all identical and deterministic, therefore independent of each other. In the following, we always illustrate our results with the mean structural vector v, unless specified.

S.3 Probability of feasibility
We study analytically the probability of feasibility P S in the framework of the random model. For moderate interactions, P S decreases exponentially with the size of the system. For weak interactions, P S is asymptotically equal to one when the growth rate vector is set to v (see Eq. (S6)). We show that for other choices of growth vector, if the interactions are weak, then P S either equals one or zero, depending on the value of the parameters. We provide estimations of P S by means of Monte Carlo simulations for the other models and observe the same phenomenon. In the case of weak interactions, P S can be estimated analytically for all unstructured models.

Mean interaction
Mean structural vector

Nested hierarchy Simulated
Supplementary Table B Models of interactions and their corresponding mean structural vector v. Note that the form of v is the same for the random, competition and mutualism models. In a competitive system, µ A < 0 and in a mutualistic one, µ A > 0. In all unstructured models (random, competition, mutualism, predation), the mean structural vector has the same direction as the all-ones vector, meaning that the structural growth rates are the same for all species. For the structured models (cascade, niche and nested hierarchy), species do have different roles in the network, while in the other models all species behave, on average, identically. The mean structural vector follows then a particular direction which is no more the all-ones vector. One sees that, for the cascade model, for appropriate choice of the parameters θ, µ A + and µ A − leads to a vector v which has negative entries for species with small indices i (top predators) and positive entries for species with large indices (basal species). This type of growth rate vector is a characteristic of such networks [69,70]. For the niche and the nested hierarchy models, the algorithmic construction of the interaction matrix makes the formula of their respective v very complicated. They are therefore obtained by Monte Carlo simulations. As for the cascade model, negative growth rates are observed.

S.3.1 Moderate interactions
We show that in the random model with moderate interactions, there exists asymptotically almost surely no feasible equilibrium of the system (S1). The lack of structure of this model brings independence of the entries of x * , the biomasses at equilibrium. Moreover, the rate of normalisation δ = 1 2 makes these entries follow asymptotically a normal law. Thus, when the number of species becomes large, there is a smaller and smaller probability that all of them have positive equilibrium, as illustrated in Suppl. Fig. A and Fig. 1, Main Text.
This phenomenon is a key feature of this model and does not depend on the choice of the intrinsic rates r. Indeed, we can show analytically that the probability of persistence at equilibrium goes towards zero for any choice of a random vector r with i.i.d. entries. This framework contains the particular case where the intrinsic rates are all equal, i.e. the same direction as v, the mean structural vector of the random model.
We extend this study by estimating P S in the other unstructured and structured models by simulations. For each model, one samples independent realisations of the interaction matrix A and sets the growth rates vector to the mean structural vector v of the model. The equilibrium is found according to Equation (S3) and P S is estimated by the proportion at which the equilibrium is admissible. The simulation results are analogous to the random model in the sense that P S decreases exponentially towards zero with S. Note that we also performed simulations for the random model, whose results correspond nicely to the analytical prediction (see Suppl. Fig. A  (a)).

S.3.1.1 Analytical results
First, we express the law of x * i for arbitrary large S. Our results ensue from the direct application, with some extensions, of Geman's work [33] on solution of random large systems. They hold under the following assumptions: Assumption S.3.1. In the Lotka-Volterra model (S1), we assume that (i) δ = 1 2 ; (ii) the interaction matrix A = (a ij ) has i.i.d. entries with common mean E(a 11 ) = 0; (iii) the intrinsic growth rates vector r = (r i ) has i.i.d. entries; (iv) r and A are independent; (v) the second moments of the laws of A and r satisfy The assumptions (v) and (vi) on the second and higher moments of A and r are very natural and not restrictive, since they are satisfied by a wide collection of laws (normal, uniform, beta, gamma, lognormal,...) as long as the chosen standard deviation is not too large. Note also that the mean structural vector satisfies these assumptions. The last condition allows us to consider that the solution to (S3) always exists. Proposition S.3.2. Under Assumptions S.3.1, the biomasses at equilibrium of the model (S1) converge in law towards Gaussian random variables: , Moreover, for every fixed 1 k S, the collection (x * 1 , ..., x * k ) has independent entries. Proof. The proof is a generalisation of a result given in [33]. We will assume without loss of generality that C = 1. Indeed, E( 1 √ C a ij ) = 0 and Var( 1 √ C a ij ) = Var(a ij |a ij = 0) so that in the case of moderate interactions the variance of the a ij is preserved when dividing any a ij by √ C. Consider the system (1), First approach the solution of the system by a Neumann progression which is given by As in [33], define which is simply the ith component of A We compute the joint moments of the α r,θ and show that they are the same as those of normal random variables. For m fixed and distinct pairs (k j , n j ) 1 j m we thus compute E m j=1 α r,θ (i j , k j , S) n j . This involves the same combinatorics as the one employed in [33] as mixing terms from (S8) arise. Indeed for asymptotical contributions, each chain a l 1 l 2 ...a l k−1 l k r l k has to be paired exactly with itself. Using moreover independence of the (a ij ) and the (r i ), this gives the joint moments are independent. Returning now to the Neumann progression and by calculating variance and expectation, we find that when S → ∞ .
This result is similar to that obtained in [71, equations (16.12) and (16.21)], where mean, variance, and coefficient of variation have been computed. The normal asymptotical distribution and the fact that the equilibria are asymptotically i.i.d. from the previous Proposition S. 3.2 allow to go one step further and to show that there exists almost surely no feasible equilibrium for species-rich systems of type (1). Theorem S.3.3. Under Assumption S.3.1, the probability that an equilibrium of model (S1) is feasible tends toward zero. That is lim S→∞ P S = 0.
Proof. Let us define P (k) S = P x * j > 0, ∀j = 1, . . . , k . With this notation, P S = P (S) S . The convergence in law in (S7), as well as the independence of (x * i ) imply for all k 1 and where Φ denotes the standard Gaussian cumulative distribution function. Since the probability on the right-hand side is strictly less than 1, this implies that lim sup S→∞ P S = 0.
When considering the mean structural vector of the random model r = v (see Suppl. Table 2), the previous proof allows to approximate the probability that an equilibrium is feasible in the following way: which is represented in Suppl. Fig. A.

S.3.1.2 Simulations
The exponential decrease of P S is not restricted to the random model, for which our analytical result holds. Indeed, we have computed numerically the probability of feasibility for the other models, as well as for the random model (see Fig. 1, Main Text). The growth rates vector is set to the mean structural vector for each model respectively (Table 2). An interaction is non-zero with probability C. We choose a ij ∼ N (0, σ) for any non-zero entries of the interaction matrix A in the random model, i.e. when the sign of the interaction does not matter. In this sense, E(a ij ) = 0 and Var(a ij ) = C · Var(a ij |a ij = 0) = C · σ 2 . In the other cases, a strictly positive interaction is randomly drawn from a folded normal distribution such that a ij ∼ |N (0, σ)|. The expectation is E(a ij |a ij = 0) = σ 2/π and the variance Var(a ij |a ij = 0) = σ 2 (1 − 2/π). A strictly negative interaction is similarly sampled such that a ij ∼ − |N (0, σ)|.
In Fig. 1 of Main Text, the connectance as been fixed to C = 0.25 for each model. To reach such a connectance in the case of the niche model, we choose the niche values uniformly on the interval [0; 1], and their breadth according to a Beta random variable with shape parameters (1, 3). In the case of the niche and nested-hierarchy models, we also simulated the equilibria feasibility when not using the mean structural vector v. In Suppl. Fig. B, the growth rates are i.i.d. gaussian random variables with standard deviation fixed to 0.15 and mean one (in Suppl. Fig. B (a-d)) or mean zero (in Suppl. Fig. B (e-f)). In the former, the probability of feasibility rapidly decreases towards zero, similarly as what is predicted with random models. When E(r i ) = 0 for all i, the equilibria abundances are centered around zero with a positive variance independently of the species index i. As such, there is no chance to observe P S > 0 for this choice of r i .

S.3.2 Weak interactions
For the case δ = 1, we provide the conditions under which the system (S1) possesses almost surely a feasible equilibrium. In this situation, the interactions between the species become extremely weak when the size of the system grows. Compared to the moderate case, the variance of the solution approaches zero allowing us to use a law of large numbers. This drives each component x * i of the solution to a constant proportional to its intrinsic growth rate r i as illustrated in Suppl. Fig. C. This leads to feasible equilibria for any positive growth rate.
Analytical results are given in the case of unstructured models, whereas structured models are explored by mean of simulations.
As for the Assumptions S.3.1, these conditions are not restrictive on the choice of the distributions for the (a ij ). The mean structural vector satisfies the Assumption S.3.4 (iii) for all the models that are considered. In (ii), E(a 11 ) = Cµ A is equivalent to E(a 11 | a 11 = 0) = µ A since C is the probability that an entry of the matrix is set to 0. Note furthermore that, in the settings of weak interactions, the entries of A are divided by CS so that E(a ij /(CS)) = µ A /S. We first introduce the analytical result for the random model. In the case of the random model, setting the growth rates vector to v = |θ| · 1 leads thus to x * ,∞ i = 1 for all i. But Theorem S.3.5 is more general and allows to consider any growth rates vector r with bounded entries. In Suppl. Fig. C for example, the entries of r are chosen to take only two values and one sees that the equilibrium vector also converges only toward two values (proportional to r i ) when S becomes large.
In Assumption S.3.4, the expectation µ A does not need to be zero. This enables the derivation of analytical results in the case of random models in which the mean interaction is positive (mutualistic and commensalistic interactions) or negative (competitive and amensalistic interactions). Indeed, a generalisation of the previous proof by allowing arbitrary µ A leads to the following: Theorem S.3.6. Under Assumption S.3.4, the asymptotic equilibrium of the model (S1) is given by However, the competition and mutualistic models introduced in Section S.2.3.1 do not completely satisfy the Assumptions S.3.4, since some dependences are introduced among the entries of the interaction matrix A. Indeed, a ij = 0 ⇔ a ji = 0, so that commensalistic/amensalistic interactions are forbidden. Moreover, the particular case of predation leads to a sign antisymmetric matrix A, i.e. a ij > 0 ⇔ a ji < 0. In the following, we show that the same type of results are obtained for these cases.
For mutualism or competition, we define the entries of A in the following way, where (w ij ) are i.i.d. random variables of mean µ A and variance σ 2 (e.g. folded normal random variables) and (b ij ) are i.i.d. Bernoulli random variables of parameter C. Note that for every i, j ∈ S, E(a ij ) = E(a ji ) = Cµ A , so that if we define M to be the matrix with all elements set to Cµ A , the matrixŴ is centred and that Cov(ŵ ij ,ŵ ji ) = µ 2 A · C(1 − C) = 0. Like in the proof [33,Thm. 2], we need to show that Ŵ /S → 0. This assertion is based on [72], which gives an upper bound for the norm of sample covariance matrix, and can be related to the work in [73], where the largest eigenvalue of random symmetric matrices is studied. Here we will show that Ŵ /S → 0 in the framework of [74], where it is demonstrated that the limiting distribution of the eigenvalues of a sample covariance matrix 1/SV V T remains the Marcenko-Pastur law, even with some dependencies among the entries of a centred random matrix V . Lemma S.3.7. Consider the random matrixŴ defined by equation (S10) and assume that there exists a constant K such that E(|ŵ ij | k ) K 2k for any k S and 1 i, j S. Then Ŵ /S → 0 almost surely.
Proof. The conditions (MP1), (MP2) and (MP3) in [74] still hold forŴ in our framework. Thus letting λ max be the largest eigenvalue of 1/SŴŴ T , the same combinatorial arguments can be used, and following [74, Section 4] we arrive to where c is a constant and Tr(·) the trace operator. Letting > 0 and using the Markov inequality, we find which goes to zero when S is large and since S P 1 √ S 1 √ SŴ > < ∞ for k 3, almost sure convergence holds. Corollary S.3.8. Under Assumption S.3.4 with a matrix A so that a ij = 0 ⇔ a ji = 0, and where E(a ij |a ij = 0) = µ A and r = γ · 1, with γ ∈ R * , the asymptotic equilibrium of the model (S1) is given by for all i = 1, . . . , S. Under these assumptions and for γ/(|θ| − µ A ) > 0, x * ,∞ is feasible with probability one.
Proof. Since each line in A, i.e. the collection (a ij ) 1 j S for arbitrary i ∈ S, still contains independent elements for arbitrary S, the proof is analog to the one of Thm S.3.6 by using Lemma S.3.7.
In the case of predation, interactions have the form (+, −), so that the previous Corollary extends in the following way.

S.3.2.2 Simulations
Unstructured networks. As in the moderate case, simulations have been performed to illustrate our analytical results. We illustrate the results for the random model in Suppl.  Structured networks. The case of the cascade model [17] is illustrated in Suppl. Fig. D (a) and (b). We represent the envelop of the equilibria for different S from 1000 simulations with C = 0.25, σ = 0.4 and θ = −1. The standard deviation of any x * i decreases towards zero when S increases, independently of the role of species i. Using the mean structural vector, the same convergence to the equilibrium 1 as in the random model is hence obtained. We hypothesise that this result is a consequence of the cascade model yielding adjacency matrices with inferior triangular parts constructed like a Erdős-Rényi network. The distribution of a top, an intermediate and a basal species are illustrated in Suppl. Fig. E (a) and shows that any species abundance is independent of its role in the web, exactly as in random models.
For the niche model [18] and the nested-hierarchy model [28], a particular phenomenon occurs that is not observed for the other models. Indeed, simulations show that the equilibrium remains random when S increases. This result is likely attributable to the construction of the models, with the standard deviation of the abundances depending on the hierarchical position of the species. Importantly, standard deviations do not decrease towards zero as the size of the network grows, as illustrated in Suppl. Fig. D (d) and (f). Indeed, these highly structured networks and their related mean structural growth rates induce correlations among the abundances at equilibrium that prevent a convergence to a deterministic value, but rather takes place on a compact support. This is illustrated by the envelop of the equilibria that is represented in Suppl. Fig. D (c) and (e) for different values of S. The empirical distribution of x * i for particular species roles i is illustrated in Suppl. Fig. E. From these results, it is apparent that basal species may rapidly converge to a deterministic constant. However, for intermediate and top species, simulations show that the convergence is very likely to occur on a non-trivial support.
We also tested i.i.d. Gaussian random growth rates with σ r = Var(r i ) = 0.15 for the niche and the nested-hierarchy models in Suppl.

S.3.3 Empirical networks
We used our method on empirical food webs to illustrate how species abundances at equilibrium depend on their intrinsic role in the web. The topology is thus fixed to the observed one and the interactions are modeled as for any random structured web (see the section Method in the Main Text). We report the results of four arbitrary networks. In Suppl.

S.3.4 Connectance and interactions strengths
In our approach, we fixed the connectance C to arbitrary values in the mathematical developments, and to fixed values in the simulations. In this section, we briefly explore the consequence of a relationship between connectance and species number on feasibility, as observed in several contributions, e.g. [75,42]. Observe first that, from the point of view of the average interaction strength, introducing zeros with a given probability 1 − C in the interaction matrix is the same as multiplying the whole matrix by C. Therefore, one has where B is a matrix with all entries i.i.d. and such that P(b ij = 0) = 0. Consider C = 1 S β , for β 0, which is a flexible function that adequately captures observed relationships between C and S. We get BC In the Main Text and above, we showed that feasibility is warranted when δ + β(1 − δ) > 1 2 , i.e. for Taking β 1 2 leads almost surely to feasible equilibria in the models considered here, independently of the exponent δ. For β < 1 2 , δ can be smaller than 0.5. This shifts the three previously described regimes for δ to the left, so that δ can be smaller to reach the same results. Consequently, the existence of a relationship between C and S does not affect our conclusions.

S.4 Consequences on local stability analysis
Any local stability analysis should be preceded by a feasibility analysis [27,22,76]. Indeed, studying the local stability of a point that is not biologically feasible is not instructive on the behaviour of the network dynamics. Moreover, the Jacobian matrix of the system (S4), therefore its eigenvalues too, depend explicitly on x * . Therefore, there is no warranty that randomly sampling directly the Jacobian may be sufficient to address local stability. Here, we explore this question, which extends the analyses of May and Allesina and co-workers [10,15].

S.4.1 Random networks
We first focus on the random model of May [10] in the case of moderate interactions in the Proposition S.4.1 below. In this case, we find that May's criterion still holds when we randomly sample the interaction matrix A (and not directly the Jacobian) under the additional condition that the equilibrium x * is feasible.

S.4.1.1 Moderate interactions
In Proposition S.4.1, we first show that May's criterion for stability, i.e. √ CSσ < |θ|, is still sufficient for any feasible equilibrium whenσ is the standard deviation of the normalised interaction matrix A/(CS) δ , and not the standard deviation of the Jacobian matrix of the system. Recall that the case of moderate interactions is here very natural, as relying on Wigner's [12] and Girko's [77] original convergence results. The criterion simply becomes Proof. Consider A/ √ CS so that σ < |θ|. This is equivalent to say that B = θI + 1 √ CS A possesses only eigenvalues with negative real parts. Define now the symmetric matrixB = DB + B T D, where D = diag 1 √ 2 . Each non-diagonal entry of the matrixB are independent of the others up to symmetry, and have mean zero and standard deviation σ. By the Wigner semicircle law and May's criterion, its eigenvalues are almost surely all negative when S → ∞. The matrix B is thus dissipative when S → ∞, which implies that any feasible equilibrium of the model (S1) is locally stable (see [76]).
Note that feasibility is required to conclude on the stability of the equilibrium. Indeed, we illustrate in Suppl. Fig. J that a stable equilibrium in the sense of J = (θI + A/ √ CS) (as in [10], i.e. so that √ SCσ < |θ|), is stable in the sense of J = diag(x * )(θI + A/ √ CS) only when it is feasible, which we showed is never the case in large systems (see Thm. S.3.3).

S.4.1.2 Weak interactions
In the case of weak interactions, May's criterion [10] and the criteria for competition and predation-prey established by Allesina and Tang [15] are asymptotically trivially satisfied if the parameters are chosen so that the (deterministic) equilibrium is admissible. Indeed, the criterion becomes σ/ √ SC < |θ| and σ/ √ SC → 0 for S → ∞. Concerning mutualistic networks, the criterion in [15] is non-trivial. We show that this criterion still holds under the additional condition that x * is feasible: Under weak interactions (Assumptions S.3.4) and for a mutualistic interaction matrix A, x * is locally stable when x * is feasible and when E(|a ij |) < |θ|. (S11) Proof. the inequality (S11) has been computed in [15] using the Geršgorin circles [78]. With analogous arguments as in the proof of Prop S.4.1, we define the symmetric matrixB = DB + B T D, where D = diag 1 √ 2 and B = θI + 1 CS A . Its largest eigenvalue has been computed in [73] and is given by √ 2E(|a ij |)+ √ 2θ, which is negative when (S11) holds. B is thus dissipative and consequently the equilibrium locally stable [76].

S.4.2 Structured models
Determining stability criteria in structured models is mathematically difficult, since the entries of the resulting random matrix J(x * ) can become highly dependent. In this situation, a direct use of the classical results in [33,79,80] is no longer possible. In Suppl. is represented in yellow and the eigenvalue with largest real part is highlighted in red. With this renormalisation, the convergence result of Girko [77] is recovered and, since the parameters are σ = 0.4, C = 1 and θ = −1, May's criterion holds. However, the equilibrium is very unlikely to be feasible with S = 1000 (see Fig. 1, Main text), so that the criterion of Prop. S.4.1 is violated and the resulting equilibrium is not stable. This is illustrated by the spectrum of the Jacobian matrix J(x * ) in cyan. The eigenvalue with largest real part is highlighted in blue. some spectra of A and J(x * ) for 50 networks of the three structured models to evidence their difference. We find that a link between feasibility and stability is very likely to exist, in a similar way as for unstructured networks.
Supplementary Figure K. Spectra of the Jacobian matrix for different values of δ in structured models. The eigenvalues of J(x * ) are plotted in blue when x * is feasible, and in red otherwise; the eigenvalues of the interactions matrix (θI + A/ √ CS) is represented in yellow. (a -c) For strong interactions, the systems are degenerate for the three models, i.e. never feasible nor stable. (d -i) For moderate interactions and weak interactions, we observe the same kind of relationship between feasibility and stability as for unstructured models. For the cascade model (d, g), x * converges to a deterministic constant; for the niche and nested-hierarchy models (e, f, h, i), x * converges on a compact support, which can include negative values. It is thus necessary to tune the parameters to obtain feasibility; here, we chose S = 150, σ = 0.4, C = 0.25 and θ = −1 in all cases.

S.5 Feasibilty with a Holling type II functional response
The question of feasibility is of course not restricted to the models studied in the present work. Instead of the Lotka-Volterra dynamics which was used in the previous and for which the interactions follows a mass action law (Holling type I [81]), we could for example explore the same question with different type of functional response.
As an illustration, we consider the following system of differential equations [82] in the predator-prey model of Section S.2.3.1, The reason for this difference with the Holling type I models comes from the saturation term. In this mean-field framework, Holling type II dynamics acts has a buffer of interactions by taking into account the number of preys and predators in the system. Asymptotically, it makes interactions weak, whatever the value of the parameter δ. Hence, it ensures the feasibility of the equilibrium.
To understand this, consider the global effect of any species j ∈ S on i given by the generic term b ij = a ij /(CS) δ h + k∈Prey(i) a ik /(CS) δ x k x j .
To simplify, assume that each a ij is equal to its conditional mean a + 0, if i is a predator or a − 0 if i is a prey. Then, considering that, in this predator-prey model, the number of preys and predators of species i is a contant time CS, that is |Prey(i)| = |Pred(i)| = O(CS) in Landau notation, b ij can be approximated by where a ± = a + if j is a predator and a ± = a − otherwise,x is the empirical mean of all considered x j , γ andγ are some positive constant. Note that the last approximation holds for 0 δ 1 and S large. Therefore, the order of an interaction b ij is 1/CS (weak interaction) no matter what regime δ was initially set. The system considered is the random predator-prey model. The 95%-confidence intervals are calculated over 1000 simulations. Holling type II reacts as if the interactions were weak in a mass-action dynamics. In particular, the probability of feasibility approaches one in this case. For comparison purposes, we set r = 1, σ = 0.4, C = 0.4, and θ = −1 in both cases and in these simulations, h = 1. The system at equilibrium becomes non-linear but has been solved numerically with the method fsolve in Matlab R2012b.