A symmetry-based approach to species-rich ecological communities

Abstract

Disordered systems theory provides powerful tools to analyze the generic behaviors of high-dimensional systems, such as species-rich ecological communities or neural networks. By assuming randomness in their interactions, universality ensures that many microscopic details are irrelevant to system-wide dynamics; but the choice of a random ensemble still limits the generality of results. We show here, in the context of ecological dynamics, that these analytical tools do not require a specific choice of ensemble, and that solutions can be found based only on a fundamental rotational symmetry in the interactions, encoding the idea that traits can be recombined into new species without altering global features. Dynamical outcomes then depend on the spectrum of the interaction matrix as a free parameter, allowing us to bridge between results found in different models of interactions, and extend beyond them to previously unidentified behaviors. The distinctive feature of ecological models is the possibility of species extinctions, which leads to an increased universality of dynamics as the fraction of extinct species increases. We expect that these findings can inform new developments in theoretical ecology as well as for other families of complex systems.

I Introduction

Ecological communities are archetypical complex systems, with a vast number of degrees of freedom and a variety of processes and scales. Central to ecology, the question of how interactions between species impact their abundances and thus lead to a variety of patterns and dynamics is made difficult by the amount of parameters required to model such a system.

Building upon tools from spin-glass theory, this problem has recently seen great theoretical progress, with a focus on reducing the number of parameters to a tractable one, that allows for analytical predictions and for a comprehensive study of the impact of each parameter. The main idea, which was pioneered by May in ecology, is to substitute as much as possible of the complexity of the system by randomness[1, 2]. This approach, rooted in the success of so-called ’disordered models’ in theoretical physics, relies on drawing species interactions from some probabilistic model with pre-defined statistics. One then finds that, despite each sample being different, collective properties (i.e. those defined at the level of the whole community) are generic, that is, dependent only on the few statistical parameters that characterize the random model. Hence one can predict important features such as Species Abundance Distributions (SAD) [3], ecosystem stability [4] and the nature of dynamical regimes, without the need to resolve interactions at the species level [5].

Whenever there is little ecological information, disordered models may thus be seen as null models that allow to test what can generically be expected from a community [6]. In other words, they are useful in the study of unstructured communities, those in which no species plays a distinctive role and no particular process can be singled out as having special importance in shaping the dynamics.

Most theoretical work so far has revolved around the so-called random generalized Lotka-Volterra model (gLV), where interactions are sampled independently from each other [3]. Despite recent qualitative agreement with some of the predictions of these models [7], many features of ecological communities remain elusive to ‘disordered ecology’[8]. Whether these mismatches are the signature of non-random structure in the interactions, and what elements should be added to the models to make them more realistic, are important research questions that have fostered developments in several directions [9, 10, 11, 12].

Yet, because current random models make stringent assumptions about the nature of the interactions, it can be hard to evaluate which of their conclusions are generic and which are model-specific, a key step to determine when structure is necessary to capture the features of real communities. The goal of this work is to propose a random model for unstructured communities that makes fewer assumptions and is therefore more general.

Inspired by physics, our approach consists in defining unstructured communities as those satisfying a natural symmetry. Intuitively, this symmetry must correspond to the idea that, in a community where no species stands out as performing any particular function, the relationship between species and interaction drivers can be randomized without altering the statistical properties of the model. For instance, if interactions are mediated by trait-matching (whereby species having similar traits compete more strongly with one another), we treat the traits of individual species as unimportant as long as correlations amongst traits are fixed (e.g. species that have one trait tend not to have another trait). By imposing this symmetry, the outcomes of the model become independent of the fate of any particular species, and are rather determined by which functions are being performed in the community[13, 14].

Surprisingly, this simple requirement is enough to simplify the dynamics up to a point where stationary regimes and stability conditions can be explicitly described. We show that our symmetry can accomodate a breadth of scenarios that is absent from current models without compromising the ability to make analytical calculations. Our results reveal that different scenarios lead to notable differences in observables such as Species Abundance Distributions (SAD) and responses to perturbations. All in all, this formalism provides a unified framework to understand what facets of species interactions shape the richness and stability of the community.

The paper is structured as follows: in Section II we recall the Lotka-Volterra equations and justify, in ecological terms, the symmetry underpinning the random matrix ensemble from which interactions will be sampled. In Section III we give an overview of the method to simplify the dynamics and analyze the properties of the Unique Fixed Point (UFP) solution. In Section IV we investigate how the extant community is related to the initial pool of species as well as the possibility of dynamical regimes beyond the UFP solution. We conclude in Section V.

II Model definition

II.1 The Lotka-Volterra equations

We consider a set of $S$ interacting species, labeled $i=1,\dots,S$ and model the temporal dynamics of their respective abundances, which we denote $x_{i}(t)$, according to the generalized Lotka-Volterra equations (gLV),

$$\frac{\mathrm{d}x_{i}}{\mathrm{d}t}=r_{i}x_{i}\left[1-u_{i}x_{i}+\sum_{j}A_{ij}x_{j}\right],$$ (1)

where $r_{i}$ and $u_{i}$ are respectively the bare growth rate and the self-regulation strength of species $i$. These parameters are intrinsic parameters of the species: $r$ depends on its adaptation to the environment, whereas $u$ encodes the competition between its members. To reduce the number of parameters, and consistently with our goal of modelling unstructured ecosystems, we will suppose that these quantities are the same for all species. In that case, $r$ can be absorbed by rescaling the time variable, so we will set $r=1$. On the other hand, we keep $u$ as a free parameter and will use it as a control parameter later on. All the following analytical work can be done without these simplifications.

Apart from its interaction with the environment, the net growth rate of the species depends on its interaction with the rest of the community, as represented by the matrix of interaction coefficients $A_{ij}$. This matrix encodes, in an abstract form, the information about all ecological processes that are relevant to understand the dynamics. As such, specifying the special characteristics of a community is tantamount to defining an interaction matrix. Tthe interaction matrix alone accounts for $S^{2}$ parameters. For $S\gg 1$, this is a prohibitively large number. Reducing the number of parameters is critical to recover a useful and tractable model.

The generalized Lotka-Volterra equations are the simplest non-linear model capturing two defining aspects of community ecology: exponential growth in the absence of interactions (due to reproduction) and stability of extinctions (that is, the impossibility for a species gone extinct to re-invade a community without migration from outside). Although seemingly simple, these two characteristics are in fact crucial to the understanding of many facets of the gLV dynamics [15, 16].

II.2 Our definition of structure

In order to define a matrix ensemble that represents unstructured communities, a fundamental question must first be tackled: how is structure reflected in the properties of the interaction matrix? This is a vast question, and many independent directions have been explored that lead to different results. Examples include hierarchical interactions [9], food webs [17] and sparse networks [11, 18]. Here, our starting point will be that structure manifests itself in the form of atypical directions in the interaction matrix. For instance, [12] has recently studied communities where interactions are constructed as the superposition of a random matrix and a non-random anisotropic matrix. The resulting non-random directions encode linear combinations of species abundances to which the dynamical system is particularly sensitive. Ecologically, these can be pictured, for instance, as ecosystemic functions (such as metabolite production or consumption) that drive the interactions and influence the dynamics.

This leads us to a definition of unstructured interactions as those that don’t possess such privileged directions. Mathematically speaking, we require the matrix ensemble to be invariant under orthogonal transformations,

$$\forall U\in\mathcal{O}_{S},\mathrm{d}\mathbb{P}\left(UAU^{T}\right)=\mathrm{d}\mathbb{P}(A),$$ (2)

where $\mathbb{P}$ is the measure associated to the matrix ensemble and $A$ is any matrix belonging to its support. Such models can be obtained by fixing an interaction matrix $A_{0}$ and rotating it by a random orthogonal matrix,

$$A=UA_{0}U^{T},~{}U\sim\mathrm{Haar}\left(\mathcal{O}_{S}\right),$$ (3)

the Haar measure being the uniform probability measure over the orthogonal group. Some ecological insight can be gained by inspecting what happens when $A_{0}$ is transformed into $A$ in such a way. Suppose that $A_{0}$ is decomposed into a sum of rank-one components that represent the different processes shaping interactions,

$$A_{0}=\sum_{\lambda}u_{\lambda}v_{\lambda}^{T},$$ (4)

which naturally arise in consumption/cross-feeding models [19, 20] (in which case $u_{\lambda}=v_{\lambda}$ is the preference vector associated to resource $\lambda$, i.e. the rate at which species feeds on/produces that resource), trait-based models [10] and function-based models (in which case $u_{\lambda}$ and $v_{\lambda}$ respectively encode the species’s response to and effect on function $\lambda$). Conjugating $A_{0}$ by the orthogonal matrix $U$ amounts to randomly rotating each of the vectors $u_{\lambda}$, $v_{\lambda}$, thereby randomizing the relationship between species and resources (resp. traits or functions). This randomization is consistent with the idea of no species playing a well-defined role in an unstructured community. Importantly, because rotations preserve scalar products, quantities such as $u_{\lambda}\cdot u_{\nu}$ are preserved, meaning that the overlaps between the important directions of the matrix are preserved, despite the directions themselves being randomized.

Mathematically, the effect of conjugating by a matrix $U$ can be pictured as a sort of mixing, whereby a new pool of species is defined whose interactions are a blend of the interactions of the previous pool. Our model is defined so that any such mixture is equally likely, as long as it preserves correlations between interaction drivers. In other words, it is statistically invariant under mixing. At the same time, due to the large number of species, typical realizations of the system will correspond to typical realizations of $\mathrm{Haar}\left(\mathcal{O}_{S}\right)$. This rules out eigenvector anisotropy (which may correspond, for instance, to a non-zero mean interaction or a division of the community into several groups), eigenvector localization (which rules out the presence of sub-communities made of strongly interacting species, as in some sparse models) and correlations between eigenvalues and eigenvectors.

Once randomization is done, the only remaining information about the original interaction matrix are its eigenvalues, which are common to the matrices of the ensemble thus defined. In particular, the spectrum encodes how the interactions were originally constructed.

In the simplest case, the spectrum is semi-circular, and the ensemble considered is the Gaussian Orthogonal Ensemble (GOE), where the $A_{ij}$ are (up to the symmetry constraint) independent and identically distributed Gaussian variables. For other spectral distributions, they display a diffuse correlation pattern, and are not Gaussian when considered as a whole. In the following section, we show how these diffuse correlation patterns have an important impact in the dynamics.

III Dynamical mean-field theory (DMFT)

III.1 General setting

In this section, we reduce the complex high-dimensional dynamics of (1) to a one-dimensional, albeit stochastic equation, assuming the prescription (3) for the interaction matrix. For readability purposes, we restrict our analysis to symmetric matrices, in which case formulas are more compact and the number of parameters is reduced. We provide a way to solve the dynamics for non-symmetric matrices, as well as other straightforward extensions of the model, in Supplementary Information (see SI Sec. LABEL:SI-sec:sparse & LABEL:SI-sec:asym).

Under the assumption that $A$ is symmetric, $A_{0}$ can be diagonalized in an orthonormal basis with real eigenvalues, so that without loss of generality $A_{0}$ can be taken to be diagonal. The matrix ensemble is then characterized by the distribution of its eigenvalues, whose probability density function (pdf) we denote $\rho$, so that the unstructured model can be divided in as many classes as there are pdfs. We further suppose that the support of $\rho$ is bounded and doesn’t scale with $S$, thus implying that our gLV model is in the so-called weak interaction regime, where interactions are weak when taken individually but strong at the collective level[21]. We emphasize that, at odds with traditional uses of orthogonal ensembles in Random Matrix Theory (RMT), the spectral distribution need not be a random object, but rather a parameter that encodes the way interactions are driven.

Since $A_{ij}$ are random variables, the dynamical trajectories will be random as well: our goal is to characterize the law of the random processes $x_{i}(t)$. To this end, we use the well-known path integral formalism [22], that starts by rewriting the probability of a given outcome $\vec{x}(t)$ as an integral over all possible realizations of the random matrix,

$$\mathrm{d}p\left(\vec{x}(0)\right)\int\mathrm{d}\mu(U)\prod_{i}\delta\left(\frac{\dot{x_{i}}}{x_{i}}-1+ux_{i}-\sum_{j}A_{ij}x_{j}\right),$$ (5)

where $\mu$ stands for the Haar measure over the orthogonal group, $p$ is the distribution of the initial conditions (which might be deterministic or random) and $\delta$ is a functional Dirac delta that enforces the trajectory to undergo the prescribed dynamics. The goal is then to simplify this expression up to a point where the probability distribution $\mathrm{d}\mathbb{P}\left[x_{i}(t)\right]$ can be recognized explicitly in a simple form. The strategy to do so is standard in disordered systems[23, 4], hence we mainly focus on the results, and relegate details to SI (see SI Sec. LABEL:SI-sec:dmft). The key step, performing the integral over $U$, relies on the following large deviation principle for the so-called HCIZ integral [24, 25],

$$\int\mathrm{d}\mu(U)\exp\left[\frac{S}{2}\mathrm{Tr}\left(UA_{0}U^{T}B\right)\right]\sim\exp\left[\frac{S}{2}\mathrm{Tr}G_{\rho}\left(B\right)\right],$$ (6)

where $A_{0}$ is the aforementioned diagonal matrix and $B$ can be any low-rank symmetric matrix of bounded spectrum. This introduces the function $G_{\rho}$, usually called free cumulant generating function, whose detailed definition can be found in Appendix B. In short, $G_{\rho}$ characterizes the shape of $\rho$. It will play a key role in determining the dynamics.

In practice, the averaging step is where the symmetry assumption is introduced in the model. This causes the path integral to no longer depend on all interaction coefficients but only on the spectrum $\rho$ and the two macroscopic quantities,

$$\begin{split}C(t,s)&=\braket{x(t)}{x(s)}=\frac{1}{S}\sum_{i}x_{i}(t)x_{i}(s)\\ K(t,s)&=\left\langle\text{\Large$\nicefrac{{\delta}}{{\delta\xi(s)}}$}\Big{|}x(t)\right\rangle=\frac{1}{S}\sum_{i}\frac{\delta x_{i}(t)}{\delta\xi_{i}(s)}.\end{split}$$ (7)

Both of them are descriptors of the dynamics. Indeed, $C$ is the auto-correlation of the abundance time-series, averaged over all species. $K$ is the average response of abundances at a given time $t$ w.r.t. an environmental perturbation at a previous time $s$ (see Eq. (8) for the definition of $\xi$). In (7) we have purposefully written them as scalar products in the space of abundances to reveal the fundamental role of orthogonal invariance: once the symmetry is averaged out, the path integral (and hence the dynamics) may only depend on quantities that are invariant under rotations. In other words, it may only depend on scalar products in the space of species abundances. By definition, these are collective quantities, defined as averages over all species. The fact that the path integral may only depend on them is what generates a regime of collective dynamics.

Upon subsequent simplifications, the path integral can finally be identified as generating the set of Stochastic Differential Equations (SDE),

$$\frac{\mathrm{d}x_{i}}{\mathrm{d}t}=x_{i}\left[1-ux_{i}+\int_{0}^{t}H(t,s)x_{i}(s)+\xi_{i}(t)\right],$$ (8)

Each species $i$ effectively follows the dynamics prescribed by this process independently of the other species. In particular, the effect of inter-species interactions has been replaced by two terms. First, a memory term driven by a memory kernel $H(t,s)$: it introduces a time-delayed self-regulation of the species on itself. Second, a Gaussian noise $\xi_{i}(t)$ characterized by a noise correlator $C_{\xi}(t,s)$. The fact that $\xi_{i}$ has Gaussian statistics can be traced back to the fact that, as was noted in the previous paragraph, scalar products are the only invariant quantity under orthogonal symmetry. In other words, the type of noise that is obtained in the effective equation is a direct consequence of the symmetries of the model.

Both $H$ and $C_{\xi}$ are deterministic functions that can be obtained from $C$ and $K$. In fact, $C$ and $K$ are also deterministic functions, since (8) defines them as averages over independent species. This implies a self-consistency relation between the effective drivers of single-species dynamics ($H$ and $C_{\xi}$) and the empirical many-species statistics ($C$ and $K$).

In principle, these self-consistency relations (which we relegate to SI Sec. LABEL:SI-sec:dmft for readability) would allow to solve for the four quantities. In practice, they include non-linear functional operators and are too intricate to be solved. Even a numerical solution would in fact prove far more costly than the already intensive solution proposed for GOE matrices in [26]. Hence, in the following we focus on stationary regimes, where most of this complexity fades away.

III.2 Stationary regimes

We define stationary regimes as those where one-point averages are independent of time and two-point averages depend only on time differences, i.e. $C(t,s)=C(t-s)$ etc. Previous works have leveraged the assumption of stationarity, even in conditions where the community doesn’t reach a fixed point [27]. Fortunately, time-translation invariance is enough to get rid of the most complex part of the general analysis, namely the fact that the self-consistent closure equations come in the form of functional operators. Indeed one can then rewrite the closure equations as scalar relationships between the Fourier modes of the different quantities. Let hatted letters denote the Fourier transforms of the corresponding quantities. Then we find (see SI Sec. LABEL:SI-sec:stat)

$$\begin{split}\hat{H}(\omega)&=G_{\rho}^{\prime}\left(\hat{K}(\omega)\right)\\ \hat{C_{\xi}}(\omega)&=\hat{C}(\omega)\frac{\mathrm{Im}\left[G_{\rho}^{\prime}\left(\hat{K}(\omega)\right)\right]}{\mathrm{Im}\left[\hat{K}(\omega)\right]}\end{split}$$ (9)

The first equation gives a mapping between the response function and the memory kernel. The second relates the fluctuations in population abundances to the DMFT noise fluctuations. In particular, the value of $\hat{K}(\omega)$ appears as a modulator to the relationship between abundance and noise fluctuations, suppressing or enhancing frequencies (but not de-phasing) depending on $G_{\rho}^{\prime}$. This is in contrast to the GOE case studied in [3], where the two are directly proportional. We interpret this as a signature of memory in the system due to the diffuse correlations present in the interaction matrix.

III.3 Fixed points

We now specialize the general solution even more to the case where the dynamics reach a fixed point. Compared to time-translation invariant regimes, this is obtained by further requiring that correlation functions are also time-independent. Loss of time dependence in the correlators ensures that all random and non-random processes converge to a fixed point. In particular, we define

$$\begin{split}x(t)&\rightarrow x\\ \xi(t)&\rightarrow\xi\\ C(t,t)&\rightarrow q\\ C_{\xi}(t,t)&\rightarrow z\\ \int_{0}^{t}\mathrm{d}sK(t,s)&\rightarrow\chi\\ \int_{0}^{t}\mathrm{d}sH(t,s)&\rightarrow\eta\end{split}$$ (10)

The first two are random quantities. In particular, $\xi$ is a centered Gaussian variable with variance $z$. On the other hand, $q$, $z$, $\chi$ and $\eta$ are defined as limits of deterministic processes and are therefore deterministic. More precisely, $q$ is the mean-square abundance at equilibrium, related to the variance of the abundance distribution; $\chi$ is the integrated response function, that is the response of the equilibrium abundance of a species to a permanent change in its carrying capacity; and $\eta$ is the integrated memory kernel, i.e. a modifier to the species self-regulation due to feedback through its interactions with other species. By setting the derivative to zero in (8), one obtains the relationship between $x$ and $\xi$,

$$x=\mathrm{max}\left(0,\frac{1+\xi}{u-\eta}\right).$$ (11)

The species abundance distribution is therefore a truncated Gaussian, extending known results for Wigner and Wishart matrices [3, 10].

The fixed point equations relating the deterministic quantities can then be obtained by taking the limit $\omega\rightarrow 0$ in the Fourier closure (9), rewriting $\eta=\hat{H}(0)$, $q=\hat{C}(0)$, $\chi=\hat{K}(0)$ and $z=\hat{C}_{\xi}(0)$, and making use of (11) to explicitly write down the resulting averages over $x$. Details of the derivation can be found in SI Sec. LABEL:SI-sec:fp. They read

$$\begin{split}\eta&=G^{\prime}_{\rho}\left(\chi\right)\\ z&=qG^{\prime\prime}_{\rho}\left(\chi\right)\\ q&=\langle x^{2}\rangle=\frac{z}{(u-\eta)^{2}}\omega_{2}\\ \chi&=\frac{\omega_{0}}{u-\eta}\end{split}$$ (12)

In writing them we have introduced the auxiliary variables $\omega_{k}=\int_{-1/\sqrt{z}}^{\infty}\mathcal{D}y\left(y+1/\sqrt{z}\right)^{k}$, where the integration is with respect to the standard Gaussian measure. In particular, $\omega_{0}$ is the fraction of surviving species at equilibrium. Eqs. (12) can also be obtained using the replica trick [28, 5], although in that case one doesn’t obtain the dynamics. Inspecting (12) immediately reveals that there are only two independent variables: for instance one may only keep the dynamical descriptors $q$ and $\chi$ or rather take the DMFT quantites $z$ and $\eta$. We keep the four equations to make them more readable.

Known results can be obtained by specializing $\rho$ to previously studied distributions. For instance, if $\rho$ is the semi-circle distribution with variance $\sigma^{2}$ then $G_{\rho}(x)=\sigma^{2}x^{2}/2$ and one obtains the well-known equations for the GOE case [5, 3, 4]. For other spectral distributions, $G_{\rho}$ contains higher-order cumulants that influence equilibrium properties. These are the signature of correlations between interactions at the level of the whole matrix, the GOE case being the very special case where those correlations are identically zero.

In the GOE case, the appearance of a Gaussian random variable in (11) can be understood from a cavity-like argument [3]: since interactions are (up to symmetry) independent from each other, and species are uncorrelated, the interaction term in (1) converges to a Gaussian in the limit of a large number of species. Moreover, the Central Limit Theorem (CLT) applies and gives the relation $z=\sigma^{2}q$. For other spectra, interactions are no longer independent, and hence the CLT cannot be applied. In spite of this, the interaction sum still falls in the basin of attraction of a Gaussian distribution, but the relation between $z$ and $q$ in Eq. (12) is modified to take correlations into account. Hence our ensemble operates in an intermediate regime between uncorrelated and strongly correlated interactions.

In the following, we illustrate the influence of $G_{\rho}$ by taking three reference classes of spectra. The first class is made of shifted and scaled Marcenko-Pastur distributions, which are associated to the Hebbian matrices studied in [10]. These matrices appear when interaction are driven by trait-matching or resource consumption. The second is made of uniform distributions over bounded intervals. The last is made of binary distributions, where the eigenvalues are evenly split between two values $\pm\lambda$. The last two examples are taken purely as theoretical tools for illustration purposes, without well-established theoretical scenarios in mind. The functions $G_{\rho}$ associated to each of these cases are given in Appendix B.

Fig. 1 shows the Species Abundance Distribution (SAD) for these three classes of spectra. Importantly, the spectra were tailored so as to all have the same variance. The differences in the results are therefore due to the higher-order cumulants encoded in $G_{\rho}$ or, in other terms, to the different shapes of each spectrum. By considering only the variance of the interactions and neglecting the shape of the spectra, a random model based on the GOE ensemble would miss these differences.

As a last remark, we note that our results build a bridge between the dynamical point of view of DMFT and the point of view of Random Matrix Theory, where one is often interested in computing the spectrum given a random matrix model. By leveraging our formalism, one can directly translate RMT outputs into Lotka-Volterra dynamics.

Figure 1: Influence of the spectral distribution on equilibrium abundances: The Species Abundance Distribution (SAD) is plotted for three different spectra. Blue: A shifted and scaled Marcenko-Pastur distribution, corresponding to Hebbian interactions. Orange: A uniform distribution. Green: A binary distribution, with only two unique values. The spectral distributions all have a variance of $0.4$ and the self-regulation is set to $u=1$, but the ensuing SADs are notably different.

IV Dynamical regimes and reduced interactions

In this section, we investigate how the dynamical process selects the interactions, and give conditions for the stability of the fixed point solution found in III.3. Finally, we discuss the possible dynamical regimes beyond the fixed point solution and their relationship to orthogonal symmetry.

IV.1 Spectrum of the reduced community

As the system reaches a fixed point, some species remain and some go extinct. We let $S^{\star}$ be the number of surviving species and $A^{\star}$ the $S^{\star}\times S^{\star}$ matrix obtained by restricting $A$ to surviving species. We call this the reduced interaction matrix. Even though species are, to leading order in $S$, independent from each other, subtle correlations exist that tend to favor species experiencing less competition. As such, $A^{\star}$ is not just a generic sub-matrix of the full interaction matrix. In particular, it has been shown that reduced matrices tend to exhibit in-row and in-column correlations [29]. On the other hand, (12) clearly shows that the properties of the fixed point, which are related to the properties of the reduced matrix, depend solely on $\rho$, thereby hinting at a simple relationship between the full spectrum and the assembled spectrum. In fact, one can check that (12) is invariant under the following transformation

$$\begin{split}G_{\rho}(x)&\longrightarrow\omega_{0}^{-1}G_{\rho}(\omega_{0}x)\\ \omega_{0}&\longrightarrow 1.\end{split}$$ (13)

Heuristically, we interpret this as the fact that the equilibrium can be seen both from the point of view of an interaction matrix with spectrum $\rho$ leading to a community with a surviving fraction $\omega_{0}$ or from the point of view of a community where all species survive but with spectrum $\nu$ such that $G_{\nu}(x)=\omega_{0}^{-1}G_{\rho}(\omega_{0}x)$. This is of course a heuristic interpretation, rather than a rigorous analytical argument, since $\omega_{0}$ is an outcome of the dynamical process, not a free parameter. Yet, numerical results under a variety of spectral distributions (semicircular, Marcenko-Pastur, uniform, Beta, and power distributions) confirm that the spectrum of $A^{\star}$ is indeed distributed according to $\nu$ so defined.

Figure 2: Influence of $S^{\star}/S$ on the assembled spectrum. The spectra of typical submatrices are shown for three classes of orthogonal models (Hebbian, Uniform and Binary, same as in Fig. 1) and different survival fractions. Blue curves correspond to the spectrum of the full matrix. Histograms correspond to numerical simulations and solid curves are predictions using Eq. (13). Both the Hebbian and the binary case have Dirac masses for $S^{\star}/S=1$, which are signalled by arrows. As $S^{\star}/S$ is decreased, spectral distributions drift away from their initial shape and towards a semi-circular distribution. In particular, the edges of the spectrum are softened to match the continuous edges of the semi-circular distribution.

Interestingly, we show (see SI Sec. LABEL:SI-sec:red-mat) that $G_{\nu}(x)=\omega_{0}^{-1}G_{\rho}(\omega_{0}x)$ also characterizes the spectral density of a generic submatrix of size $S^{\star}=\omega_{0}S$ of the original matrix. In other words, even though the reduced matrix has intricate correlations and doesn’t belong to an orthogonal ensemble anymore, the bulk of its spectrum matches what one would get by randomly sampling a set of $S^{\star}$ species and looking at their interactions, therefore bypassing the whole dynamical process. This was known for GOE matrices, where the reduced matrix also displays a semi-circular bulk up to a rescaling of the variance to $\sigma^{2}\omega_{0}$.

Yet, for general spectral distributions, the transformation induced by (13) is more subtle than a simple rescaling of the variance, and $A^{*}$ will not necessarily belong to the same class of matrices as $A$. In particular, given that higher-order cumulants are scaled down, communities with a lower survival rate tend to have a spectrum that ’interpolates’ between the original spectrum and the semi-circle distribution. This is shown in Fig. 2 for our three reference classes of spectra.

IV.2 Linear instability

The results in III.3 and IV.1 have been derived under the assumption that a unique stable fixed point (UFP) exists. This is trivially true when species don’t interact (that is, $\rho(x)=\delta(x)$) due to the self-regulation $u$ of the species. As the spectrum grows and the interactions become stronger, the stable fixed point solution remains until an instability threshold is crossed [3, 23]. In this section, we characterize the location at which this transition happens and the dynamics beyond the transition.

We probe the stability of the fixed point by introducing a small perturbation to the fixed point abundances (11) and studying the stability of the linearized SDE. This is equivalent to checking the stability of the Jacobian.
Let $\{x_{i}^{\infty}\}$ be an equilibrium configuration where each species is an independent realization of (11). Let $\delta x_{i}(0)$ be a small non-random perturbation to those equilibrium abundances. We let $x_{i}(t)=x_{i}^{\infty}+\delta x_{i}(t)$ be the ensuing dynamics. Because the system is no longer at equilibrium, the DMFT noise will also acquire a fluctuating part, which self-consistently stems from the fluctuations in the abundances. We denote it $\xi_{i}(t)=\xi_{i}^{\infty}+\delta\xi_{i}(t)$. Assuming that $\delta x_{i}(t)$ and $\delta\xi_{i}(t)$ remain small, we linearize the full dynamics (8) to get

$$\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}\delta x_{i}&=x_{i}^{\infty}\Big{[}-(u-\eta)\delta x_{i}(t)+\delta\xi_{i}(t)\Big{]}\\ &+\delta x_{i}(t)\Big{[}1-(u-\eta)x_{i}^{\infty}+\xi_{i}^{\infty}\Big{]},\end{split}$$ (14)

As per the equilibrium conditions, the second term in the r.h.s is zero for extant species and negative for extinct species. It can therefore never cause an instability. The first term is zero for extinct species. Therefore we may focus on this term and keep track of the surviving community only. In particular, instabilities may never originate in the extinct community.

Just as for the full dynamics, $\delta\xi_{i}(t)$ is a Gaussian process. Because the dynamical equation is now linear, this implies that $\delta x_{i}(t)$ is also a Gaussian process. In fact, it is convenient to take the Fourier transform of (14),

$$\left[i\omega+x_{i}^{\infty}\left(u-\eta\right)\right]\hat{\delta x_{i}}(\omega)=x_{i}^{\infty}\hat{\delta\xi_{i}}(\omega)+\delta x_{i}(0).$$ (15)

Instabilities are signalled by resonances in this expression. In particular, given that $\hat{\delta x_{i}}(\omega)$ is a Gaussian variable, we may look for divergences of its mean and its variance. After some manipulations (see SI Sec. LABEL:SI-sec:instab) we find that the mean is always stable and that the variance diverges when the following condition is met:

$$\left(u-\eta\right)^{2}-\omega_{0}G_{\rho}^{\prime\prime}\left(\frac{\omega_{0}}{u-\eta}\right)=0$$ (16)

Note how this condition is consistent with the spectral transformation (13). Leveraging this, we show (see SI Sec. LABEL:SI-sec:mat-inst) that the onset of instability happens precisely when

$$u-\lambda_{+}\left(A^{\star}\right)=0,$$ (17)

where $\lambda_{+}$ is the upper edge of the bulk of $A^{*}$. That is, the stability of the fixed point is equivalent to the stability of the reduced interaction matrix, similar to what happens in the Wigner case.

Interestingly, (16) not only implies a linear instability of the fixed point, but also a divergence of the mean square abundance $q$ (see SI Sec. LABEL:SI-sec:is-div). In other words, all fixed point instabilities lead to a Diverging Abundance (DA) phase, where a finite fraction of species’s abundances go to infinity in finite time. The relevance of such dynamical phase is discussed in Sec. IV.3.

Despite always leading to a DA phase, the transitions signalled by (16) can be of two different types, depending on the behavior of $u-\eta$. In the first case, $u-\eta\rightarrow 0$, which is associated to a divergence of the response function $\chi$, as defined in (12). In the second, $\eta<u$ and $\chi$ remains finite. The second case corresponds to the usual divergence observed in the Gaussian case, as it is indeed the only possible one under the assumption that $\rho$ is semi-circular.

Phenomenologically, finite $\chi$ transitions occur when the variance of the DMFT noise becomes too strong to be countered by the self-regulation of the species. At this point, some extensive subset of the community has strong-enough mutualistic interactions to produce a runaway of their abundances to infinity. In other words, the transition is caused by the feedback loop between $q$ and $z$ and causes $q,z\rightarrow\infty$. This always happens with a survival fraction $\omega_{0}=1/2$. On the other hand, infinite $\chi$ transitions are characterized by a vanishing effective self-regulation $u_{\mathrm{eff}}=u-\eta$. Since species are thus no longer subjected to competition from their own kind, they become able to grow exponentially even at large abundances, hence causing $q\rightarrow\infty$ again. Here, however, $z$ remains finite, i.e. the runaway is not caused by the feedback between $q$ and $z$. When this happens, more than half of the species survive, $\omega_{0}>1/2$.

This difference in behavior can be traced back to the shape of the reduced interaction matrix at the onset of instability. Whenever the spectrum of $A^{*}$ vanishes fast enough near the edge, a finite $\chi$ transition happens. This is the case, for instance, when the spectral density goes as $\nu(\lambda)\sim(\lambda_{+}-\lambda)^{\delta}$ near the edge $\lambda+$. Most continuous spectra (including the semi-circle law, for which $\delta=1/2$) will satisfy this requirement. On the other hand, if $\nu$ doesn’t vanish near the edge, an infinite $\chi$ transition occurs. The simplest case scenario here is that of a spectrum having a Dirac mass at the edge. This happens, for instance, if the reduced matrix has a Marcenko-Pastur spectrum with more surviving species than degrees of freedom [10]. The impact of the shape of the spectrum on the kind of transition stresses the importance of the transformation given by Eq. (13). Indeed, as the surviving fraction is decreased from one, the edges of the spectrum are softened, as was noted in Fig. 2. Hence $A^{*}$ may not have the same edge properties as $A$. In particular, if $A$ already has soft edges then so does $A^{*}$ and $\chi$ always remains finite at the transition (e.g. Fig. 2, center panel). For $A$ having hard edges, the question is whether these have been softened enough at the onset of instability. This can elucidated by going back to Eq. (12). It is clear that for a divergent $\chi$ transition to occur $G_{\rho}^{\prime}(x)\xrightarrow[x\rightarrow\infty]{}u$ needs to hold, in which case Eq. (16) automatically holds in the limit. This is a condition for the transition to be plausible. For it to occur, we also need to check that a finite $\chi$ transition isn’t reached beforehand. Simple manipulations (see SI Sec. LABEL:SI-sec:div-type) lead to the criterion

$$\lim_{x\rightarrow\infty}x^{2}G_{\rho}^{\prime\prime}(x)>\frac{1}{2}$$ (18)

When this is satisfied, the finite $\chi$ transition takes place, otherwise the divergent $\chi$ transition occurs.

IV.3 Global dynamical regimes

For ecological applications, the Diverging Abundances (DA) phase is meaningless. It stems from the fact that, as a phenomenological model, the Lotka-Volterra equations break down when there are sufficiently many positive interactions in the system. Further terms are needed to fix this pathological behavior. For example, a negative coupling to the total abundance (which would correspond to a non-zero mean in the interaction coefficients, biasing them towards competition) can displace (although not fully remove) the divergence of abundances. The interaction matrix then reads

$$A_{ij}=\frac{\mu}{S}+\tilde{A}_{ij},$$ (19)

for some $\mu<0$ and $\tilde{A}$ satisfying orthogonal invariance. In that case, even though the fixed point becomes unstable, the increased competition as abundances grows prevents the runaway to infinity beyond the linear instability. This region, where the fixed point solution is unstable but the system is non-divergent, has been proven to be a ’Multiple Attractors’ (MA) phase by previous works [5] in the particular case of a Wigner spectrum. There, the unique fixed point is replaced by an exponentially large number of fixed points[30].

Figure 3: Influence of the type of transition on the survival fraction. Biasing the community towards competition by means of a non-zero average interaction (see (19), here $\mu=-5$) removes the divergence of abundances but keeps the linear instability of the fixed point. This leads to a ’Multiple Attractors’ (MA) phase that is absent from the orthogonally-invariant model. Here we show the transition from the Unique Fixed Point (UFP) phase to the MA phase as the second moment $\sigma^{2}$ of the spectra is increased. Whether $\chi$ is finite or infinite at the instability can be determined from the behavior of the fraction of surviving species $\omega_{0}$ at and beyond the transition in this regularized model. Finite $\chi$ transitions (as for Wigner spectra, blue) are characterized by a continuous decrease of $\omega_{0}$. Conversely, $\chi\rightarrow\infty$ leads to a discontinuous drop (as for binary spectra, green). Dots are obtained from numerical simulations with $S=1000$ and solid lines are predictions adapting the fixed point equations (see III.3) to accomodate $\mu\neq 0$, in the region where the UFP is stable. By extending these predictions to the beginning of the MA phase (dashed lines), we see that the now-unstable UFP solution still captures the drop in $\omega_{0}$.

Orthogonal invariance is the reason why such MA regimes cannot be observed in our model. Indeed, a non-zero mean in the interactions corresponds to an anisotropy along the direction of the total abundance, given by the vector $(1,\dots,1)$. More generally, other low-dimensional anisotropies, such as a division of the community into groups, may lead to interesting phases such as bi-stability or, in the context of non-symmetric interactions, persistent fluctuations [15]. In any case, our model shows that these interesting dynamical regimes are a feature of structured communities, in the sense that some sort of symmetry breaking is required to get them. Even though we do not pursue this path in the present work, we note that all our calculations can be extended straightforwardly to account for any low-dimensional structure in the interaction matrix. Instead, we restrict ourselves to the simplest symmetry breaking, given by (19), and check how the two transitions identified in IV.2 influence the MA phase beyond the linear instability. Given that the difference between finite and infinite $\chi$ is a local feature (i.e., based on a local stability analysis), we expect the main features identified in the following to remain meaningful for more complicated structures.

The two types of transitions can notably be distinguished by the behavior of the fraction of surviving species $\omega_{0}$. As shown in Fig. 3, this fraction is continuous at the instability for finite $\chi$ transitions, whereas infinite $\chi$ transitions are characterized by a sudden drop at the onset of instability. This qualitative difference can be related to the spectral properties of $A^{\star}$, discussed in IV.2. Indeed, as the reduced interaction matrix becomes unstable, the only possibility for the system to recover a new stable fixed point is to kill enough species so as to recover a stable $A^{\star}$. Whenever the spectral density is singular at the edge (that is, whenever $\chi\rightarrow\infty$), a non-zero fraction of the eigenvalues vanish at once. Hence the need to kill at least as many species as there are zero eigenvalues. This explains the discontinuous drop in $\omega_{0}$. On the contrary, if the spectrum of $A^{\star}$ vanishes at the edge (that is, if $\chi$ remains finite), stability can be achieved by killing one single species at a time, since the eigenvalues of the matrix become unstable one at a time. In that case, $\omega_{0}$ decreases continuously.

V Discussion

Since their inception in theoretical physics, models with quenched disorder have proven a useful tool to understand complex systems in a wide range of fields. In ecology, they allow to bypass the need for finely-detailed biological information while providing non-trivial phenomenology. Yet, defining a random model requires deciding at which point in its construction randomness is added, and in which form. In turn, this implies strong assumptions that may both hinder the generality of the conclusions and limit the interpretability of the model.

In this work, we have provided a random framework that goes beyond previously studied cases and shown that a fundamental symmetry is enough to use the powerful tools of Dynamical Mean Field Theory and solve the dynamics in this more general scenario. Our symmetry, which requires the ensemble of interaction matrices to be invariant under orthogonal conjugation, is rooted in the idea that, in an unstructured ecosystem, the categorization of individuals into species is unrelated to the collective properties of the community. In other words, the community is defined by the macroscopic relations amongst its collective functions, rather than the association of specific species to specific functions. Due to this, the dynamics may not depend on individual interactions or species, but only in a few macroscopic quantities. In this collective regime, the trajectories of species statistically decouple, and each of them can be seen as an independent realization of an effective stochastic process.

In more fundamental terms, orthogonally-invariant ensembles can also be seen as a solid and more general class of matrices to recover some of the results that have stood out as being shared across previous disordered models. In particular, the fact that SADs are truncated Gaussians can be directly related to the fact that scalar products are the only quantities that are preserved by orthogonal conjugation. Similarly, we have shown that the spectrum of the reduced interaction matrix and that of the full matrix are related by a transformation rule that generalizes previous results for Wigner matrices. This transformation is relevant for at least two reasons. First, it tells us that the surviving community, despite being dependent on the realization of the interactions, interacts in a way that is not so dissimilar from any typical sub-community of the same size, i.e. any subset of species sampled at random rather than obtained from the ecological dynamics. In other words, even though the extinction process might create subtle correlations and biases in the reduced interactions, these are not strong enough to yield a reduced spectrum that differs from the typical expectation. Even though our analysis focuses on the bulk, our numerical simulations and analytical calculations show no trace of outlier eigenvalues in the dynamics . Second, it tells us that smaller assembled communities (that is, communities with a smaller survival rate), tend to be more similar to the ones obtained from the Wigner model of interactions. Indeed, we have shown that the reduced spectrum interpolates between the full one and the semi-circle distribution, due to a faster decay of higher-order moments.

Within our framework, only two dynamical phases are possible. Either the system possesses a unique fixed point (UFP phase), which is reached regardless of initial conditions, or abundances diverge in finite time (DA phase). The analytical fixed point equations then generalize the findings of previous models, and make explicit the relationship of the spectral distribution with the properties of the equilibrium and the location of the instability. As for the diverging phase, even though it is ecologically meaningless, we argue that the addition of low-dimensional structure can replace it with more interesting regimes such as the Multiple Attractors phase described in previous works. In any case, the added structure must break the symmetry of the model. This paves the way for further work towards understanding how different ways of breaking the orthogonal symmetry are related to different types of ecological structure and in turn generate different types of dynamical regimes.

Furthermore, we believe that understanding how the orthogonally-invariant model transitions from the UFP phase to the DA phase can still be relevant to capture the local characteristics of the transitions that take place in more structured models. In particular, we show that the two types of diverging transitions present in our model (identified by the behavior of the response function $\chi$) preserve their distinctive features once a competitive bias in the interactions is added to regularize the divergences. This is because, as the mean abundance increases, the interaction bias decreases the growth rate of all the species at the same time. This curbs the divergence in the non-linear regime, but doesn’t change the nature of the linear instability. Interestingly, we find that each kind of transition can thereby be mapped to a different behavior of the fraction of surviving species, which ties to the important ecological concept of coexistence. Infinite $\chi$ transitions lead to a sudden drop in the survival fractions, which could have a substantial impact in real ecological scenarios.

Crucially, all our calculations are dependent only on the spectrum of the interaction matrix. First, this establishes a powerful bridge with Random Matrix Theory, since spectra are typical outputs of RMT calculations. Additionally, by showing the equivalence between the instability of the fixed point solution and the instability of typical submatrices (which so far had only been addressed for Gaussian matrices [9]), we provide a solid justification for the use of RMT stability calculations. Second, and perhaps more importantly, our work paves the way towards a more thorough analysis of what characteristics of the spectrum are important to understand the dynamics. For example, we have shown that the behaviour of the spectrum at its edges plays a crucial role in determining the instabilities undergone by the system.

Finally, we wish to stress the potential benefits of grounding random models for community ecology on general considerations about the symmetries in species interactions. Here, orthogonal invariance has been used to translate the qualitative intuition that the dynamics of unstructured communities do not depend on the species’s identities. More generally, other symmetries may translate qualitative understandings of more complicated ecological scenarios, by differentiating between the characteristics that are relevant to the dynamics (and hence, must not be averaged out) and those that are not (and can be averaged out). The ensemble of matrices considered in Sec. IV.3 can be seen as a very simple breaking of orthogonal symmetry into a slightly less stringent symmetry, whereby the total abundance of the community is singled out as a direction of the interaction matrix, and hence becomes relevant to the understanding of the dynamics. Likewise, further breaking the $\mathcal{O}_{N}$ symmetry would give rise to partially structured matrices generalizing those in [12]. More generally, starting from an altogether different symmetry group could translate a different ecological intuition and lead to a different pattern in the dynamics. To what extent the link made here between symmetries and patterns can be generalized to other symmetry groups is an open question that we believe is worth addressing in future work.

Appendix A Numerical simulations

All numerical calculations where performed in Python. Simulations where done with Scipy’s solve_ivp function with the RK45 method. Interaction matrices were sampled by fixing the spectral distribution and sampling an orthogonal matrix from the Haar measure. Initial conditions were sampled uniformly between zero and 1. The solver was stopped after a fixed amount of time, generally $T=1000$. We used no immigration threshold. Classification between the two possible phases is simplified by the fact that in the divergence phase the divergence happens in finite time. Numerically, this causes Scipy’s integrator to stop before $T$ is reached. Hence, we classify a trajectory as divergent if the solver stops before $T$ is reached.

Numerical solutions for the fixed point equations where obtained by using an iterative scheme. The equations where rewritten in vector form as $\vec{X}=f(\vec{X})$ and the iteration ran as $\vec{X}_{n+1}=(1-r)\vec{X}_{n}+rf(\vec{X}_{n})$ for a fixed number of iterations. We used $r=0.25$. The initial condition was taken to be $\vec{X}_{0}=\vec{1}$. Except very close to the divergence transition, we find that $200$ iterations were more than enough to converge to the solution, although the examples of the main text were solved with $20000$ iterations.

Appendix B The function $G_{\rho}$

The reader will find a more comprehensive review of the following notions in [31]. For any random matrix $A$ with spectral distribution $\rho$, we define the Stieltjes transform of $\rho$ as

$$\mathfrak{g}(x)=\int\frac{\mathrm{d}\lambda\rho(\lambda)}{x-\lambda}.$$ (20)

$\mathfrak{g}$ is defined outside the support of the spectral distribution. For $x>\sup\mathrm{Supp}(\rho)$, $\mathfrak{g}$ is monotonically decreasing to zero, and is therefore invertible. We can then define the $R$-transform as

$$R_{\rho}(x)=\mathfrak{g}^{-1}(x)-\frac{1}{x}$$ (21)

The $R$-transform is defined over some interval $[0,\mathfrak{g}_{\mathrm{max}})$, where $\mathfrak{g}_{\mathrm{max}}=\lim_{x\rightarrow\sup\mathrm{Supp}(\rho)^{+}}\mathfrak{g}(x)$. One always has $R_{\rho}(0)=0$. We will make extensive use of the integral of the $R$-transform,

$$G_{\rho}(x)=\int_{0}^{x}\mathrm{d}yR_{\rho}(y)$$ (22)

In particular, $G_{\rho}$ is the ’free probabilites analogue’ of the cumulant generating function from free probabilities. Indeed, its series expansion is formally given by

$$G_{\rho}(x)=\sum_{n}\frac{\kappa_{n}}{n!}x^{n},$$ (23)

where $\kappa_{n}$ are the so-called free cumulants of $\rho$. For our needs the existence of the cumulants is guaranteed by the fact that we’ll only consider cases where $\rho$ is bounded. In particular, one may take the previous series expansion and use it to extend the definition or $R$ and $G$ to the complex plane. For the examples of the main text this gives

Type	$G_{\rho}^{\prime}$	$G_{\rho}^{\prime\prime}$
Wigner	$\sigma^{2}x$	$\sigma^{2}$
Hebbian	$\frac{\alpha x}{1+x}$	$\frac{\alpha}{(1+x)^{2}}$
Uniform	$\frac{\sqrt{\epsilon}}{\tanh(\sqrt{\epsilon}x)}-\frac{1}{x}$	$\frac{1}{x^{2}}-\frac{\epsilon}{\sinh^{2}(\sqrt{\epsilon}x)}$
Binary	$\frac{\sqrt{1+4\lambda x^{2}}-1}{2x}$	$\frac{\sqrt{1+4\lambda x^{2}}-1}{2x^{2}\sqrt{1+4\lambda x^{2}}}$