Cooperation promotes biodiversity and stability in a model ecosystem

# Cooperation promotes biodiversity and stability in a model ecosystem

###### Abstract

Empirical observations show that natural communities of species with mutualistic interactions such microbes, can have a number of coexisting species of the order of hundreds. However standard models in population dynamics predict that ecosystem stability should decrease as the number of species in the community increases and that cooperative systems are less stable than communities with only competitive and/or exploitative interactions. Here we propose a stochastic model which is appropriate for species communities with mutualistic/commensalistic and exploitative interactions and that can be exactly solved at the leading order in the system size. We obtain results for relevant macro-ecological patterns, such as species abundance distributions and correlation functions. We also find that, in the large system size limit, any number of species can coexist for a very general class of interaction networks. For pure mutualistic/commensalistic interactions we analytically find the topological properties of the network that guarantee species coexistence. We also show that the stationary state is globally stable and its stability increases as the number of species grows, as empirical observations suggest.

Keywords: Voter Model, Population Dynamics, Ecological Networks, Cooperation, Stability-Complexity

Research in population dynamics has a long history dating back to almost one thousand year ago with Fibonacci modeling of rabbits population. Nevertheless it is still unknown which are the mechanisms allowing the coexistence of many interacting species in the same environment [47, 52, 45, 29]. The current loss of earth biodiversity [46] makes this open question of great relevance today more than ever, and physicists may decisively contribute to tackle this challenge [15, 18]. Historically, the Lotka and Volterra (LV) equations [51, 72] have provided much theoretical guidance and several microscopic derivation of these equations have been proposed [14, 55, 7]. Furthermore, these equations are the core of most of the multi-species deterministic population dynamics models based on the ecological concept of niche partitioning: competing species in order to coexist need to interact with the environment differently and to rely on not-overlapping resources [47, 52].

While prey-predator and competitive interactions have been extensively studied [40, 19, 55], mutualistic/commensalistic interactions, which are beneficial to one or both the involved species, have historically received less attention. The current approach to mutualistic population dynamics is a mere generalisation of the LV types of models, which does not change the functional form of the phenomenological equations with a two-body interaction, but utilizes beneficial (+ +) instead of exploitative (+ –) interactions among individuals of different species [43, 52, 36]. In particular, a microscopic derivations of the phenomenological equations for the population dynamics in mutualistic communities is still missing. Moreover, a generalization of the stability-complexity theorem [52, 54] has revealed that mutualism is even more detrimental to stability as the product increases [16, 68, 66, 29], where is the number of species and , the connectivity, is the fraction of non-zero pairwise interactions between species. This prediction clashes with the observation of widespread mutualistic interactions (or other facilitating interactions) in many natural communities where the biodiversity is very high [20, 13].

An alternative theoretical approach to niche-based multi-species deterministic modeling is the Neutral Theory (NT) of Biodiversity [45, 69, 70, 17, 57, 22, 18, 44]. In NT organisms of a community have identical per-capita probabilities of giving birth, dying, migrating, and speciating, regardless of the species they belong to. In this sense NT is symmetric and aims to model only species on the same trophic level - species therefore competing for the same pool of resources. For instance, plants and trees in a forest compete for resources like carbon, nitrate and light. An important example of neutral model is the Voter model (VM) [42, 30, 65]. In its simplest version, at every time step a randomly selected individual dies and the corresponding resources are freed up to be colonized by an individual randomly sampled within the community [22, 18]. An important limitation of this modeling is that it does not consider explicitly species interactions (e.g. mutualism/commensalism).

Although it has been already shown that niche based and neutral approaches are only apparently contrasting [10, 11], two crucial issues in the current literature are: (i) the lack of a general framework specifically developed to model mutualistic and commensalistic interactions where species interactions are added on neutral models and can modify the species birth-death rates; and (ii) understand the role of cooperation in determining species coexistence and ecosystem stability[15, 18]. In this Letter we present a theoretical framework, where starting from a VM, we add interactions among species and we properly account the effect of cooperation. These interactions affect neutrality and lead, in their mean field formulation, to an alternative to LV models. Reconciling apparently contrasting observations and previous results [52, 54, 16, 68, 29], we show that in our model ecosystem cooperation promotes biodiversity and diversity increases its stability.

In details, be the species label at spatial position , where and . The state at time of the system is given by . We also set to be the fraction of individuals of the -species. We now introduce a directed graph on the set , where the nodes correspond to species and directed links represent the network of ecological interactions. Such a graph is defined through two matrices (cooperation matrix) and (exploitation matrix) satisfying the following conditions: for all , ; For all it must be or ; For all we have . Both intra and inter-species competition is indirectly accounted by fixing the total number of individuals in the community [45, 18].

In ecological terms, given two species and , a directed link of strength from to means that the -th species receives a beneficial effect from the interaction with the -th species, while and denotes that the -th species exploits (is exploited by) the -th species. Just as an example: in the former case we can think to a microbial community where the presence of a certain species creates an environment which favors the growth of other bacteria [13, 9]; in the latter one, we may think to host-parasite symbiosis. Typically it is very difficult to measure the strength of the interactions among two species, but adopting a standard approach [16, 41], we draw the matrix entries from a given probability distribution.

We deal with a well mixed system where spatial effects can be neglected and the dynamics is described by a continuum time stochastic Markov process: a randomly chosen individual is removed and substituted by an individual of the -th species at a rate

 ω(j,η,M,L)=¯ηj+ϵ1∑k¯ηkMkjθ(¯ηj)+ϵ2∑k¯ηkLkj¯ηj, (1)

where and give the cooperation and exploitation intensity, and is the Heaviside step function, i.e., when and 0 otherwise. The presence of the -function in the mutualistic contribution, guarantees that the transition rate is zero if the -th species is extinct. For we recover the standard VM. When the species is favored by the presence of the other species ( in the summation) to which it is connected and by their population; on the other hand allows the possibility that a species exploits (or is exploited by) one or more other species.

It is important to highlight the differences of the contribution on eq. (18) between exploitative and mutualistic/commensalistic interactions. In the first case, the interaction term is quadratic in , (i.e. ) as exploitative interactions can be derived using the law of mass-action used to describe chemical reactions [14, 55, 7]: species must physically get in contact and the chance of interaction is, in the simplest hypothesis, proportional to both species concentrations. On the other hand, in mutualistic/commensalistic relationships, the contribution to the birth rate is linear in , (i.e. ). Indeed, mutualistic interactions are typically mediated by some resource for the species produced by the species (e.g. pollen, faecal pellets, metabolic waste) and we assume that these resources are always fully utilized in the community. In this setting, the benefit that a species receives does not depend on its own abundance (that is limited by that resource), but only on the abundance of the mutualistic partner. In the Supplementary Material, Section 1, we derive this result from a model taking explicitly into account resource dynamics.

The microscopic dynamics given by rates (18) induce a Markovian evolution on the relative abundance of each species. Standard techniques [32] can be used to prove that as , the process weakly converges to the solution of the system of ordinary differential (mean field) equations:

 ddt¯ηs(t) =ϵ1S∑k=1¯ηk(t)Mksθ(¯ηs(t))+ϵ2S∑k=1¯ηk(t)Lks¯ηs(t) −¯ηs(t)S∑i,k=1(ϵ1¯ηk(t)Mkiθ(¯ηi(t))+ϵ2¯ηk(t)Lki¯ηi(t)) (2)

for , where and is conserved by the dynamics.

All presented results do not change when the hard constraint of total fixed population size is relaxed by introducing the possibility for a site to become empty (see Supplementary Material, Section 2.1). We will show below that, under suitable hypothesis, a stationary solution of eq.(9) exists and it will be denoted .

Through eq. (9) we can study many ecosystem properties of interest. One important emergent pattern in ecology, which we can determine within our model, is the relative species abundance (RSA) [45, 69, 70, 18]. It describes commonness and rarity of species, thus characterizing the biodiversity of an ecological community. In our model, the RSA is exactly given by the mean field stationary solution , which in turn depends on the species interaction matrix and . The cumulative RSA is defined as the fraction of species with population greater that a certain value, , where we have fixed when all species coexist. We numerically find that the shape of the stationary RSA resembles the one found in many real ecosystems [18]. It weakly depends on the distribution from which the interaction strengths (the elements of and ) are drawn, but only on its coefficient of variation (see Supplementary Material Section 3.5). This is important as it allows to reduce the model free parameters.

Another relevant quantity characterizing the ecosystem biodiversity is the covariance matrix , , describing the correlations in the population abundance fluctuations between pairs of species population abundances [71]. In our setting we can compute analytically this quantity in the limit of normal fluctuations. We define species abundance fluctuations as for The stochastic process converges in distribution to a Gaussian Markov process which solves the stochastic differential equation where is a -dimensional Brownian motion, which corresponds to a -dimensional Ornstein-Uhlenbeck process [32, 37]. The analytical expressions for the matrices and in terms of the interaction matrices and , and of the equilibria, , of eq. (9), are given in the Supplementary Material, Section 3.3. The covariance matrix, , can be obtained by solving the following Lyapunov matrix equation .

This quantity is typically measured from species population time series, through the Pearson (or other type of) correlations [33]. Moreover, in many studies once opportunely thresholded, it is used as an empirical proxy of the species interactions matrix [33, 50]. In other words many works assume that can be approximated through . Other works, applying maximum entropy approach, use as the quantity to describe the species interactions network [71]. However we find that both and are not good proxies of the species interactions matrix (see Supplementary Material, Section 3.3). This result highlights the importance to properly infer interaction networks from data [33, 8] by considering a suitable model, which explicitly takes into account species interactions.

We now show how our shift in the assumptions behind mutualistic/commensalistic species interactions could resolve the problematic aspect of stability in ecosystem dynamics. In particular, for (voter model with cooperation and indirect competition, but no exploitation) we are able to analytically relate key dynamical features of eq. (9) to the topology of the interaction matrix and prove various results of ecological importance.

First we show that the presence of non-supported species – the i-th species is non-supported if – inhibits coexistence equilibria of the whole ecological community. More precisely if species is non-supported by other species then at stationarity eq. (9) implies that . The extinction of the -th may create new unsupported species that go to zero in the large time limit. Such a cascade of extinctions may eventually end only when for all nodes of the network (see Supplementary Material, Section 2.2). The elimination of nodes of the interaction network corresponding to all non-supported species will be called pruning in the following. Furthermore we have found sufficient conditions on the topology of the mutualistic interaction matrix for the existence of stable stationary states of eq. (9) . In fact, if is irreducible, i.e. if for any node we can reach any other node through a path of oriented links such that , then the Perron-Frobenius (PF) theorem holds [63] and a non-trivial stationary state, , exists with all positive entries and it is unique. It is proportional to the left eigenvector, , of corresponding to the eigenvalue of with the largest modulus, which turns out to be non-degenerate, real and positive, denoted by in the following (for brevity we will refer to it as PF eigenvalue in the following). The corresponding right eigenvector will be denoted by . All components of both and are strictly positive and so . An example of irreducible matrix occurs when implies and the network has a single connected component.

If PF holds and the initial condition , then the time dependent solution of eq. (9) is

 ¯η(t)=¯η(0)TeϵMt∑i(¯η(0)TeϵMt)i . (3)

Since for any eigenvalue, , of we have the dominant term in both numerator and denominator in eq. (14) is leading to . This is an easy computation when has a basis of eigenvectors and in general can be derived using the Jordan decomposition. As a corollary of the derivation above we have also that the stationary solution is globally stable in the region , for all . Therefore, within our framework, we can analytically study the impact of the species interaction network architecture on system stability and species extinction. The results of the mean field predictions are shown in Figure 1. Two simple examples are shown corresponding to an ecosystem with no extinction (panels A-B) and with extinction (panel C-D). Our results can be applied to study the effect of the interaction network topology to species coexistence in real mutualistic ecological communities. In particular, we found that nested architecture [21], observed in plant - pollinators ecological communities [20, 68], where specialist species, with only few mutualistic links, tend to interact with a proper subset of the many mutualistic partners of any of the generalist species, (see Figure 1 E) satisfies the hypothesis of the PF theorem and thus favour species coexistence (Figure 1 F). We have also numerically explored the effect of adding exploitation, i.e. . Specifically we find that adding exploitations does not change the main conclusions of our results, as long as a mutualistic network of interactions is present, corresponding to an irreducible matrix, , and the transition rates given by (18) never become negative during the time evolution of the mean field equation (9) (otherwise it would invalidate the derivation of the equation itself, see Supplementary Material, Section 3.2).

More generally, we can study analytically the stability of the equilibria as a function of ecological complexity, by analyzing the eigenvalues of the linearization of eq. (9), i.e. the Jacobian matrix , around the equilibria, , of the system. We set equal to zero the diagonal of whereas the the off-diagonal pair is equal to with probability and with probability it is drawn from a bivariate gaussian distribution of means and interaction covariance matrix . This guarantees that, for a connected cluster, coexistence of all species occurs. We define , and as mean, variance and correlation of the elements of matrix . The case in which each element of is assigned independently of simply correspond to the case (notice that even if we can have ). Similarly, when considering also exploitative interactions, we can sample randomly the off-diagonal pairs , obtaining a given mean , variance and correlation . If , the leading eigenvalue and the corresponding eigenvector has positive components [16]. Moreover, the components of the leading eigenvector are approximately constant, i.e. the equilibria of system given by eq. (9) can be written as for with . Using the fact that , and taking into account that all the terms involving are sub-leading in , we obtain that the leading term of the system Jacobian does not depend on (see Supplementary Material, Section 3.4) and it is equal to:

 Aij=−δijSμM+(Mij−μM)=−δijSμM+M′ij , (4)

where is a random matrix with zero mean variance and correlation . This implies that the eigenvalues are uniformly distributed in an ellipse centered around with semi-axis and [39, 64]. The largest eigenvalue of the Jacobian is therefore given by . Thus, for fixed connectivity, , the system stability increases with , whereas becomes independent of if the connectivity scale as , that is just above the percolation threshold of a random network [12]. This also implies that in our framework, even moderate mutualistic interactions can stabilize the dynamics and if present the stability increases with the ecosystem complexity (see Figure 2).

Therefore, in the proposed approach cooperation promotes ecosystem biodiversity, that in turn increases its stability without any fine tuning of the species interaction strengths or of the self-interactions [67].

Our framework proposes an alternative approach to model cooperation in species population dynamics starting from an individual based neutral stochastic model. We have developed a generalisation of the classic voter model, adding the effect of species interactions on birth rates. The mean field approximation give us effective equation to model multi-species population dynamics. In particular we show that a shift in the assumptions behind species interactions resolve long-standing open theoretical question on the relation between stability and complexity and provides a unifying modeling approach between neutral stochastic individual based model, useful to describe emergent patterns in ecology and interacting deterministic systems. In particular we highlight that, when properly accounted in the dynamics, mutualistic relationship are crucial in order to have coexistence of species in the communities, as observed very recently in real microbial communities [9].

## 1 Mechanistic interpretation of linear growth rates

The mutualistic dynamics introduced in the main text assumes that the benefit that a species receives from other species is independent of its own abundance. This assumption is radically different from the typical form of growth rates for exploitative (e.g. predator-prey) interactions, where some sort of mass-action law, typical of chemical reaction, is usually invoked . Here we provide a simple mechanistic interpretation on how the form we consider could emerge. While we explicitly consider the case of microbial communities, this could be extended to other systems, too. Here we consider mutualism/commensalism as the presence of certain species is able to create an environment (e.g. by producing some public good or nutrient) or to release some substances (e.g. faecal pellets, metabolic waste), which favor the growth of some others species.

Let be the concentration of a given resource used by the species . This resource is provided, at a rate , by certain species (e.g. through metabolic waste [9, 5]) and related to their populations in a linear way, that is . The kinetic of nutrient concentration is then [3, 5]

 dc(t)dt=1τR(s(t)−ηj(t)r(c(t))), (5)

where is the consumption rate per individual whose specific form is irrelevant for the purpose of this example (e.g., one can consider the Monod function , with and some suitable constants). The constant is the timescale of the dynamics of resources.

The contribution to the growth rate of the -th species, due to this nutrient, is

 δω=ϵr(c)ηj , (6)

where is a conversion factor measuring how the nutrient contributes to the biomass of the -th species. If the nutrient concentration is in quasi-steady state [3], that is , which occurs if it relaxes much faster than populations (i.e. for small ), then, from the above equation, we get

 r(c(t))=s(t)ηj(t). (7)

 δω=ϵs=ϵ∑kηkMkj (8)

if . This is what we have assumed in eq. (1) of our main text when considering mutualistic interactions.

## 2 Proofs and further results for mutualistic/commensalistic ecological communities.

### 2.1 Mean field analysis for the voter model with empty sites.

If we turn off exploitation (), the mean field equation without empty site reads ()

 ddt¯ηs=ϵS∑k=1¯ηkMksθ(¯ηs)−ϵ¯ηsS∑i,j=1¯ηiMijθ(¯ηj) (9)

where represents different species, is the average fraction of individuals of the -th species, is the interaction matrix whose non-zero entries define the network of ecological interactions, is the Heaviside step function ( for ) and is the cooperation intensity (the average of the non-zero is fixed to ). For simplicity we have omitted time dependence of .
An intuitive derivation is as follows. The key point is that for large the evolution of the quantity becomes deterministic because the noise is cancelled in the macroscopic regime and in the thermodynamics limit the relative abundance converges to its mean. Then, observe that the dynamics of the relative abundance in the infinitesimal time is simple as it can only decrease by when a site of kind change type or can increase by when the new symbol of a certain site is .

We now extend the model presented in the main text introducing the possibility for a site to be empty. In our setting empty sites do not interact with species. Thus the species rates remain unchanged after the introduction of empty sites. Thus the species rates are the same as before whereas non-empty sites become empty with rate . In the case , the rate has to be less than 1 otherwise empty sites will cover all the available space. The mean field equations become now:

 ddt¯ηs =¯ηs¯η0−¯ηsλ+ϵS∑k=1¯ηkMksθ(¯ηs)−ϵ¯ηsS∑i,j=1¯ηiMijθ(¯ηj) (10) ddt¯η0 =(1−¯η0)(λ−¯η0)−ϵ¯η0S∑i,j=1¯ηiMijθ(¯ηj). (11)

Let us analyze the stationary mean-field equations for . In this case the stable equilibrium for the empty sites is . Substituting in the equations for , we obtain

 ϵS∑k=1¯ηkMksθ(¯ηs)−ϵ(11−λ)¯ηsS∑i,j=1¯ηiMijθ(¯ηj)+O(ϵ2)=0 (12)

where . After the change of variable , the above eq. (12) reduces to the same equation as one would get for , i.e. in absence of empty sites. In other words, when is small, the introduction of empty sites leads to stationary abundances which are trivially rescaled with respect to the case in absence of empty sites, as a consequence of the reduction of the available space.

### 2.2 Topology of the Interaction Networks, Coexistence and Stationary States

In this section, we discuss some features of the topology of the mutualistic interaction matrix and how they relate to stationary states of the system. The main concept in this section is the one of pruned graph and the operation of pruning a network. A node with in-degree equals to zero and out-degree different from zero is called a dead leaf of the network. The operation of pruning consists in eliminating one by one the dead leaves of a given network together with their outbound links. After a first pruning, we will obtain a new network (that is a subnetwork of the starting one) that may still have dead leaves - the elimination of dead leaves may create new dead leaves. The pruning process end when the resulting network has no more dead leaves. The latter network is called stable or pruned. It is easy to see that the minimal pruned network (i.e. with the smallest number of links) that can be constructed with nodes is the cyclic graph. More in general, we have:

Proposition: The pruned network is a union of isolated nodes and graphs that contain at least one cycle each.

Indeed, pruning stops when the obtained graph is a union of isolated nodes and graphs where all nodes have at least an ancestor (i.e. the in-degree of each node is positive). Now a finite graph where each node has a least one incoming link contains at least a cycle. In fact, starting from one node it is possible to walk through the ancestors and never stop. Since the graph is finite, soon or later, the walker will visit twice the same node - so the walk contains a cycle - at most after a number of steps that equals the size of the graph.

The pruned network has at least one cycle but when not simply union of isolated cycles it can be very complex. Fig. 3 shows an example of the pruning procedure and of a non-trivial pruned network.

As we anticipated at the beginning of this section, the dynamics of species sitting on dead leaves of the interaction network is trivial as their relative abundance goes to zero. This is a simple consequence of the fact that a dead leaf has no incoming bond. Thus, when is a dead leaf, the first term on the right of (9) is zero and simple estimate gives . The previous simple remark leads to the following:

Limiting dynamics of dead leaves: Start the dynamics from a point with for all . If is a dead leaf then .

Thus the presence of a dead leaf inhibits coexistence equilibria on the whole graph. More precisely, if are dead leaves (at some step of the pruning), the stable equilibria must have .

### 2.3 Stability of the equilibria of networks with M irreducible.

Let with the stationary solution of equation (9). If all the components are positive, they are solutions of

 ∑kmkMki=mi∑jkmkMkj . (13)

We consider now the case where the interaction matrix, , is irreducible. This means to require that a path of oriented links (a link is present from to is ) exists joining each pairs of nodes, say and . In this case the interaction network is pruned. Furthermore the Perron-Frobenius theorem holds (see Theorem 1.5 in [4]) implying that a positive solution, i.e. for all , exists. This solution is unique and it is proportional to the left eigenvector, , of corresponding to the eigenvalue of with the largest real part (see Theorem 1.7 in [4]), which turns out to be non-degenerate, positive and will be denoted in the following. All components of are strictly positive and so . Notice that when an irreducible matrix is also aperiodic (i.e. returns to any state can occur in any number of steps) then it is primitive (see Theorem 1.4 in [4]), that is there exists a positive integer such that for all pairs of nodes and . This condition allows for a a stronger version of the Perron-Frobenius theorem (compare Theorem 1.1 with Theorem 1.5 [4]).

We now show that the stationary solution is reached at infinite time as far as the initial condition . It is immediate to see that the time dependent solution of eq. (9) is given by

 ¯η(t)=¯η(0)TeMt∑i(¯η(0)TeMt)i . (14)

If has a basis of eigenvectors then one can expand using the left eigenvectors. Since for any eigenvalue, , of we have the dominant term in both numerator and denominator in eq. (14) is , where is the right eigenvector of corresponding to , leading to

 limt→∞¯η(t)=v∑ivi=m . (15)

As a byproduct of the derivation above we have also that the stationary solution is stable in the whole domain where initial abundances are strictly positive.

Now we want to quantify the degree of stability of the stationary solution, , for the case of random matrices [52, 16]. This is done by analyzing the eigenvalues of the Jacobian matrix resulting from the linearization of eq.((9)) around . This derivation assumes the existence of a feasible solution (i.e., with only positive components) of equation ((13)) exists, but does not require that the hypothesis of Perron-Frobenius theorem hold. This argument can therefore be applied also when some entries of the interaction matrix are negative, including therefore non-mutualistic interactions (e.g., predator-prey), as long as a feasible fixed-point exists.

Standard theory of dynamical systems assures that when the real parts of the eigenvalues are all negative then equilibria are stable. The Jacobian of (9), evaluated at the fixed point, reads

 Aij =Mji−miS∑k=1Mjk−δij∑klMlkml =Mji−miS∑k=1Mjk−δijα . (16)

We set equal to zero the diagonal of whereas the the off-diagonal pair is equal to with probability and with probability it is drawn from a bivariate gaussian distribution of means and covariance matrix :

 Σ=(σ2ρσ2ρσ2σ2)

The global mean, standard deviation and correlation coefficient can be straightforwardly obtained [67]

 μM=CMμ,
 σM=√CM(σ2+(1−CM)μ2)

and

 ρM=ρσ2+(1−CM)μ2σ2+(1−CM)μ2.

We want to prove the following result. The average eigenvalue distribution of the matrix , in the large limit, is uniform in an ellipse in the complex plane, centered at , with real semi-axis and imaginary semi-axis . The last part of the proof is heuristic. First we show that if is an eigenvalue of with , then is an eigenvalue of , where, as in the main text, , i.e. is the left eigenvalue of with eigenvalue , whose existence is assumed (guaranteed if is irreducible as seen above). Indeed if is a right eigenvector of corresponding to the eigenvalue of (i.e. ) then and . Thus the spectrum of is simply the one of shifted by apart from the eigenvalue itself, which is transformed in since .
Let us now define the new matrix . All the eigenvalues are also eigenvalues of whereas the eigenvalue itself is transformed in . In fact, with the same notation as above, it is trivial to see that and . Furthermore the ensemble average of matrix elements . In summary, apart for the eigenvalue , the other eigenvalues of are those of translated by . Thus we have reduced the problem of calculating the spectrum of . For large the law of large number implies that apart from corrections that are gaussian distributed of zero average and variance of order . This implies that and , which is a random matrix of the same kind as itself with zero mean and covariance matrix given above. Thus we can apply the results of references (see also [16, 39, 64, 62]) to to deduce that its average eigenvalue distribution, in the large limit, is uniform in an ellipse, in the complex plane, centered at , with real (imaginary) semi-axis , which leads to the claimed results about the spectrum of .

Substituting the formulas given above for , and we get that the eigenvalues of the Jacobian matrix are uniformly distributed inside the following ellipse

 center =(−CMSμ, 0) horizontal semiaxis =√CMS((ρ+1)σ2+2(1−CM)μ2)√(σ2+(1−CM)μ2) vertical semiaxis =√CMS(1−ρ)σ2√(σ2+(1−CM)μ2) eigenvalue with the largest real part =−CMSμ+√CMS((ρ+1)σ2+2(1−CM)μ2)√(σ2+(1−CM)μ2) . (17)

Thus if is large enough (), the system is always stable (see Figure 3 in the main text).

## 3 Generalisation of the presented results when adding exploitative interactions

### 3.1 Mean Field Equations, Birth Rates and Species Coexistence

In the main text we have presented the mean field equations also for the case of exploitative interactions (described by the matrix ), in addition to the mutualistic/commensalistic ones (given by the matrix ). In particular, each species is characterised by a birth rate that depends on the species concentration and on the species interactions as:

 ω(j,η,M,L)=¯ηj+ϵ1∑k¯ηkMkjθ(¯ηj)+ϵ2∑k¯ηkLkj¯ηj, (18)

where and give the cooperation and exploitation intensity, and is the Heaviside step function, i.e., when and 0 otherwise. From the microscopic dynamics given by rates (18), the evolution of the relative abundance of each species in the limit is described by the mean field equations:

 ddt¯ηs(t) =ϵ1S∑k=1¯ηk(t)Mksθ(¯ηs(t))+ϵ2S∑k=1¯ηk(t)Lks¯ηs(t) −ϵ1¯ηs(t)S∑i,k=1¯ηk(t)Mkiθ(¯ηi(t))−ϵ2¯ηs(t)S∑i,k=1¯ηk(t)Lki¯ηi(t) (19)

for , where and it is conserved by the dynamics. The derivation of the above equations follow the same steps as described in sec. 2.1.

Note that networks and are non overlapping, i.e. the pair of species cannot interact in both mutualistic and exploitative ways. If the matrix is irreducible - and thus the results presented in section hold - and the transition rates given by Eq. ((18)) are positive during the time evolution - a necessary condition in order that the derivation of the mean field is justified - then we find numerically that, even in presence of a large concentrations of exploitative interactions, at stationarity the system still admits an high biodiversity and full coexistence is observed (see Figure 4). Indeed, we have numerically and systematically investigated the number of extinctions in ecological systems with both mutualistic and exploitative species interactions, as a function of different parameters: the average interaction strengths , the connectance , , the network size , etc. . In all these cases we found that, as long as the rates given by Eq. ((18)) remain positive during the evolution, extinctions are not observed (see Figures 5 6).

### 3.2 Analytical justification of the coexistence condition

Here we want to heuristically justify what we have observed numerically. Namely that adding exploitative interactions does not lead to extinctions, as long as the mutualistic network of interactions is present, corresponding to an irreducible matrix, . We argue that, under this hypothesis, when is positive but close to zero the complete mean field equations - where both and are positive - are perturbation of the mean field equation where only mutualistic interaction are present. Since we have proved that a pure mutualistic system has no extinction as long as the matrix is irreducible. Following the notation in the main text, our continuous time Markov process is defined by the rule: a randomly chosen individual is removed and substituted by an individual of the -th species at a rate

 ωj:=ω(j,η,M,L)=¯ηj+ϵ1∑k¯ηkMkjθ(¯ηj):=ωMj+ϵ2∑k¯ηkLkj¯ηj:=ωLj, (20)

where and give the cooperation and exploitation intensity, and is the Heaviside step function, i.e., when and 0 otherwise. As the relative abundance converges to the solution of the system of ordinary differential equation for . Equation for , when is positive but close to zero, can be written in the following form

 ddt¯ηs(t)=ωMs−¯ηs(t)∑iωMi≃δ>0+O(¯ηs) (21)

The first two terms in (21) are the vector fields corresponding to mean field equation for irreducible and no exploitation (i.e. ). We know that such a system has no extinctions and its vector field is typically greater than out of equilibrium when . The last term in (21) contains terms which are linear dependent of which is Thus is positive for close to zero. The requested transition rates never become negative during the time evolution of the mean field equation. This is a necessary condition otherwise the derivation itself of the mean filed equation would be meaningless.

### 3.3 Covariance matrix and Species Interaction Networks

In this section, we consider the normal fluctuations around the deterministic limit of eq. (3.1). This allows us to calculate the matrix describing the correlation between pairs of species population abundances [6]. As highlighted in the main text, this quantity, once opportunely thresholded, is used as an empirical proxy of the species interactions network [1, 34, 2]. Other works, applying maximum entropy approach, use as the quantity to describe species interactions [6, 60]. The aim of this section is to test how well or approximate the true interactions described by in our model.

For sake of simplicity, we assume that the limiting dynamics start at the equilibrium with , . Thus, we define the fluctuation process as

 xiN(t)=√N(¯ηiN(t)−mi) for i=1,…,S. (22)

One can apply standard techniques of convergence of generators to get weak convergence to the thermodynamic limiting evolution. Indeed, the stochastic process converges in distribution to a Gaussian Markov process which solves the stochastic differential equation

 dX=AXdt+ΦdBt (23)

where is a -dimensional Brownian motion, and

 Aij =ϵ1(MTij−δijS∑h,k=1mhMhk−miS∑k=1Mjk) +ϵ2(Ljimi+δijS∑h=1mhLhi−δijS∑h,k=1mhLhkmk −miS∑k=1Ljkmk−miS∑k=1Lkjmk) for i, j=1,…,S

and satisfies the following constraint equation

 (ΦΦT)ij= −2(mimj(1+ϵ1miS∑h,k=1mkMkh+ϵ2S∑h,k=1mkLkhmh))(1−δij) +2(1−mi)(mi+ϵ1S∑k=1mkMki+ϵ2S∑k=1mkLkimi)δij for i, j=1,…,S,

where is the Kronecker delta

From eq. (23), it is then possible to derive the dynamics of the covariance matrix (see [35] for details):

 Vij(t)=⟨Xi(t)Xj(t)⟩−⟨Xi(t)⟩⟨Xj(t)⟩. (24)

Therefore, we have

 dV(t)dt=AV(t)+V(t)AT+ΦΦT, (25)

and at stationarity the covariance matrix, , resolves the following equation

 AV+VAT+ΦΦT=0. (26)

Eq. (26) is a Lyapunov equation, so we could apply standard algorithms to solve it numerically [59].

We have determined from the solution of eq. (26) and determined . If one assume that the population fluctuations around their means are gaussian distributed, then represents the species interaction matrix [49, 6]. Indeed, within a maximum entropy approach, is typically used to infer species interactions based on the available information of the system [60]. In our framework and as shown by eqs.(25) and (26), the relation between the interaction matrix and the matrix or is highly non-linear. Moreover, because of the constraint, , is not invertible, and thus in order to compute we apply a pseudo-inverse scheme, i.e. we invert is the subspace of spanned by the eigenvectors corresponding to non-zero eigenvalues. As shown in Figures 7, even for very simple structure of matrix and , and are not good proxies of the species interactions. The results are shown for the model without empty sites, but there is no qualitatively difference with the model including empty sites. This result highlights the importance to properly infer interaction networks from data.

### 3.4 Stability of the equilibria when ϵ2≠0

In the case of , the entries of the Jacobian read

 Aij =ϵ1(MTij−δijS∑h,k=1mhMhk−miS∑k=1Mjk) +ϵ2(Ljimi+δijS∑k=1mhLhi−δijS∑h,k=1mhLhkmk −miS∑k=1Ljkmk−miS∑k=1Lkjmk) .

The diagonal entries of the Jacobian are

 Aii=−ϵ1S∑h,k=1mhMhk+ϵ2S∑h=1mhLhi−ϵ2S∑h,k=1mhLhkmk . (27)

Since , it is simple to observe that the term proportional to is of order (plus subleading fluctuations). On the other hand, the leading order of the terms proportional to , is of order and therefore always subleading if . A similar argument applies to the off-diagonal elements. In that case, the terms proportional to are of order , while the ones proportional to are of order .

Similarly to what found in the case , we have that the following relations hold: