Mutualism supports biodiversity when the direct competition is weak

A key question of theoretical ecology is which properties of ecosystems favour their stability and help maintaining biodiversity. This question recently reconsidered mutualistic systems, generating intense controversy about the role of mutualistic interactions and their network architecture. Here we show analytically and verify with simulations that reducing the effective interspecific competition and the propagation of perturbations positively influences structural stability against environmental perturbations, enhancing persistence. Noteworthy, mutualism reduces the effective interspecific competition only when the direct interspecific competition is weaker than a critical value. This critical competition is in almost all cases larger in pollinator networks than in random networks with the same connectance. Highly connected mutualistic networks reduce the propagation of environmental perturbations, a mechanism reminiscent of MacArthur’s proposal that ecosystem complexity enhances stability. Our analytic framework rationalizes previous contradictory results, and it gives valuable insight on the complex relationship between mutualism and biodiversity.

where N (P) i denotes the abundance of plant species i, α (P) i is its intrinsic growth rate in the absence of other species, β (P) ij is the direct competition matrix and γ (P) ik is the mutualistic matrix. The equations and parameters for animals are obtained interchanging the superscripts P and A and they will be presented only when necessary. We will sometimes use the superscript X to indicate either group.
We adopt the so-called soft mean field model, which assumes that all equivalent interaction parameters are identically distributed. The competition network is fully connected, while the mutualistic network is sparsely connected, with adjacency matrix a ik , and its interactions are parameterized as γ ik . The parameter γ 0 , which represents the strength of the mutualistic interactions with respect to competitive ones, is chosen in a range such that the dynamical system is globally stable for all the studied mutualistic networks (see subsection Equivalent Lotka-Volterra system and dynamical stability).
Given a realization of the ecological interactions β ij and γ ik , if we randomly extract the intrinsic growth rates α i we obtain with high probability equilibrium abundances that are not positive (i.e. the equilibrium is not feasible), in particular if the degrees of the mutualistic network are broadly distributed and we extract the α i independent of these degrees. Therefore, instead of randomly extracting the intrinsic growth rates, we impose that the equilibrium is feasible and dynamically stable and study how structurally stable it is against random perturbations of the growth rates.
The first step in our numerical experiments consists in choosing positive equilibrium abundances N i > 0 for all species, i.e. we guarantee the feasibility of the equilibrium. We then determine the unperturbed growth rates α (A) i in such a way that these N i are a fixed point of the dynamical system: ( Note that, if the abundances are narrowly distributed, the unperturbed growth rates are negatively correlated with the number of non-zero mutualistic interactions γ ik , which can be interpreted as a trade-off between the number of interactions and their metabolic cost. The effective interaction and growth rates parameters are obtained by differentiating the full dynamical equations Supp.Eq.(1) at the equilibrium point, and they are (5) We see from this equation that for each species there are two regimes of parameters: weak mutualism, in which the equilibrium mutualistic benefit is far from saturation (z i ≪ 1) and strong mutualism, in which the saturation is reached (z i ≫ 1). Thus, it is possible that, for given meta-parameters, species with large degree are in the strong mutualistic regime while species with small degree are in the weak regime. The effective mutualistic strength increases with γ 0 in the weak mutualistic regime and decreases in the strong regime, reaching a maximum in between. The local stability of the equilibrium point can be tested from the linearized system Supp.Eq.(3), which we rewrite as 1 The equilibrium is locally stable if and only if all eigenvalues of the Jacobian matrix J ik = −N i A ik have negative real part. This requires that the elements of the matrices γ LV are small. From Supp.Eq.(4), it is clear that this happens both for small and for large values of γ 0 , i.e. the equilibrium is stable if γ 0 < γ (1) 0 or γ 0 > γ (2) 0 , which define the lower and upper critical mutualistic strengths and the weak and strong mutualistic regimes, respectively. We determine numerically the critical strengths for each given network and given realization of the interaction parameters by diagonalizing the matrix J ik .
We can obtain greater analytic insight on dynamical stability using the effective competition matrix defined below.

The effective competition framework
The effective competition matrix represents the interactions between species in the same group, either P or A, both due to their direct interaction (in this case, competition) and to their interaction with species in the other group (in this case mutualism). It allows to decouple the fixed point equations as where the effective competition matrices C and the effective productivity vectors p are defined as and analogous for plants [4]. The usefulness of the effective competition matrix stems from the fact that, while the full interaction matrix A has both positive and negative components, we expect that most elements of the matrix C (P) and C (A) are positive in the regime in which the system is dynamically stable (for simplicity, we omit here the superscripts since the properties that we describe are common to both matrices). In particular, we assume that the components of C are non-negative and that the matrix is irreducible, i.e. the complete space is the only spaces that is invariant under the action of C, so that we can apply the Perron-Frobenius theorem [5] that states that each matrix C has a dominant eigenvalue λ 1 of order S (the dimension of the space) associated with left and right eigenvectors u 1 and v 1 whose elements are all positive, while all other eigenvectors have at least one negative or complex component.
The main eigenvalue λ 1 can be interpreted as the weighted average of the matrix C with weights given by the main left and right eigenvectors, [4], it is convenient to change abundance units in such a way that the intraspecific competition is equal to one, i.e.
This transformation preserves the fixed point equations, Supp.Eq. (6). From the main eigenvalue of the transformed effective competition matrix we can compute the effective interspecific competition parameter as [4] SinceC ii = 1, we can write ρ eff = i =j C ij v 1 i v 1 j /(S − 1), which shows that ρ eff is proportional to the weighted average of the interspecific competition, justifying its name. Since the trace of the matrix iC ii = S is equal to the sum of the eigenvalues, 1 − ρ eff represents the average value of the minor eigenvalues: α>1 λ α /(S − 1) = 1 − ρ eff . This is the key property that we used for deriving Supp.Eq.(11), based on which we predict structural stability (see below).
Here we adopted an alternative and more general way of determining the parameter ρ eff that does not require rescaling the competition matrix. We have found that this procedure yields slightly more accurate predictions of the structural stability. The formulas reported above are recovered as a special case when C ii = 1. We define ρ eff as the ratio between interspecific and intraspecific competition, Computing the trace, we recover the result that 1 − ρ eff is proportional to the average value of the minor eigenvalues, However, note that with the above definition ρ eff is not anymore proportional to the weighted average of the interspecific competition. The effective competition parameter ρ eff plays a central role in determining the dynamical stability, structural stability and abundance of the model ecosystem, as shown in previous works [4,6,7] and recalled below for completeness.

Types of stability
The stability of the dynamical system Supp.Eq.(1) can be defined in several ways. Dynamical stability refers to perturbations of the dynamical variables N i , and it can be can be local or global. Local stability against small perturbations of the dynamical variables around their equilibrium values is analytically studied considering the linearized dynamical system Supp.Eq.(3). Global stability refers to perturbations of whatever size of the dynamical variables. Goh has shown that, if an interaction matrix A is diagonally positive, meaning that there is a diagonal and positive matrix D such that the symmetric matrix DA + A T D is positive definite, then the Lotka-Volterra system defined by this matrix, 1 [8] 1 . This implies that, for whatever perturbation of the dynamical variable, no species will get extinct provided that the parameters of the system allow feasibility. In this situation, it is possible to analytically investigate a more general type of stability, namely structural stability.
A dynamical system is said to be structurally stable if its qualitative properties do not change when its parameters suffer a perturbation. If the interaction matrix is diagonally stable, then structural stability against perturbations of the growth rates can be defined as the maximum perturbation that maintains a feasible equilibrium in which all species coexist. Since the interaction matrix does not change, its global stability remains guaranteed. We can define in an analogous way the structural stability with respect to changes in the interaction matrix, but in this case we have also to test that the matrix remains diagonally positive.

Effective competition and dynamical stability
The dynamical stability of LV systems is usually addressed in terms of the complete interaction matrix A, which involves both groups of plants and animals. Recently, we have addressed dynamical stability through the effective competition matrix of each group of species separately [6]. For this purpose, we used Goh's theorem mentioned above [8] that states that a LV system is globally stable if its interaction matrix A is diagonally positive (with the definition of A that we adopted; with the conventional definition it must be diagonally negative, see previous section).
In particular, we have shown that diagonal positivity of the effective competition matrix C (P) and the conjugated direct competition matrix β (A) is a necessary condition for diagonal positivity of the complete interaction matrix A, which in turn is sufficient for global stability. If A is symmetric, then the above condition is necessary and sufficient for diagonal positivity of A. If it is not, mild symmetry conditions on the main principal vectors of the intergroup interaction matrices must be also imposed to guarantee the diagonal positivity of A [6]. Furthermore, we have shown in [6] that diagonal positivity of the LV interaction matrix A is also a sufficient condition for global stability of the dynamical equations Supp.Eq.(1) with their saturation terms.
Diagonal positivity is not simple to test numerically, because we must find a suitable matrix D or rule out its existence. However, for the model that we simulated the interaction matrix is almost symmetric, since the direct competition matrices are approximately symmetric and the mutualistic interaction coefficients are the same for plants and animals. Therefore, we conjecture that for our model positivity of C (P) and β (A) is an almost necessary and sufficient for global stability of the linearized system, which implies that no species will get extinct provided that the equilibrium is feasible, as Supp.Eq.(2) guarantees by construction. Since also β (A) is positive definite and almost symmetric by construction, it is enough to test that C (P) is positive definite.
We tested numerically that the threshold value of γ 0 obtained through the condition that C (P) is positive definite, and which is almost sufficient for global stability, approximately coincides with the threshold value obtained from the local stability condition that Re (λ (J ik )) < 0 (see Supp. Fig.1). In all cases that we investigated numerically, local stability also guaranteed global stability.
Because of Supp.Eq.(10), the smaller is ρ eff(P) , the larger is the average of the minor eigenvalues of C (P) and the less likely it is that the minimum eigenvalue of C (P) is negative and the system is dynamically unstable. Therefore, we expect that dynamical stability is inversely related with ρ eff(P) , in particular that the lower mutualistic threshold γ (1) 0 is inversely related with ρ eff(P) , consistent with our numerical results reported in Supp. Fig.2.
This relationship also predicts that the critical mutualistic strengths are influenced by the same network properties that influence ρ eff(P) , i.e. nestedness in the weak mutualistic regime (Supp. Fig.3 left plot) and connectance in the strong mutualistic regime (Supp. Fig.3 right plot).

Effective competition and structural stability
In this section, to simplify the notation we omit the superscripts P and A whenever no ambiguity arises. The structural stability with respect to changes in the productivities is in large extent determined by the effective interspecific competition ρ eff [4,7]. Sufficient condition for all species having positive abundance larger than n c (a threshold abundance below which a species is not viable) is that all productivities fulfill the inequality We call η i the vulnerability factor. Here v 1 is the main eigenvector of the effective competition matrix C, p 1 = i p i v 1 i is the projection of the productivity vector onto v 1 , and S eff is a natural biodiversity scale set by the effective competition,   If S eff /S is small (either large ρ eff or large S) all of the η i must be close to zero (note that the weighted average of η i (p) is i η i (p)(v 1 i ) 2 = 0 for all vectors p), i.e. the productivity vector must be almost parallel to v 1 , so that the feasibility condition Supp.Eq.(11) is very demanding and even small perturbations of the productivities can violate it.
The feasibility condition Supp.Eq.(11) is a sufficient condition. Nevertheless, we have observed with numerical experiments that it becomes an approximately necesary condition for feasibility when the productivity is exposed to random perturbations [7], and we will use it for predicting the critical perturbation ∆ c .
Structural stability and dynamical stability are positively related As a corollary of the above results, in the regime A in which ∆ c is negatively influenced by ρ eff we expect a negative correlation between the dynamical stability with respect to fluctuations of the species abundances, measured through γ (1) 0 , and the structural stability with respect to fluctuations of growth rates, measured through ∆ c . Conversely, in the regimes in which ∆ c is mainly influenced by the propagation of perturbations, both ∆ c and γ (1) 0 increase with the connectance of the mutualistic network and we expect that they are positively related. This implies that the structural stability with respect to variation in the mutualistic interactions and with respect to variations in the intrinsic growth rates go in the same direction in these regimes. This prediction is supported by our numerical results reported in Supp. Fig.4.

Effective competition and equilibrium abundances
In the limit of vanishing S eff /S, the feasibility condition implies that the productivity vector must be directed along v 1 , p i ≈ p 1 v 1 i . Through the equilibrium equation N = (C) −1 p, this in turn implies that the equilibrium abundances are also directed along v 1 : Thus for large systems with large S/S eff whose productivities are strongly constrained, the equilibrium abundances are expected to be inversely related with ρ eff , with the consequence that a decrease of ρ eff positively affects structural stability, dynamical stability and equilibrium abundances.

Connectance and nestedness of the random networks
For each set of meta-parameters we considered 125 mutualistic networks with different combinations of connectance and nestedness. The properties of the simulated networks are depicted in Supp. Fig.5, where one can see that almost all networks have nestedness larger than connectance, a limit that is very difficult to overcome with the Monte Carlo schemes that we implemented.

Regimes not shown in the main text
We show in Supp. Fig.6 how structural stability ∆ c , effective interspecific competition ρ eff and propagation of perturbations η ′ depend on network properties for some regimes of parameters that are not represented in the main text.
Supplementary Figure 6: Structural stability and related quantities versus network architecture (next page). ∆ c quantitfies structural stability, ρ eff quantifies the effective competition and η ′ quantifies the propagation of perturbations. Each panel depicts a different regime versus connectance (left column) and nestedness (right column). Top panels: facultative mutualism, large equilibrium abundance, strong (left) and weak (right) direct competition. Intermediate panels: obligatory mutualism, weak mutualistic strengths, strong (left) and weak (right) direct competition. Bottom panels: obligatory mutualism, large mutualistic strength, strong (left) and weak (right) direct competition. Error bars represent the standard error of the mean over 50 realizations of the parameters and the ecological dynamics. The horizontal lines refers to absence of inter-species interactions, i.e. pure direct competition. Obligatory, ρ=0.05, γ 0 =100, N * =1

Effective competition parameter versus nestedness
We show in Supp. Fig.7 the effective competition parameter ρ eff versus nestedness for different values of γ 0 and ρ. In the weak mutualism regime ρ eff decreases with nestedness, in agreement with the results of the previous section. The critical competition increases with nestedness, and it is given by the point where ρ eff = ρ. Each pair of curves of the same colour refers to a different value of the direct competition parameter ρ. The dashed curve is for pure competition and the solid curve for mutualism. Note that the critical interspecific competition at which ρ eff = ρ depends on the nestedness but not on γ 0 , which only influences the absolute value of the difference ρ eff − ρ. Error bars represent the standard error of the mean over 50 realizations of the parameters and the ecological dynamics.

Supplementary Note 3: Influence of the unperturbed abundance distribution on structural stability
As we explained in the Methods section, we extract the unperturbed equilibrium abundances from distributions with mean value one and coefficient of variation CV (standard deviation divided by mean). We have presented in the main text results obtained with an unperturbed system with a narrow abundance distribution. Unperturbed abundances are drawn from an uniform distribution in [0.85, 1.15] and have small coefficient of variation, i.e. ratio between standard deviation and mean, CV=0.087. This is in line with our previous study [3] and with Ref. [9], where the unperturbed system was constructed such that p 0 ∝ v 1 , from which we expect that the vulnerability of the unperturbed system is η v (p 0 ) ≈ 0, corresponding to a narrow distribution of equilibrium abundances (see below).
The assumption that the distribution of abundances is narrow is justified by the results that we report in this Note, where we find that structural stability tends to vanish when the distribution of the unperturbed abundances approaches a critical value. To put this result in context, we have to consider that here we adopt species-specific units such that the intraspecific competition coefficients β ii , inversely related with the carrying capacities, are equal to one for all species of animals or plants. Starting from the original units of abundance, for instance number of individuals per square kilometre, we obtain a mathematically equivalent system of equations withβ ii ≡ 1 through the change of units Thus, it is possible that the abundances N i are broadly distributed in the original units, but they are narrowly distributed in the rescaled unitsÑ i = N i √ β ii that enforceβ ii = 1. We can restate the result presented in this Note by stating that, if the distribution of abundances is large when measured in species-specific units such thatβ ii ≡ 1, then the vulnerability of the unperturbed system is large and its structural stability is very small. In this Note, we present the analytic computation of the influence of the unperturbed system on structural stability, tested with simulations performed with equilibrium abundances with lognormal distributions with large CV equal to 0.15, 0.30 and 0.95. These numerical experiments were performed in the weak facultative mutualism regime (γ 0 = 0.15 < γ (1) 0 , N = 1), but the analytic predictions are valid for all regimes.

Effective abundances and structural stability
For completeness, we recall here the computation that shows how the unperturbed abundances influence structural stability. Again we omit the superscripts P and A unless there is some ambiguity. The structural stability of the system with respect to changes in the intrinsic growth rates α i depends on how these changes propagate onto the productivities and, through them, onto changes of the vulnerability factors η i , Supp.Eq.(11). This computation is complicated by the fact that a change in α i modifies the equilibrium abundances, and consequently the saturation factors z i and, through them, the effective growth rates and mutualistic interactions, Supp.Eq.(5) and Supp.Eq.(4). Therefore, we would need to explicitly compute the perturbed abundances and z i , which would make the feasibility condition Supp.Eq.(11) useless. Nevertheless, we can get analytic insight by noticing that, when mutualistic interactions are far or close to saturation (z i = h i l γ il N l either small or large), the change in γ LV and α LV due to a change in equilibrium abundances is small and it can be neglected, except for obligatory mutualism (see main text). Thus, we assume that the perturbation of growth rates does not modify γ LV and the effective competition matrices C and that the perturbed α LV is simply given by α LV i (∆) = α LV i (∆ = 0) (1 + ∆r i ), where r i are random numbers. From the α, using Supp.Eq. (7) we compute the productivity vector as p k , where G (P) = γ LV(P) β (A) −1 , and we obtain the perturbed productivity p i (∆) = p 0 i + ∆q i , where p 0 i represents the unperturbed productivity, and q k r k is a random variable. We then obtain p 1 (∆) by projecting the vector p onto the main eigenvector of the unperturbed effective competition matrix v 1 and we compute η i (p(∆)) as where we neglected terms of order 2 or higher in ∆. The feasibility condition depends on the most vulnerable species, which is different in the unperturbed system ∆ = 0 (the species i = v′ that maximizes η i (p 0 )) and for large ∆ (the species i = v that maximizes η i (∆)). If ∆ is large enough, the most vulnerable species coincides with the species that maximizes η i (q) − η i (p 0 ) and it does not change for larger ∆, so that the maximum value of η i (∆) incrases linearly as We determine ∆ c as the maximum value for which the feasibility condition η(∆) ≤ η c ≡ S eff /(S + S eff ) = (1 − ρ eff )/(ρ eff (S − 1) + 1) holds, i.e.
where we have reintroduced the superscript X (either plants or animals) for clarity. We call η v (p 0 ) the vulnerability of the unperturbed system and η′ the propagation of perturbations. Numerically, we compute them as where η(∆) = max i η i (p(∆)) (the overline denotes the average with respect to perturbations) and ∆ 0 , ∆ 1 are in the range where η(∆) is linear (in the main text, we present results with ∆ 0 = η c − 0.25 and ∆ 1 η c + 0.25, since we verified that ∆ is linear in this range). From Eq.(14) we see that η′ and η v are inversely related, since where v is the species that maximizes η v (q) − η v (p 0 ). Supp.Eq.(16) clearly shows that the structural stability vanishes when the unperturbed vulnerability η v approaches the critical value η c = 1 − ρ eff(X) / ρ eff(X) (S (X) − 1)

Numerical results
We present in Supp. For ρ = 0.23 and large CV=0.95 it holds η c < η v (p 0 ) for all networks, thus the predicted ∆ c is negative. In this situation, we observe very small values of ∆ c smaller than 0.01 for all networks, in the limits of the numerical precision. These results support the prediction of our theory that ∆ c vanishes when η v (p 0 ) approaches η c . We then show in Supp. Fig.9 the vulnerability of the unperturbed system, η v (p 0 ), averaged over all of the networks, as a function of the variance of the distribution of the abundances. Clearly, the variance has a strong enhancing influence, i.e. the larger the variance is, the more vulnerable to extinctions is the system and the smaller is its structural stability, justifying our choice to focus our analysis on systems with narrowly distributed unperturbed abundances (we recall one again that the abundances must be narrowly distributed in species-specific units of carrying capacity, not in units of number of individuals per square kilometer).
Finally, we study how η v (p 0 ) depends on the properties of the mutualistic network such as the connectance and the nestedness. One can see from Supp. Fig.10 that η v (p 0 ) increases with the nestedness, as expected, since increasing nestedness the distribution of the mutualistic degree becomes broader, and consequently the unperturbed productivities p 0 are more broadly distributed. For fixed nestedness, η v (p 0 ) increases with the connectance only when the nestedness is large. From the figure, one can see that for CV≈ 0.1 the scale of η v (p 0 ) is one order of magnitude smaller than η ′ , justifying that its effect can be neglected in the results reported in the main text.  Supplementary Figure 10: Vulnerability of the unperturbed system versus network properties.η v (p 0 ) as a function of the connectance (left) and nestedness (right) of the mutualistic networks for CV=0.097 with ρ = 0.05 (top) and ρ = 0.23 (bottom). Each point represents one network with different connectance and nestedness, and for each of them we performed 50 independent simulations. Error bars represent the standard error of the mean over the realizations.

Supplementary Note 4: Other ecological interactions
The computations presented in the paper for mutualistic systems can be easily generalized to other types of ecological interactions. If the two groups of species compete with each other, the sign of the inter-group interaction parameter γ 0 is negative, but γ 2 0 does not change. For competitive interactions, the effective competition parameter ρ eff , which depends on γ 2 0 , would have the same expression as for mutualism. The qualitative interpretation is that shared competitors of the other group tend to reduce the effective competition between species of the same group exactly as shared mutualistic partners do. Also in this case we predict that there is a critical value of the direct competition parameter ρ below which competition with the other group tends to decrease the effective competition, with a beneficial effect on the structural stability, while above the critical value the opposite holds. On the other hand, competitive inter-group interactions reduce the productivity and therefore they are expected to increase the propagation of perturbations, opposite to what mutualistic interactions do, so that we expect that the connectance of the inter-group competition network tends to increase the propagation of perturbations and decrease structural stability. Also in this case, we expect that the net effect of inter-group competition on structural stability depends on the balance between the possible positive effect on ρ eff and the negative effect on η′. For the case of predatory interactions, the interaction is positive for one group and negative for the other group, so that we have to substitute γ 2 0 with −γ 2 0 in the expressions for the effective competition, recovering the known qualitative result that shared preys increase the effective competition between predators. The same happens for shared predators, at least if the functional response of predators is linear. The computation of the effective competition predicts that predatory inter-group interactions decrease the effective competition when the within-group competition ρ is large, with a beneficial effect on structural stability, and increase the effective competition when the within-group competition is low, opposite to what happens with mutualistic and competitive interactions. Interestingly, the predicted value of this critical competition would be the same as for mutualistic and competitive interactions. On the other hand, predatory inter-group interactions increase the productivity of predators but decrease the productivity of preys, thus they are expected to have opposite effects on the propagation of perturbations of predators and preys. We expect that the net effect on structural stability of the connectance of the inter-group predatory network depends on the balance between these two effects, and changes depending on the regime of parameters. In summary, we expect that for weak within-group competition, inter-group interactions of mutualistic and competitive type tend to decrease the effective competition, while intergroup interactions of predatory type tend to increase it, and that the opposite holds for strong direct competition.