Diversity increases the stability of ecosystems

In 1972, Robert May showed that diversity is detrimental to an ecosystem since, as the number of species increases, the ecosystem is less stable. This is the so-called diversity-stability paradox, which has been derived by considering a mathematical model with linear interactions between the species. Despite being in contradiction with empirical evidence, the diversity-stability paradox has survived the test of time for over 40+ years. In this paper we first show that this paradox is a conclusion driven solely by the linearity of the model employed in its derivation which allows for the neglection of the fixed point solution in the stability analysis. The linear model leads to an ill-posed solution and along with it, its paradoxical stability predictions. We then consider a model ecosystem with nonlinear interactions between species, which leads to a stable ecosystem when the number of species is increased. The saturating non linear term in the species interaction is analogous to a Hill function appearing in systems like gene regulation, neurons, diffusion of information and ecosystems The exact fixed point solution of this model is based on k-core percolation and shows that the paradox disappears. This theoretical result, which is exact and non-perturbative, shows that diversity is beneficial to the ecosystem in agreement with analyzed experimental evidence.


Introduction
The relationship between species diversity and stability in ecosystems has been extensively studied in the literature [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The pioneering study led by Sir Robert May [1] predicted that the more diverse an ecosystem is, the more unstable it is. May's claim resonated powerfully among ecologists as it contradicted the biological principle that great variety of species (and genes) promote ecosystem stability in the face of external stress, and this foundation turned May's claim into a paradox, referred to as the diversity-stability paradox. For almost 40 years this paradox was not able to be refuted, despite evidence showing that ecosystems in nature that have a high degree of diversity tend to be more stable [14]. It was only until recently that the concept of high diversity linked to stability started to emerge; supporting the idea that increasing species diversity is positively correlated with increasing stability at the ecosystemlevel [13,15,17,18] and negatively correlated with species-level stability due to declining population sizes of individual species [14]. However, so far there has not been a theoretical proof that demonstrates mathematically the reason why this occurs. In this article, we show that, using a nonlinear interactions model, the system becomes more stable when there is more species diversity, a statement that differs from the results of the linear model derived from May. In section II we first derive the diversity-stability paradox explicitly by solving the linear model studied by [9][10][11][12][13], which follows the same reasoning as [1]. We show that the solution of the linear model diverges for certain values of the interaction species and thus, it's ill posed. In Section III we propose a nonlinear solution based on a model proposed by [19][20][21][22] and developed in [16] by analyzing the solution of this model, we illustrate that when the interaction strength between mutualistic species is positive and strong, more species in the ecosystem survive. Both the solutions of the linear and nonlinear model are applied to real world ecosystems with positive mutualistic interaction terms between species, so to give a practical example of the two different conditions for stability. In Section IV we present a discussion of the results. We will see that the experimental evidence will support the use of the nonlinear model as a more accurate description of the ecosystems' stability.

Solution of the linear model for arbitrary adjacency matrix
We will first show the solution of the linear model diverges for given values of the interaction species.
In general the evolution of species abundances x i (t) in a ecosystem can be described by dynamical equations of the form: The linear model for ecological networks is described by an adjacency matrix A ij (with A ij = 1 if i and j are connected by a network link, and A ij = 0 otherwise), and linear interactions between species given by: g i ðx 1 ; ::::; x n Þ ¼ The dynamics of species densities x i is then described by the following dynamical system of equations: where b i > 0 is the growth rate of species i, s is the self limitation term representing the selfinteraction of species, that we set equal for all species, γ ij is the strength of the interaction between species i and j, and N is the total number of interacting species. The fixed point equations dx i ðtÞ dt ¼ 0 admit a trivial solution x � = 0, which represent the extinction of all species, and a non-trivial solution x � 6 ¼ 0 which, in implicit form, is given by the following linear system: where Γ ij = γ ij A ij . If s= 2SpecðĜÞ, then the matrix (sI − Γ) is invertible and we can write the solution x � as: where μ a are the eigenvalues ofĜ and L a , R a the corresponding left and right eigenvectors.
Eq (5) shows that the fixed point solution of the linear model has a singularity whenever there is an eigenvalue of Γμ a = s, for some a 2 {1, . . ., N}. In particular, we can think of a situation where the ecosystem is going through a period where the strength of the interactions γ ij is increasing. In this case the linear model becomes ill-defined when the largest eigenvalue of the matrixĜ equals s, μ max = s, because the species densities diverge, x � i ! 1 (see Fig 1). Through the analysis of the solution of the fixed point we can also find that the condition of stability; Eq (5) is feasible if and only if the densities x � i are all positive. This is certainly true when μ max is smaller than s. But when it is close to it (i.e. 0 < s − μ max ⪡ 1), the sum on the r.h. s of Eq (5) is dominated by the term containing μ max , thus giving SinceĜ is an irreducible matrix with non-negative entries, then, by the Perron-Fronebius theorem, the right and left eigenvectors R (max) and L (max) have all positive components, so the vector x � does have strictly positive components, too. On the contrary, when μ max becomes larger than s, all densities turn negative, x � i < 0, and the solution x � becomes unfeasible hence the stability condition is given by s − μ > 0. Next we will show that this condition for stability can also be found via a local stability analysis of the dynamical system.

Local stability of the fixed point solution
The criterion for ecosystem stability is given by the sign of the largest eigenvalue of the stability matrix calculated at the fixed point x � for the dynamical system of Eq (3), which is expressed by the Jacobian J ij ðx � Þ ¼ @ _ x i @x j j x � . Negative eigenvalues indicate that the system is stable. That is, if one of the eigenvalues of the Jacobian is positive, the average may be positive, or zero, and in that case, the system is not stable.
The Jacobian is expressed by: The eigenvalues of J ij (x � ) are not simply related to those of Γ ij , due to the multiplicative term x � i in Eq (7). However, whenĜ is symmetric, we can use the following strategy to infer the crucial properties of the eigenvalues ofĴ . First, we define the matrix Next, we observe thatĴ is similar to the symmetric matrix J ¼X 1=2MX 1=2 , soĴ andĴ have the same eigenvalues. The crucial point is thatĴ andM are congruent matrices, and therefore, by Sylvester's law of inertia, they have the same number of positive, negative, and zero eigenvalues. Therefore, ifM has all negative eigenvalues,Ĵ also has all negative eigenvalues, hence, by similarity, also J has all negative eigenvalues. On the other hand, when μ max = s, we knowM has a zero eigenvalue, but then also the JacobianĴ must have a zero eigenvalue, which means that the solution x � is not stable anymore, as we anticipated in the previous section. It is interesting to observe for which cases of the ratio s g the system is not stable. Fig 2 shows a plot of the sign of the maximum eigenvalue as a function of the interaction term taken for the real networks numbered 1 to 9 in Table 1.
To explicitly re-derive the validity of the solution of the fixed point equation we analyzed the spectrum of the matrixM instead of the spectrum of the Jacobian directly in order to avoid incurring into computation problems at the singularity of x � . The condition to be verified is the sign of the real part of the maximum eigenvalue ofĴ . If this sign is negative the system is stable. If zero or positive, the system is unstable. In Fig 2 we fix the cooperation value to the average of all interactions of the system γ ij = γ and plot the sign of the maximum eigenvalue of matrixM (l M max ) for 9 different ecosystems shown in Table 1 as a function of s g : The figure shows that l M max can be positive, negative or zero when γ is varied, which in this case is a scalar, hence, according to what we previously noted,Ĵ will also have a zero maximum eigenvalue which occurs at the critical condition when μ max = s and changes sign for varying γ. According to Eq 7, the maximum eigenvalue ofĴ (l J max ) will have the same sign of l M max , hence the condition for stability is given by l J max < s g . Note that l J max 6 ¼ m max where μ max is the maximum eigenvalue of Γ but the two conditions of stability are equivalent since the point in which the sign of the eigenvalues ofĴ andĜ change are equivalent, as shown in   It is worth mentioning that, as shown in Fig 2, that l J max changes sign as a function of the interaction term γ. The point at which the real part of the eigenvalue l J max becomes negative indicates the position of the transition from the stable phase (x > 0), to the unstable phase (x < 0) (as shown in Fig 1). In particular, we note that the condition for stability inferred from the analysis of the jacobian coincides exactly with the critical point obtained directly from the analysis of the required positivity of the density of the fixed point x � brought out in the previous section.

Condition of stability through Wigner's law
May's approach, which is usually adopted also from more recent studies of linear model [9,10] is the application of the analog of Wigner's semicircle law for asymmetric matrices, the circular law [1]. This law states that for self regulating systems where the diagonal elements are such that J ii = s < 0, and the off-diagonal elements J ij are independent and identically distributed random variables, with zero mean and variance σ 2 , the eigenvalues ofĴ lie in a disk of radius r J ¼ s ffi ffi ffi ffi N p for N ! 1 centered in −s. Ecological systems are usually only sparsely connected. Hence, both May [1] and [9,10] introduce the connectance C in their calculations; C measures the probability that species interact, consequently the probability of no interaction is given by 1 − C. In this case the circular law states that r J ¼ s ffi ffi ffi ffi ffi ffi ffi NC p . Applying the condition for stability that Reðl J max < 0Þ gives: For N ! 1, the radius of the disk and hence the maximum eigenvalue is l J max � ffi ffi ffi ffi ffi ffi ffi NC p and σ can be seen as the average interaction strength γ hence l J max < s=g is the condition for the local stability of the feasible equilibria x � 6 ¼ 0. In other words, when l J max > s=g, the nontrivial fixed point x � 6 ¼ 0 is unstable (and unfeasible since the average species density hx � i would be negative: hx � i < 0). This stability of the average hx � i < 0 is shown in Fig 1. These results lead to the paradox, since when N increases the system becomes more unstable.
Thus, in the stable feasible region of the linear model Eq (3), the condition: l J max < s=g ðcondition of stability in the linear modelÞ; ð9Þ holds true.
We have then shown that all three methods of study of stability for the linear model produce the same stability condition: the so called diversity-stability paradox.
We test the stability condition by analyzing 9 real mutualistic networks (with positive interactions) compiled from available online resources and detailed in Table 1.
We are able to test the stability phase diagram since for these networks the parameters of the model are provided, in particular the strength of interactions which is the parameter that control the stability of the linear ecosystem via Eq (9). Fig 3 shows the phase diagram for stability predicted by Eq (9) in terms of the values of ðs=g; l A max Þ for each of the 9 analyzed mutualistic ecosystems.
We find that the stability condition of the linear model, Eq (9), is not satisfied by these real ecosystems. That is, all real mutualistic networks are located in the unstable region l A max > s=g, as seen in Fig 3, and thus, according to the linear model (3), all systems should collapse. Below we will explain in more detail the nonlinear dynamical model, which predicts opposite results for the condition of stability with respect to the linear model and explains the existence of the 9 real mutualistic networks, suggesting the nonlinear model as a more adequate study of ecological systems.

Stability for the nonlinear functional response
Most of the studies on stability for ecosystems have been done using the linear model explained in Section, mainly because one can find an analytical solution to the fixed point  Table 1. All the networks lie in the unstable region μ max > s/γ, and hence they do not satisfy the condition of the linear model μ max < s/γ, which is necessary to have a feasible (i.e. hx � i > 0) and stable solution. Hence the linear model of Eq (3) predicts that these 9 existing ecosystems should indeed collapse (i.e. hx � i = 0 for all of them) in contradiction to the fact that they are real feasible ecosystems present in nature. https://doi.org/10.1371/journal.pone.0228692.g003

PLOS ONE
Diversity increases stability equation and the stability condition is directly related to the eigenvalues of the adjacency or jacobian matrix. On the other hand, our group has previously analyzed mutualistic ecosystems using network theory and found an exact solution of the nonlinear Type II functional responses ( [16]), where the ratio of species consumed as a function of the species' population is expressed by a term that saturates featuring a more realistic situation when, even if the size of the species is increased, the number of species depleted remains constant at saturation. This behavior is common for the description of gene regulation, neurons, diffusion of information and ecosystems as presented in their article. The dynamics of species densities, x i (t), interacting via the network A ij , is described by the following set of nonlinear differential equations [19][20][21][22]32]: where d > 0 is the death rate of the species, s > 0 is a self limitation parameter directed to the limitation of species' growth when x i exceeds a certain value, α is the half saturation constant, and γ > 0 is the average interaction strength of the nonlinear interaction term. All these dynamical parameters (γ, d, s, α) are discussed in depth in previous works [19][20][21][22]32].
To study the stability of the solution, one has to first find the nontrivial fixed point x � 6 ¼ 0, which has been obtained in [16]. Using a simple logic approximation on the saturating term the solution of the dynamical equations for constant interaction term γ is given by: where y � i represents the reduced density, the Heaviside function Θ(y) = 1 if y > 0 and zero otherwise, and K γ is the threshold on the mutualistic benefit, where the subscript emphasizes its main dependence on γ, the interaction term. The interaction term can be rewritten as a Hill function of degree n = 1. For degree n ! 1, the Hill function can be replaced by a Heaviside function. Even if one approximates the interaction term to the Heaviside function, it is possible to compare the solution given by exact numerical simulations to the solution of the approximated method of Eq (11) and show consistency within a 12.5% error.
We have said that K γ in Eq (11) represents the threshold of the Θ-function; species i interact with species j just if the reduced density y � j is above K γ . When γ is small, which means that the interaction is weak, K γ is large and a smaller number of species j survive (e.g a small number of species have densities such that x � j > K g for weak γ). For γ = γ c , K g c is too large so that no mutualistic benefit between species can be exchanged; at this point the system collapses to the fixed point solution x � = 0. At this stage the K g c is given by In contrast, when the interaction is strong the threshold for the mutualistic benefit K γ is low and a large number of interacting species j survive. After a trimming process of the species in the network, the solution of Eq (11) is shown to be: where K γ -core of a network is the subset of nodes in a network that have degree of at least K γ (integer number), therefore its the most connected subgroup of the graph. According to Eq (13), the tipping point of a mutualistic ecosystem, whose motion is described by Eq (10), is given by the extinction of the species that belong to the K g c -core of a network, expressed in Eq (12) which allows a relation between the dynamical properties of the mutualistic network and a topological invariant of the system, the k-core. The density y � i is positive and therefore N i ðK g Þ is also positive. This condition can only be satisfied if K g < k max core , where k max core is the maximum k-core of the network. Consequently, if this condition is not satisfied the system collapses to the fixed point solution x � i ¼ 0, and the tipping point, described in Eq (12), occurs when which relates the dynamical parameters of (10) to the structural properties of the networks. As a consequence, the stability condition is given by: k max core > K g / s=g ðcondition of stability in the nonlinear modelÞ; ð15Þ According to the solution of the nonlinear model, the larger the k-core number k max core (i.e., the more k-shells in the network) the larger the resilience of the system against external global shocks that reduce the interaction strength γ. With the solution Eqs (11) and (13) of the fixed point equations, we can now study the local stability of the type II dynamic equations by analyzing the Jacobian of the stability matrix To guarantee the stability of the fixed point one has to verify that the real part of the eigenvalues of (7) are all negative. The eigenvalues l J i of J are: and we can easily see that all eigenvalues are negative. The maximum eigenvalue, which is obtained when the nodes (i.e species) of the network have fewest number of edges with the K γcore, is given by l J max ¼ À gðK g þ 1Þ À 1 , is evidently always negative; therefore when the solution is feasible, according to Eq (15), it is always stable. Fig 4 plots the phase diagram of ecosystem stability in the space ðK g ; k max core Þ predicted by the nonlinear model and features the 'tipping line', which separates the feasible stable phase of the nonlinear model (15) from the collapsed phase k max core < K g . Here we plot the values of ðK g ; k max core Þ obtained from the 9 real mutualistic networks illustrated in Table 1, for which these two parameters have been measured in the literature [9,21]. All mutualistic networks lie in the stable feasible region situated below the tipping line, in agreement with the dynamical theory of the nonlinear model, and in contrast to the result of the linear model. We have then shown that by taking into account the actual fixed point solution in the stability analysis, along with the saturation effect of the interaction term, this resolves the diversity-stability paradox.

Discussion
We have presented two different approaches for the study of the stability condition in ecosystems and have seen that in the case of a fully connected network, the linear approach [1, 9, 10, (11) of the nonlinear dynamical system for the same 9 empirical mutualistic networks used in Fig 3. The tipping line is plotted according to the nonlinear model which predicts that the tipping point is given for k max core ¼ K g . All the networks lie in the region k max core > K g and thus they are stable, as they should. https://doi.org/10.1371/journal.pone.0228692.g004 33] leads to counterintuitive predictions which are in contrast with the exact solution of the full nonlinear model (Eq (10)).
The linear model of Eq (3) contains the diversity-stability paradox [1], for which, a more diverse ecosystem is closer to the point of turning unstable [1]. According to the linear model, increasing the diversity has, in general, a destabilizing effect on the mutualistic ecosystem, since it requires the interacting term γ to be smaller. We can then state that the linear stability analysis of the ecosystem features two main features: first, it cannot detect the tipping point of the system collapse Eq (11) and second, the stability analysis [1,9,10,33] is in contrast with the exact solution of the full nonlinear model (Eq (10)) described by a saturating function. The evidence for such controversy is provided in Fig 3, which plots the maximum eigenvalues of 9 real networks which are known to be stable as lying in the unstable regime. According to the linear model then, all real mutualistic networks of Table 1 would be collapsed, and, for the 9 networks studied, this is not the case.
On the other hand, the effect of considering a saturation Hill function in the interaction term, (Eq 10) leads to the opposite stability condition: in a fully connected mutualistic network, Eq (15) read N > s/γ. In this case, by increasing the species diversity N, the condition N > s/γ is easier to satisfy. Similarly, by increasing the mutualistic interaction γ, the stability condition N > s/γ is also easier to satisfy. In conclusion, the nonlinear model predicts both diversity and mutualism to have a stabilizing effect on the whole ecosystem and correctly predicts that the analyzed ecosystems should be feasible as shown in Fig 4. Thus, the effect of the nonlinear model is then crucial to predict the stability and feasibility of the ecosystem.
Furthermore, the linear approximation predicts that more diverse systems (i.e. systems with larger λ max due to either larger connectivity k or larger number of species), are closer to collapse. Analytically, the origin of May's paradox can be traced back to a mathematical singularity in the linear model at the tipping point (Fig 1): the density x � i diverges at s/γ c = μ max , and then collapses instantaneously to the state x � i ¼ 0. This singularity is absent in the nonlinear model due to the saturation effect of the nonlinear interaction term, thus resolving the paradox of the linear model for two main reasons. First, because Eq (15) predicts that the larger the mutualistic strength γ, the more stable the system is. Second, increasing diversity by the number of connections or the number of species, increases the maximum k-core, (or at least leaves it unchanged), thus increasing the robustness of the system. Therefore, stronger mutualistic interactions and augmented diversity stabilize the system, as confirmed by real ecosystems in Fig 3. Thus, all these reasonings indicate the importance of considering the full set of nonlinear interactions when reaching conclusions on the stability of ecosystems. For instance, recent studies [10] have used the linear model to analyze the stability of the microbiota, and have concluded that mutualism in bacteria species is detrimental to the ecosystem. Such a conclusion would be reversed if one were to use the nonlinear model to analyze the data.