Breakdown of Random-Matrix Universality in Persistent Lotka-Volterra Communities

The theory of disordered systems enables one to deduce the behavior of collections of many interacting constituents, whose interactions are assumed to be random, but fixed in time [1]. A related discipline, random matrix theory (RMT), is concerned with the eigenvalue spectra of matrices with entries drawn from a joint probability distribution. Both fields have found numerous applications in physics [2,3] (the study of spin glasses in particular [1]), and in other disciplines such as neural networks [4–9], economics [10,11], and theoretical ecology [12–19].

It is frequently assumed that the distribution of the randomness in RMT or disordered systems is Gaussian, possibly with correlations between different interaction coefficients or matrix entries. Reasons cited for this assumption include analytical convenience, maximum-entropy arguments, and the observation that higher-order moments often do not contribute to the results of calculations [1,20,21].

In random matrix theory, this latter observation is referred to as the principle of universality [22–24]. The principle states that results obtained for the spectra of Gaussian random matrices frequently also apply to matrix ensembles with non-Gaussian distributions. The conditions for universality to apply are usually mild (higher-order moments of the distribution must fall off sufficiently quickly with the matrix size [22,23]), and it is often tacitly assumed that these conditions will hold.

In this Letter, we offer a pertinent counterexample to the universality principle in RMT. We focus on the ecological community resulting from the dynamics of the generalized Lotka-Volterra equations with random interaction coefficients. The stability of this community is governed by the interactions between species that survive in the long run [25,26]. This is a submatrix of the original interactions, which we will refer to as the “reduced interaction matrix.”

Firstly, using dynamic mean-field theory [27], we obtain the statistics of the elements in the reduced interaction matrix. These turn out to be non-Gaussian (even when the original interaction matrix is Gaussian). Secondly, we analytically calculate the leading eigenvalue of this non-Gaussian ensemble of random matrices. We show that this eigenvalue is different from the one that we would obtain from a Gaussian ensemble with the same first and second moments as in the reduced interaction matrix. This demonstrates that the principle of universality fails, and it indicates that the Gaussian assumption should not be made lightly.

Our findings have relevance to the random matrix approach to ecosystem stability, introduced by Robert May [12,13]. This approach assumes a random interaction structure between species in the community. One line of criticism of May’s model is the observation that such interactions do not necessarily describe a feasible equilibrium (that is, an equilibrium for which all species’ abundances are positive) [25,28–31]. The community of surviving species in the generalized Lotka-Volterra model on the other hand is feasible by construction, and we derive the statistics of the emergent random matrix ensemble that describes this community [26,32–34]. From this ensemble, we then recover the stability criteria that have previously been derived from the dynamic Lotka-Volterra model [16,35]. We thus show that one can construct a random matrix ensemble (in the sense of May) that correctly reflects the stability of a feasible community of coexistent species. This ensemble is non-Gaussian and quite intricate. In May’s words, our work contributes to “elucidating the devious strategies of nature which make for stability in enduring natural systems” [36].

We start from the generalized Lotka-Volterra equations (GLVEs) [16,35,37,38]

where the 𝑥𝑖 ≥0 describe the abundances of species 𝑖 =1,…,𝑁. The interaction matrix elements in Eq. (1) 𝑎𝑖⁢𝑗 are quenched random variables. We refer to these as the “original interaction matrix” elements. We assume that the mean of each matrix element is ‾𝑎𝑖⁢𝑗=𝜇/𝑁 (we use an overbar to denote averages over the ensemble of interaction matrices), and that they have variance Var⁡(𝑎𝑖⁢𝑗) =𝜎2/𝑁. We also allow for correlations between diagonally opposed matrix elements, Corr⁡(𝑎𝑖⁢𝑗,𝑎𝑗⁢𝑖) =Γ, ( −1 ≤Γ ≤1) where Corr⁡(𝑎,𝑏)=(‾𝑎⁢𝑏−¯𝑎⁢¯𝑏)/√Var⁡(𝑎)⁢Var⁡(𝑏).

The scaling with 𝑁 of the moments of 𝑎𝑖⁢𝑗 follows the standard conventions in disordered systems [1] and guarantees a well-defined thermodynamic limit 𝑁 →∞. All our results are independent of the higher moments of 𝑎𝑖⁢𝑗 as long as these moments decay sufficiently quickly with 𝑁. Further details can be found in Sec. S1 of the Supplemental Material [39].

Previous analyses of this system [16,35] in the thermodynamic limit have shown that there is a range of parameter combinations 𝜇,𝜎2, and Γ for which the dynamics reaches a unique stable fixed point, independently of the starting conditions. This is the case in the region to the left and below the instability lines in the phase diagram in Fig. 1.

FIG. 1.

When a fixed-point solution is reached, not all species survive, i.e., there are some species for which 𝑥⋆𝑖 >0 and others with 𝑥⋆𝑖 =0 (we use an asterisk to denote the fixed point). Using dynamic mean-field theory (DMFT) [15,40–42], one can deduce these statistics of the species’ abundances at the fixed point.

From the DMFT analysis, one can also find the combinations of system parameters at which the system is no longer able to support a unique stable fixed point. There are two types of transitions: (1) the average species abundance can diverge [i.e., 𝑀 →∞], or (2) the fixed-point solution can become linearly unstable to perturbations. Closed-form expressions for the critical sets of parameters ( 𝜎, Γ, and 𝜇) at which each of these transitions occurs were derived in Refs. [16,35]. A selection of phase lines for different values of the correlation parameter Γ are shown in Fig. 1.

We now examine an alternative approach to analyzing the stability of the GLVEs in Eq. (1). Namely, we consider the reduced interaction matrix (the interaction matrix between the species in the surviving subcommunity). More precisely, this is defined by

where 𝑖,𝑗 ∈𝒮 (with 𝒮 the set of surviving species), and where the shift in the diagonal elements reflects the −𝑥𝑖 term inside the brackets of Eq. (1). It can be shown that a fixed point of the GLVEs is stable if and only if all of the eigenvalues of the reduced interaction matrix have negative real parts [17,25,26] (see also Sec. S2 in the Supplemental Material).

We note that the statistics of the reduced interaction matrix elements are determined by the extinction dynamics in the GLVE system, and are consequently vastly different from those of the original interaction matrix [33,43]. For instance, they are non-Gaussian (even when the 𝑎𝑖⁢𝑗 are Gaussian), and there are correlations between elements sharing only one index (see the Supplemental Material [39], Sec. S6). This makes the calculation of the eigenvalue spectrum of the reduced interaction matrix a nontrivial task.

As is illustrated in Fig. 2, the spectrum of the reduced interaction matrix consists of a bulk set of eigenvalues and a single outlier. Writing 𝑧𝑖⁢𝑗 =𝑎𝑖⁢𝑗 −𝜇⁢𝑁−1 (where once again 𝑖,𝑗 ∈𝒮), both the outlier eigenvalue 𝜆outlier and the bulk spectral density 𝜌bulk⁡(𝜆) can be obtained from the resolvent matrix 𝐺__=[𝜔⁢𝟙__−𝑧__]−1. The bulk density is calculated from the trace of 𝐺__ via well-known relations [44]. The outlier eigenvalue in turn fulfills [45–47]

where 𝒢⁡(𝜔)≡(𝑁⁢𝜙)−1⁢‾∑𝑖,𝑗∈𝒮𝐺𝑖⁢𝑗⁡(𝜔), and where 𝜙 is the fraction of surviving species at the fixed point.

FIG. 2.

We first briefly discuss the bulk spectrum, for which the results do not run counter to the universality principle. We use a series expansion for a Hermitized version of the resolvent of the reduced interaction matrix. This standard approach accounts for the nonanalytic nature of the resolvent in the bulk region [48–50].

We find that the resulting series for the trace of the resolvent matrix is identical to that of a Gaussian random matrix in the limit 𝑁 →∞. That is, we show that the higher-order statistics of the reduced interaction matrix do not contribute to this series and, therefore, that the universality principle holds for the bulk region. The only statistics of the reduced interaction matrix that contribute are (𝜎′)2 ≡𝑁𝒮⁢Var⁡(𝑎′𝑖⁢𝑗) =𝜙⁢𝜎2 and Γ′ ≡Corr⁡(𝑎′𝑖⁢𝑗,𝑎′𝑗⁢𝑖) =Γ where 𝑁𝒮 is the number of surviving species (we calculate these statistics in Sec. S6 of the Supplemental Material). One obtains the familiar elliptic law [51]

where 𝜆 =𝑥 +𝑖⁢𝑦. We can show (see the Supplemental Material [39], Sec. S5C) that the bulk of the eigenvalue spectrum crossing the imaginary axis corresponds to the linear instability of the GLVEs, represented by the dashed horizontal lines in Fig. 1. This is verified in Fig. 3(a).

FIG. 3.

We now move on to the outlier eigenvalue, which is a far less trivial matter. We first discuss two candidate expressions for the outlier eigenvalue based upon calculations for Gaussian random matrix ensembles. We show that neither of these expressions is accurate, and that the universality principle fails to predict the outlier eigenvalue. We subsequently derive an accurate expression for the outlier, which we show correctly predicts stability.

Noting previous work [14,24,45,52], one might perhaps expect that 𝜇′ =𝑁𝒮⁢‾𝑎′𝑖⁢𝑗 ( 𝑖 ≠𝑗), together with (𝜎′)2 and Γ′ would be sufficient to predict the outlier eigenvalue of the reduced interaction matrix. Using an established formula for the outlier eigenvalues of Gaussian random matrices [45,52], one then obtains 𝜆0 =−1 +𝜇′ +Γ′⁢𝜎′2/𝜇′.

If we also include the effects of correlations between elements sharing only one index 𝛾′ =𝑁2⁢Corr⁡(𝑎′𝑖⁢𝑗,𝑎′𝑘⁢𝑖) (where 𝑘 ≠𝑖), we arrive at (using results from our previous work [47])

The approach leading to Eq. (5) takes into account all possible correlations for a Gaussian random matrix with statistical symmetry between different species. We note that correlations between elements in the same row or column also exist in the reduced interaction matrix (see the Supplemental Material [39], Sec. S6A), but these do not affect the location of the outlier [47].

If the universality principle were to apply to the reduced interaction matrix, then the Gaussian prediction 𝜆1 and the true outlier eigenvalue would coincide, whether or not the elements of the reduced interaction matrix were also Gaussian distributed. As can be seen in Fig. 4, 𝜆1 is a better approximation than 𝜆0, but neither expression correctly predicts the outlier.

FIG. 4.

We now take into account the full statistics of the matrix elements 𝑎′𝑖⁢𝑗, as we did when calculating the bulk eigenvalue spectrum, and deduce the correct expression for the outlier eigenvalue. In the region of the complex plane outside the bulk (where the outlier resides), the resolvent can be expanded as a series in 1/𝜔 [Eq. (S36) in the Supplemental Material]. We evaluate each term in this series in terms of the statistics of species abundances, which are available to us via DMFT. This is accomplished via a generating-functional approach (Supplemental Material [39], Sec. S4).

Using diagrammatic techniques to recognize the self-similarity of the resulting series [48–50,53], we arrive at a compact formula for the resolvent [Supplemental Material [39], Eq. (S69)]. Using Eq. (3), we then obtain an implicit set of equations for the outlier eigenvalue in terms of the statistics of the surviving species abundances [see Eqs. (S71)–(S73) in the Supplemental Material]. We emphasize that in finding our final expression for the outlier, no approximations have been made other than assuming the thermodynamic limit. The simulation data in Figs. 3 and 4 verify that the expression in Eqs. (S71)–(S73) accurately predicts the outlier eigenvalue.

We also demonstrate analytically (see the Supplemental Material [39], Sec. S4D) that this prediction for the outlier eigenvalue correctly predicts instability of the fixed point of the GLVE system. That is, 𝜆outlier crosses the imaginary axis precisely at locations in parameter space where the 𝑀 →∞ transition occurs in the GLVEs. This is also verified in Figs. 3 and 4.

We thus conclude that stability cannot be predicted from the reduced interaction matrix using Gaussian random matrix results, even if all correlations are accounted for. This indicates that the extinction dynamics leads to some more intricate structure to the interactions in the surviving community.

Advancing ideas in Refs. [26,33,54], we show in the Supplemental Material [39] (Sec. S10) how one can generate the ensemble of reduced interaction matrices “from scratch” (i.e., without running the Lotka-Volterra dynamics and eliminating extinct species). This is achieved by first drawing a set of mock abundances from the known distribution of GLVE fixed-point abundances [16,33]. Subsequently, one then draws interaction matrices from a carefully constructed distribution, which is dependent on the mock abundances. We verify in the Supplemental Material [39] that this bottom-up construction leads to non-Gaussian matrices with the same statistical properties and leading eigenvalue as the ensemble of true reduced interaction matrices.

Having constructed the reduced interaction matrix ensemble in this way, we can thus see more clearly why universality fails to capture stability. The ensemble is manifestly non-Gaussian with complex interdependencies between matrix elements. By making a simple Gaussian assumption and ignoring the higher-order moments, one does not correctly take into account this intricate underlying structure.

Finally, we perform some additional tests of our results to demonstrate their robustness. For example, realistic ecological communities might be composed of only a relatively small number of species. We have verified that our expression for the outlier in Eqs. (S71)–(S73) of the Supplemental Material [39] is also a better predictor of stability than the more naive theories when 𝑁 =50, leading to communities of about 25 surviving species (Fig. S4 in the Supplemental Material). It has also been pointed out that heterogeneity of carrying capacities across species can significantly affect ecological equilibria [33,55]. We show in Sec. S9 of the Supplemental Material [39] that our conclusions continue to hold in such situations.

To conclude, we have deduced the stability of the generalized Lotka-Volterra system by calculating the eigenvalue spectrum of the interaction matrix of the surviving species. We have shown that results that are derived for Gaussian random matrices, which are often assumed also to apply to non-Gaussian ensembles, fail in this case. Instead, higher-order statistics of the reduced interaction matrix must be taken into account. We have therefore found a noncontrived class of random matrices for which the universality principle of RMT is not applicable. This demonstrates that there are limitations to results in RMT that are derived making an assumption of Gaussian interactions. Universality should therefore not be invoked without careful consideration.

Our results also have immediate relevance for the field of theoretical ecology. In the widely used approach pioneered by Robert May [12,13], one supposes that the Jacobian governing small deviations of species abundances about a fixed point can be represented by a random matrix. May does not say what the dynamics are that lead to this Jacobian. One particular objection to this approach is hence that the statistics of May’s random matrices do not necessarily correspond to “feasible” equilibria [24,25,28,30].

The fixed point of the GLVEs is feasible by construction. Therefore, our work shows that the stability of a feasible equilibrium in a complex ecosystem can be found by studying the eigenvalues of a random interaction matrix. Feasibility is reflected in the higher-order statistics of the interactions between species. Crucially, we find that these intricate statistics cannot be ignored if one is to correctly predict stability.