# Fluctuation spectra of large random dynamical systems reveal
hidden structure in ecological networks

Yvonne Krumbeck  Centre for Networks and Collective Behaviour, Department of
Mathematical Sciences, University of Bath, Bath, BA2 7AY, UK  Qian Yang
Beijing Institute of Radiation Medicine, Beijing 100850, P. R. China  George
W. A. Constable  Department of Mathematics, University of York, Heslington,
York, YO10 5DD, UK  Tim Rogers<EMAIL_ADDRESS>Centre for Networks and
Collective Behaviour, Department of Mathematical Sciences, University of Bath,
Bath, BA2 7AY, UK

###### Abstract

Understanding the relationship between complexity and stability in large
dynamical systems —such as ecosystems— remains a key open question in
complexity theory which has inspired a rich body of work developed over more
than fifty years. The vast majority of this theory addresses asymptotic linear
stability around equilibrium points, but the idea of ‘stability’ in fact has
other uses in the empirical ecological literature. The important notion of
_temporal stability_ describes the character of fluctuations in population
dynamics, driven by intrinsic or extrinsic noise. Here we apply tools from
random matrix theory to the problem of temporal stability, deriving analytical
predictions for the fluctuation spectra of complex ecological networks.
Excitingly, we are able to show how different network structures leave
distinct signatures in the spectrum of fluctuations, opening the door to a new
avenue of research in the age-old complexity/stability debate.

##  I Introduction

“Will a large complex system be stable?” asks the title of Robert May’s
seminal 1972 paper may_will_1972  that threw fuel on the fire of the
complexity-stability debate and popularised the use of random matrix theory
(RMT) in theoretical ecology. At first sight, answering this question with
mathematics seems impossible. The huge number of interactions in real-world
ecosystems hampers any attempt to create a precisely calibrated model, as the
challenge of measuring all necessary parameters seems insurmountable. What May
pointed out was that it might in fact not be necessary to know exact parameter
values; knowledge of their statistical distribution could be sufficient.
Combining the random model ecosystems proposed by Gardner and Ashby
gardner_connectance_1970 , with results of Ginibre ginibre1965statistical  in
RMT, May showed how complexity —measured in terms of the number of species and
the connectance of their interaction network— could decrease ecosystem
stability.

Although modelling ecosystems as using random community matrices has been
criticised james_constructing_2015 ; jacquet2016no  —with some arguing that
these serve best as a null model for ecosystem structure jr_unified_2010 —
this growing field has continued to provide insights into the mechanisms that
promote ecosystem stability. For instance, Allesina and Tang
allesina_stability_2012 ; allesina_stabilitycomplexity_2015  generalised the
community matrix model to account for different interaction types, elucidating
the important stabilising role of predator-prey interactions. We now have
quite a detailed view on the extent to which high-level ecosystem information
(such as trophic allesina2015predicting  or community grilli2016modularity
structures) can be incorporated into the RMT framework to give more accurate
predictions of the stability boundary.

The notion of stability referred to by May and these later works is that of
asymptotic linear stability of an equilibrium point. While this definition is
a natural mathematical choice, it belies the rich array of interpretations of
“ecological stability” present in the empirical ecological literature
donohue_2016 . In order to make clear the differences between these
interpretations, Grimm and Wissel grimm_babel_1997  created an inventory for
different types of stability measures used in ecology. Some of these in
particular are more attuned to the measures favoured by empirical ecologists.
One such measure is _temporal stability_ , often described as the constancy of
ecological variables relative to their mean, which is commonly used as an
indicator for ecological stability levins_coexistence_1979 ;
ives_stability_1999 ; lehman_biodiversity_2000 ; tilman_biodiversity_2006 ;
jiang_different_2009 ; loreau_species_2008 ; campbell_experimental_2011 ;
donohue_dimensionality_2013 . Recently Arnoldi et al. arnoldi_resilience_2016
employed the term _variability_ to describe the inverse of temporal stability
in a random community matrix model. In that work, they consider the scale of
response to persistent external (environmental) noise applied to an ecosystem.
While this is an important and useful measure, it does not capture anything of
the temporal characteristics of fluctuations in ecosystems which determine
temporal stability and are important precursors to linear and nonlinear
instabilities wiesenfeld1985noisy .

Figure 1: Fluctuation spectra in predator-prey systems. Left: a mixed
community of $N=1000$ species with randomly assigned predation relationships.
Right: a model ecosystem with two trophic levels— 200 predators and 800 prey
species. In both cases fluctuations are illustrated via a typical single-
species time series ($\xi(t)$, upper) and the mean power spectra
($\phi(\omega)$, lower), where circles are simulation results, and solid lines
our theory. Full details of all simulations are found in the Methods.

In this paper, we seek to bridge empirical and theoretical measures of
stability by developing a theoretical framework for the analysis of temporal
stability of ecosystems. Our key object of study is the _power spectral
density_ , a statistical measure that captures the frequency and amplitude of
noisy fluctuations in time series (see Fig. 1 for examples). Power spectra are
readily computed from empirical data and provide detailed information about
intrinsic fluctuations and (via the fluctuation dissipation theorem
kubo1966fluctuation ) response to external perturbations. Previous theoretical
studies of power spectra in low-dimensional systems have yielded important and
sometimes surprising results in fields including epidemiology, game theory,
and ecology alonso2007stochastic ; galla2009intrinsic ; mckane_predator-
prey_2005 . The method is particularly powerful in explaining the emergence of
persistent quasi-cyclic oscillations driven by noise. Until now, however, a
major limitation of this theory has been its restriction to models with very
small numbers of interacting elements for which the approach is analytically
tractable with existing methods. Here, by applying techniques from the
statistical physics of complex systems, we demonstrate the possibility to
derive exact analytic formulae for the power spectra of large random and noisy
dynamical systems.

We apply our new method to characterise the stochastic fluctuations of species
abundances in random Lotka-Volterra type ecosystem models. As a result, we
find that their temporal stability is universally characterised by a few key
parameters, including the proportion of predator-prey interactions and the
rate of population turnover. This result is a temporal analogue of the famous
Winger semi-circle law for random matrix eigenvalue distributions
wigner1958distribution  and points to the wide applicability of the theory we
develop.

Just as May’s RMT calculations are open to generalisations and refinements, so
too is our approach to temporal stability. We illustrate this flexibility of
the theory by incorporating trophic structure to our ecosystem models.
Subsequently, we discover a distinct signature of this type of structure: the
confinement of fluctuations to a fixed band of frequencies. Taken together,
these results raise the exciting prospect of being able to draw conclusions
about the internal structure of an ecosystem through the analysis of its
fluctuations. This approach could in the future be applied to the type of
dataset being generated in the emerging high-frequency temporal ecosystem
monitoring studies dornelas_2018 .

##  II Results

###  II.1 Interaction types determine fluctuation spectra

Our approach enables the computation of the power spectral density of
fluctuations in large random systems of a very general class; a full and
detailed derivation is given in the Methods. In the case of ecosystem
stability, the dynamical system in question is that describing the
interactions of different species. Many modelling choices are possible in this
context. For clarity we will focus here on an established modelling paradigm —
large Lotka-Volterra type ecosystems — and explore the extent to which the
nature of the species interactions affects the shape of the fluctuation
spectrum.

Following classic models of ecosystem dynamics, we consider $N$ species
occupying a domain of size $V$, writing $x_{i}(t)$ for the density of
individuals of species $i$ at time $t$. For large but finite $V$, standard
techniques (see Methods) allow us to describe the change of the species
densities by a set of stochastic differential equations (SDEs) :

$\frac{dx_{i}}{dt}=x_{i}\left(b_{i}+\sum_{j}^{N}\alpha_{ij}x_{j}\right)+\frac{1}{\sqrt{V}}\eta_{i}(t).$
(1)

Here, the coefficients $\alpha_{ij}$ for $i\neq j$ describe the interaction
between species $i$ and $j$, and the $\eta_{i}(t)$ are Gaussian noise term
with correlations
${\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\delta(t-t^{\prime})B_{ij}(\bm{x})}$.

We parameterise the model as follows. For simplicity (and to isolate the
effect of interaction types) we model each species as having the same birth
rate $b_{i}\equiv b$ and density dependent mortality rate
$\alpha_{ii}\equiv-b$. The other interaction coefficients $\alpha_{ij}$ are
chosen at random so that (i) each species interacts with an average of $c$
others, (ii) interactions have mean strength $\mathbb{E}|\alpha_{ij}|=\mu$ and
second moment $\mathbb{E}\alpha_{ij}^{2}=\sigma^{2}$, (iii) the correlation is
controlled by the symmetry parameter
$\gamma=\mathbb{E}[\alpha_{ij}\alpha_{ji}]/\sigma^{2}\in[-1,1]$. Crucially,
the full details of the distribution of the parameters $\alpha_{ij}$ are not
required, thanks to the universality of property of large random matrices
tao2008random ; tao2010random .

In the methods we show how these rates can be derived from a simple model of
pairwise species interactions which can be mutualistic, competitive, or
predatory. The frequency of predator-prey type interactions is tied to the
symmetry parameter $\gamma$. At $\gamma=-1$ all interactions are of the
predator-prey type ($\alpha_{ij}=-\alpha_{ji}$), at $\gamma=1$ only purely
mutualistic ($\alpha_{ij}=\alpha_{ji}>0$) or competitive
($\alpha_{ij}=\alpha_{ji}<0$) are present, and between these extremes there is
a random mix of interaction types.

With this choice of (random) parameters, each species density will fluctuate
around the scaled carrying capacity $x^{*}_{i}=1$, which, following
allesina_stabilitycomplexity_2015 , is stable provided
$b>\sqrt{c\sigma^{2}}(1+\gamma)$ (we refer to stone_feasibility_2018 ;
gibbs_effect_2018  for stability and feasibility of equilibrium states with
heterogeneous species abundance distributions). Around this fixed point,
species in the stochastic system (1) will exhibit approximately linear
fluctuations $\xi_{i}$, described by an Ornstein-Uhlenbeck process of the form

$\frac{d\bm{\xi}}{dt}=\bm{A}\bm{\xi}+\bm{\zeta}(t)\,.$ (2)

Here $\bm{A}$ is the Jacobian of (1), known as the community matrix in the
context of theoretical ecology, and $\bm{\zeta}$ is an $N$-vector of Gaussian
white noise with correlation matrix $\bm{B}=\bm{B}(\bm{x}^{*})$.

The power spectral density of fluctuations $\bm{\Phi}(\omega)$ is defined as
the Fourier transform of the covariance
${\mathbb{E}[\bm{\xi}(t)\bm{\xi}(t+\tau)^{T}]}$. For multivariate Ornstein-
Uhlenbeck processes one can show (see e.g. gardiner_stochastic_2009 ) that

$\begin{split}\bm{\Phi}(\omega)&:=\int_{-\infty}^{\infty}\mathrm{e}^{-\text{i}\omega\tau}\mathbb{E}[\bm{\xi}(t)\bm{\xi}(t+\tau)]d\tau\\\
&=(\bm{A}-i\omega\bm{I})^{-1}\bm{B}(\bm{A}^{T}+i\omega\bm{I})^{-1}.\end{split}$
(3)

In the Methods we show how to apply random matrix theory techniques to compute
the power spectral density, via a complex Gaussian integral representation of
the above matrix equation. The approach provides a general framework for
computing the fluctuations in large systems specified by random matrices
$\bm{A}$ and $\bm{B}$. We derive an expression for the mean-field power
spectral density $\phi(\omega)=\mathbb{E}[\bm{\Phi}_{ii}]$ in terms of the
resolvent function $r$. Specifically,

$\phi=|r|^{2}\frac{\mathbb{E}[B_{ii}]+2\mathrm{Re}(r)c\mathbb{E}[A_{ij}B_{ij}]}{1-|r|^{2}c\mathbb{E}[A_{ij}^{2}]}\,,$
(4)

where expectation is taken only over the non-zero entries of $\bm{A}$ and
$\bm{B}$, and $r\in\mathbb{C}$ solves the self-consistent equation

$\frac{1}{r}=-\mathbb{E}[A_{ii}]+i\omega-rc\mathbb{E}[A_{ij}A_{ji}]\,.$

This result holds for general random matrix models in which interaction
parameters are drawn from the same distribution for all species pairs. For the
present case of our Lotka-Volterra model, the community matrix coincides with
the interaction matrix (that is $A_{ij}=\alpha_{ij}$). In the methods we
derive the rules $\mathbb{E}[B_{ii}]=2b+c\mu$, $\mathbb{E}[A_{ij}B_{ij}]=0$
for the statistics of the noise correlation matrix. To get a sense for the
information contained in equation (4), we explore the result for several cases
with varying interaction structures.

Figure 2: Left: stable systems with $\gamma=0$ (purple $\diamond$) and
$\gamma=1$ (green $\circ$) have similar fluctuation spectra, with a higher
proportion of predator prey interactions associated with higher over all
excitation. Right upper: near instability a pole emerges at $\omega=0$,
corresponding to baseline drift in marginally stable systems. Right lower: In
log-log axes we see the different nature of the pole for mixed ($\gamma=0$,
purple $\diamond$) communities compared with those with only symmetric
interactions ($\gamma=1$, purple $\diamond$).

First consider a weak interaction limit where the diversity of species is
rather small, so that $\sigma^{2}\ll 1$. In this case we find a simple
Lorentzian spectral density:

$\begin{split}\phi(\omega)&=\frac{2b+c\mu}{b^{2}+\omega^{2}}+\mathcal{O}(\sigma^{2}).\end{split}$
(5)

Fluctuations of this type are indicative of a highly stable system in which
the balance of interaction types $\gamma$ has no influence. Next let us
consider a limit where the power spectral density shows significant
differences depending on the proportion of predator-prey interactions within
the community, in particular focussing on ecosystems that are near the
stability boundary.

In the case of a system with predator-prey interactions only, we have
$\gamma=-1$ and the ecosystem is stable for all positive birth rates $b$.
Expanding in small $b$ we find that fluctuations are of order $1/b$, but are
almost completely confined to a low-frequency window. If
$\omega^{2}<4c\sigma^{2}$ then

$\begin{split}\phi(\omega)&=\frac{2b+c\mu}{2c\sigma^{2}}\left[\frac{1}{b}\sqrt{4c\sigma^{2}-\omega^{2}}-\frac{c\mu}{2b+c\mu}\right]+\mathcal{O}(b),\end{split}$
(6)

with an order $1/\omega$ tail outside this range. Note that equation (6) has
the shape of a quarter-circle, to be viewed as a natural counterpart to the
Wigner semi-circle law in classical random matrix theory
wigner1958distribution . The result is illustrated in Fig. 1.

For a random mixture of interaction types with $\gamma=0$ no approximations
are necessary as (4) simplifies to

$\begin{split}\phi(\omega)&=\frac{2b+c\mu}{b^{2}-c\sigma^{2}+\omega^{2}}.\end{split}$
(7)

The stability boundary here is given by $b^{2}=c\sigma^{2}$. The above result
therefore implies the emergence of a $1/\omega^{2}$ divergence in the power
spectrum at low frequencies when such a system is close to instability (see
Fig. 2).

When only mutualistic or competitive interactions are present (i.e.
$\gamma=+1$), the full solution to (6) in this case is complicated, but for
stable systems appears qualitatively similar to the result (7) above. Near the
stability boundary, however, we find another behaviour. When
$b^{2}=4c\sigma^{2}$, we find

$\phi(\omega)=\frac{\sqrt{2}(4\sqrt{c\sigma^{2}}+c\mu)}{\sqrt{c\sigma^{2}(\sqrt{16c\sigma^{2}\omega^{2}+\omega^{4}}-\omega^{2}})}-\frac{(4\sqrt{c\sigma^{2}}+c\mu)}{2c\sigma^{2}}.$
(8)

In contrast to the previous case, this power spectrum exhibits a pole of order
$1/\sqrt{\omega}$ at low-frequency, followed by a $1/\omega^{2}$ tail at high
frequency (see Fig. 2).

Between these results, we are able to see how the proportion of predator-prey
interactions in an ecosystem leaves a signature in the fluctuation spectrum.
When predator-prey interactions are dominant, the shape of the spectrum is
pulled towards a quarter circle law (Fig. 1, left panel); when they are rare,
the low-frequency pole near instability changes its character (Fig. 2).

###  II.2 Species fluctuations exhibit strong heterogeneity

So far, we have considered only the mean power spectral density of
fluctuations. The cavity method technology employed in the derivation of (4)
can also yield detailed information about the fluctuation spectra of
individual species in an ecosystem model. Suppose one is interested in a focal
species $i$, and has data on the type and strength of interactions this
species has with others in its ecosystem, as well as an estimate of the large
scale ecosystem parameters such appearing in (4). It is possible to make use
of this data in a _single defect approximation_ (SDA) scheme in which one
considers the fluctuations of species $i$ when embedded in a large unknown
ecosystem.

In the Methods we show how to derive an SDA approximation
$\phi_{i}^{\text{SDA}}$ to the spectral density of fluctuations for species
$i$, given by the expression

$\phi_{i}^{\text{SDA}}=\frac{\phi^{\text{MF}}\sum_{i\sim
j}A_{ij}^{2}+2\mathrm{Re}(r^{\text{MF}})\sum_{i\sim
j}A_{ij}B_{ij}+B_{ii}}{|A_{ii}+i\omega+\bar{r}^{\text{MF}}\sum_{i\sim
j}A_{ij}A_{ji}|^{2}}\,,$ (9)

where $\phi^{\text{MF}}$ and $r^{\text{MF}}$ are the mean-field power spectrum
and resolvent obeying (4).

In Fig. 3 we compare the average power spectral density of all species with
the spectra of individual species as computed directly and via the SDA
approximation. It is immediate to notice that the mean-field power spectral
density is often not representative of individual species, which show
surprisingly strong heterogeneity in their fluctuation spectra. Another
interesting feature of these results is the presence of peaks in the power
spectral density away from zero for some species — this implies quasi-periodic
fluctuations in the population that are not observed in the ecosystem as a
whole. Finally we observe the curious feature that (for this model at least)
the total power of fluctuations appears approximately conserved, meaning that
those species which do not have large fluctuations at low frequencies are the
same as those with unusually large fluctuations at higher frequencies. At
present we do not have an intuitive explanation for this behaviour,
highlighting the richness of non-obvious information present in these complex
power spectral densities.

Figure 3: Heterogeneity in species fluctuations. Thin lines show power
spectral densities for $N=1000$ individual species in a predator-prey
ecosystem model, computed using the single-defect approximation (9). Thick
lines show the power spectral density for comparison the two species with
extremal fluctuations at low frequency, computed directly via (16); the dash
line shows the mean power spectral density.

###  II.3 Trophic structure induces fluctuation band-gap

In the above investigations, we have employed a simple ecosystem model in
which species interactions are assigned completely at random. In the past
fifty years of research into random matrix ecosystem models, far more
sophisticated and realistic models have been developed. Let us now illustrate
how our methods may be applied to more detailed models using the example of
ecosystems with explicit trophic structure.

Consider a large model ecosystem composed of $N_{x}$ predator species and
$N_{y}$ prey species, writing $x_{i}$ and $y_{j}$ for the density of predator
species $i$ and prey species $j$, respectively. With no prey-prey or predator-
predator interactions, the interaction structure is bipartite. Each predator
species has an extrinsic death rate $d$ and depends upon the consumption of
prey for reproduction. This consumption may come from a selection of $c_{x}$
different prey species for each predator, with $R_{ij}>0$ giving the predation
rate of predator $i$ on prey $j$. Conversely, each prey has birth rate $b$,
but is hunted by $c_{y}$ predators, where $N_{x}c_{x}=N_{y}c_{y}$. The SDEs
for the predator and prey densities are given by

$\begin{split}\frac{dx_{i}}{dt}&=x_{i}\left(-d-x_{i}+\sum_{j}^{N_{y}}R_{ij}y_{j}\right)+\frac{1}{\sqrt{V}}\eta_{i}(t),\\\
\frac{dy_{j}}{dt}&=y_{j}\left(b-y_{j}-\sum_{i}^{N_{x}}R_{ji}x_{i}\right)+\frac{1}{\sqrt{V}}\eta_{j}(t),\\\
\end{split}$ (10)

where $\eta_{i,j}(t)$ are Gaussian noise with
$\langle\eta_{i}(t),\eta_{j}(t^{\prime})\rangle=\delta_{ij}(t-t^{\prime})B_{ij}(\bm{x},\bm{y})$.
In the Methods we show how these equations (and the specific form of $B_{ij}$)
are derived from an individual-based model.

This model has an equilibrium state $(\bm{x}^{*},\bm{y}^{*})$, around which
linear-order fluctuations will occur, analogously to equation (2) above. We
compute a community matrix of the form

$\bm{A}=\left(\begin{matrix}-x^{*}\bm{I}&x^{*}\bm{R}\\\
-y^{*}\bm{R}^{T}&-y^{*}\bm{I}\end{matrix}\right),$ (11)

where the first ${i=1,\dots,N_{x}}$ rows and columns represent the predator
species, and the remaining ${j=N_{x}+1,\dots,N_{x}+N_{y}}$ rows and columns
correspond to the prey species. The noise matrix is given by

$\bm{B}=\left(\begin{matrix}2x^{*}(x^{*}+d)\bm{I}&-x^{*}y^{*}\bm{R}\\\
-x^{*}y^{*}\bm{R}^{T}&2y^{*}b\bm{I}\end{matrix}\right).$ (12)

In the Methods we develop a general approach to computing the power spectral
density of large random systems with bipartite structure such as this. The
method requires explicitly keeping track of the contributions associated to
each species group and their interactions. In the mean-field, this approach
delivers a set of equations (60) to be solved for the mean contributions to
the resolvent $r_{x},r_{y}$, and to the power spectrum, $\phi_{x},\phi_{y}$.
Fig 4 shows the shape of the power spectrum for predator and prey species in
this bipartite ecosystem. Surprisingly, we find that fluctuations are mainly
confined to a narrow window of frequencies, with a band-gap around zero.
Examination of the system (60) allows us to determine the window of excited
frequencies to be bounded by the critical frequencies

$\omega_{\pm}=\sqrt{bd\left(\frac{1}{c_{x}}+\frac{1}{c_{y}}\pm\frac{2}{\sqrt{c_{x}c_{y}}}\right)}.$
(13)

Figure 4: The exact power spectral densities (dots) for predator species in a
two trophic level model, computed numerically from Eq. (16), and the
corresponding mean power spectral density (solid line) obtained by solving Eq.
(60) (inset shows the prey species). Excited frequencies are confined to a
band (shaded) between critical frequences given in Eq. (13). The order $1/N$
peak at $\omega=1$ in the simulation result relates to the high-level
bipartite structure. Its location is predicted by a corresponding 2D system,
also shown here for comparison (dashed line).

The contrast between the power spectral density of this two-trophic-level
model to that of a mixed ecosystem with predator-prey interactions was
illustrated in Fig. 1. In Fig. 4 we show the spectrum in more detail,
highlighting the band of excited frequencies predicted by (13). The
suppression of fluctuations at low frequencies is analogous to spectral band-
gaps occurring in solid state physics. In the present context, it means that
observed time series will not exhibit baseline wander and can therefore be
considered to have a higher long-term temporal stability than the mixed
ecosystems explored above.

Comparisons between simulations and our analytical results shows another
interesting feature: an order $1/N$ disagreement at frequency $\omega=1$,
which is outside of the excited range. This can be explained by considering an
effective two-species model in which we consider only a single “average”
predator and prey pair. This 2D system has an eigenvalue pair with unit
imaginary part, giving rise to quasi-cycle behaviour as documented in
mckane_predator-prey_2005 . It is important to note that this contribution is
small relative to the rest of the spectral density, meaning that the bulk of
fluctuations of a structured ecosystem cannot be inferred from considering a
low-dimensional representative model.

##  III Discussion

In this study we have revisited the complexity–stability question in
theoretical ecology with a fresh perspective that develops a random matrix
theory approach to temporal stability as captured by the power spectrum of
fluctuations. We have applied our techniques to calculate analytic formulae
describing the mean power spectra of large Lotka-Volterra ecosystems. We find
the fluctuations are described by just a few key parameters: the mean,
variance and correlation of entries of the community matrix and noise
correlator. We further expanded the method to investigate the role of trophic
structures in determining temporal stability. In short, our approach allows us
to link the large scale statistical properties of interaction parameters with
the emergent fluctuations in species dynamics.

Amongst the many new results that this promising technique grants access to,
several findings from our investigation are worth recapping here. Part of the
power of random matrix theory is that it uncovers universal properties of
large classes of systems of a certain type. In the present case we find that,
in analogy to the famous Wigner semi-circle law, the details of the
distributions of matrix elements are unimportant beyond the handful of key
parameters identified. Our parameter $\gamma$, which controls the proportion
of predator-prey interactions (and hence the correlation of off-diagonal
elements in the interaction matrix) is found to be of crucial importance. At
one extreme, we find a semi-circular spectral profile, at the other we find a
pole at zero frequency which has either $1/\omega^{2}$ or $1/\sqrt{\omega}$
divergence, depending on the symmetry of interactions. When an explicit
trophic structure is incorporated into the model, it was necessary to adapt
our method to general bipartite networks. Here, we found a band-gap in the
power spectral density, implying that this high-level structure leads to
greater long-term temporal stability. Finally, going beyond these universal
results for the mean spectrum, we find a huge variability in fluctuations at
the individual species level. These are not visible within the bulk but are
captured by a single defect approximation, showing that some species may
exhibit quasi-cyclic oscillations even when no such signal is present in the
larger system.

In each model investigated in this paper, we have characterised stochastic
behaviour emerging on a macroscopic scale from the statistical properties of
the underlying microscopic interactions. We emphasise that, as illustrated in
our section on trophic levels, the fluctuations observed in large scale
systems with a certain structure are likely to be substantially different from
those of small scale models previously investigated. The models presented here
have been chosen for simplicity and clarity, and they only scratch the surface
of what can be achieved with this method. More realistic models might include
a consideration of e.g. the dynamical assembling process of ecosystems
galla2018dynamically , heterogeneous turnover rates gibbs2018effect , or
explicitly spatial models where spatio-temporal patterning may persist
mckane2014stochastic .

From an ecological perspective, it would be desirable to connect our
theoretical work to empirical investigations into ecosystem stability. In
contrast to the traditional viewpoint of asymptotic linear stability, our
methods directly address a fundamental empirical quantity – timeseries of
species abundance. It is exciting to consider the possibility of identifying
the signature of a certain interaction structure in the power spectrum of
oscillations in data gathered in the field. Before this vision can be
realised, some important further work is needed. Real ecosystems do not exist
in a vacuum – we must consider the role of the surrounding environment,
including interactions with external factors such as seasonal variation or
changing climate.

Finally, we wish to emphasise that —despite the ecological focus in this
paper— the models of the kind we analysed are ubiquitous in many different
fields, and the methods we use throughout the paper offer a general framework
for large dynamical systems with random variables. Models of large interaction
networks are also used in fields as varied as deep learning
pennington2017nonlinear , finance moran2019will , biochemistry
luo2007constructing  and neuroscience almog2019uncovering . All of these
systems depend on a high number of parameters that are often difficult to
measure empirically. Our method provides a possibility to compute the power
spectral density and gain new insight into the model, which relies only on
statistical meta parameters.

### Acknowledgements

T.R. and Y.K. gratefully acknowledge the support of the Royal Society.
G.W.A.C. thanks Leverhulme Trust for support through the Leverhulme Early
Career Fellowship.

##  IV Methods

###  IV.1 Power Spectral Density for a General Ornstein-Uhlenbeck Process

In the following we develop a new method to compute the power spectral density
of $N$-dimensional Ornstein-Uhlenbeck processes,

$\frac{d\bm{\xi}}{dt}=\bm{A}\bm{\xi}+\bm{\zeta}(t),$ (14)

where $\bm{\zeta}(t)$ is an $N$-vector of Gaussian white noise with
correlations
$\mathbb{E}[\bm{\zeta}(t)\bm{\zeta}(t^{\prime})^{T}]=\delta(t-t^{\prime})\bm{B}$.
The matrix $\bm{A}$ determines the mean behaviour of $\bm{\xi}$ and is
considered to be locally stable, i.e. all eigenvalues of $\bm{A}$ have
negative real part. Using the matrices $\bm{A}$ and $\bm{B}$ one can fully
determine the power spectral density of fluctuations for the Ornstein-
Uhlenbeck process.

We are interested in the case that the coefficients $A_{ij}$ and $B_{ij}$ are
derived from a complex network of interactions with weights drawn at random,
possibly with correlations. This framework encompasses a very general class of
models with a wealth of real-world applications including but not limited to
the ecological focus we have here. The method we describe exploits the
underlying network structure of $\bm{A}$ and $\bm{B}$ to deduce as self-
consistent scheme of equations whose solution contains information on the
power spectral density.

We start with the definition of the power spectral density $\bm{\Phi}(\omega)$
as the Fourier transform of the covariance
${\mathbb{E}[\bm{\xi}(t)\bm{\xi}(t+\tau)^{T}]}$ at equilibrium,

$\bm{\Phi}(\omega)=\int_{-\infty}^{\infty}\mathrm{e}^{-\text{i}\omega\tau}\mathbb{E}[\bm{\xi}(t)\bm{\xi}(t+\tau)]d\tau.$
(15)

From gardiner_stochastic_2009  on multivariate Ornstein-Uhlenbeck processes,
we know that the power spectral density can also be written in the form of the
matrix equation,

$\bm{\Phi}(\omega)=(\bm{A}-i\omega\bm{I})^{-1}\bm{B}(\bm{A}^{T}+i\omega\bm{I})^{-1}.$
(16)

In practice, this equation is difficult to use for large systems as large
matrix inversion is analytically intractable and numerical schemes are slow
and sometimes unstable. We take an alternative route by recasting (16) as a
complex Gaussian integral reminiscent of problems appearing in the statistical
physics of disordered systems. Our approach in the following is to treat
$\omega$ as a fixed parameter and drop the explicit dependence from our
notation. We begin by writing

$\begin{split}\bm{\Phi}(\omega)=&\frac{|\bm{A}-i\omega\bm{I}|^{2}}{\pi^{N}|\bm{B}|}\int_{\mathbb{C}}e^{-\bm{u}^{\dagger}\bm{\Phi}^{-1}\bm{u}}\bm{u}\bm{u}^{\dagger}\prod_{i=1}^{N}du_{i}\,.\end{split}$
(17)

Simplification of the integrand is achieved by unpicking the matrix inversion
in the exponent via a Hubbard-Stratonovich transformation stratonovich1957 ;
hubbard1959 . To this end we recast the system in the language of statistical
mechanics by introducing $N$ complex-valued ‘spins’ $u_{i}$ and $N$ auxiliary
variables $v_{i}$, with the ‘Hamiltonian’

$\mathcal{H}(\bm{u},\bm{v})=-\bm{u}^{\dagger}(\bm{A}-\text{i}\omega)\bm{v}+\bm{v}^{\dagger}(\bm{A}-\text{i}\omega)^{\dagger}\bm{u}+\bm{v}^{\dagger}\bm{B}\bm{v}\,.$
(18)

Introducing a bracket operator

$\langle\cdots\rangle:=\frac{\int_{\mathbb{C}}e^{-\mathcal{H}(\bm{u},\bm{v})}(\cdots)d\bm{u}d\bm{v}}{\int_{\mathbb{C}}e^{-\mathcal{H}(\bm{u},\bm{v})}d\bm{u}d\bm{v}}\,,$
(19)

we can obtain succinct expressions for the power spectral density
$\bm{\Phi}=\langle\bm{u}\bm{u}^{\dagger}\rangle$ as well as the resolvent
matrix
$\bm{\mathcal{R}}=(\text{i}\omega-\bm{A})^{-1}=\langle\bm{u}\bm{v}^{\dagger}\rangle$.
Thus we may write,

$\bm{\Phi}=\frac{1}{\mathcal{Z}}\int_{\mathbb{C}}e^{-\mathcal{H}(\bm{u},\bm{v})}\bm{u}\bm{u}^{\dagger}\prod_{i=1}^{N}du_{i}dv_{i}\,,$
(20)

where ${\mathcal{Z}=|\bm{A}-i\omega\bm{I}|^{2}/\pi^{2N}}$.

This construction may seem laborious at first, but it unlocks a powerful
collection of statistical mechanics tools, including the _cavity method_.
Originally, the cavity method has been introduced in order to analyse a model
for spin glass systems mezard_sk_1986 ; mezard_spin_1987 . Further
applications of the method include the analysis of the eigenvalue distribution
in sparse matrices rogers_cavity_2008 ; rogers_cavity_2009 ; metz2019spectral
. We will exploit the network structure in a similar fashion in order to
compute the power spectral density.

In our analysis, we find that it is convenient to split the Hamiltonian in Eq.
(18) into the sum of its local contributions at site $i$, $\mathcal{H}_{i}$,
and contributions from interactions between $i$ and $j$, $\mathcal{H}_{ij}$,

$\mathcal{H}=\sum_{i}\mathcal{H}_{i}+\sum_{i\sim j}\mathcal{H}_{ij}\,.$ (21)

These terms can be decomposed as
$\mathcal{H}_{i}=\bm{w}_{i}^{\dagger}\bm{\chi}_{i}\bm{w}_{i}$ and
$\mathcal{H}_{ij}=\bm{w}_{i}^{\dagger}\bm{\chi}_{ij}\bm{w}_{j}$, where we
introduce the compound spins $\bm{w}_{i}=(u_{i},v_{i})^{T}$ and transfer
matrices,

$\begin{split}\bm{\chi}_{i}&=\left(\begin{matrix}0&A_{ii}+i\omega\\\
-A_{ii}+i\omega&B_{ii}\end{matrix}\right)\,,\\\
\bm{\chi}_{ij}&=\left(\begin{matrix}0&A_{ji}\\\
-A_{ij}&B_{ij}\end{matrix}\right)\,.\end{split}$ (22)

Let us focus on the power spectral density of a particular variable $\xi_{i}$,
obtained from the diagonal element $\phi_{i}=\Phi_{ii}$. For this we compute
the single-site marginal $f_{i}$ by integrating over all other variables,

$f_{i}(\bm{w}_{i})=\frac{1}{\mathcal{Z}}\int_{\mathbb{C}}e^{-\mathcal{H}}\prod_{j\neq
i}d\bm{w}_{j}.$ (23)

Alternatively, $\phi_{i}$ can be obtained as the top left entry of the
covariance matrix
$\bm{\Psi}_{i}=\langle\bm{w}_{i}\bm{w}_{i}^{\dagger}\rangle$. We write the
covariance matrix as the integral,

$\bm{\Psi}_{i}=\int_{\mathbb{C}}f_{i}(\bm{w}_{i})\bm{w}_{i}\bm{w}_{i}^{\dagger}d\bm{w}_{i}\,,$
(24)

which could also be expressed in terms of a Gaussian integral,

$\bm{\Psi}_{i}=\frac{1}{\pi^{2}|\bm{\Psi}_{i}|}\int_{\mathbb{C}}e^{-\bm{w}_{i}^{\dagger}\bm{\Psi}_{i}^{-1}\bm{w}_{i}}\bm{w}_{i}\bm{w}_{i}^{\dagger}d\bm{w}_{i}\,.$
(25)

By comparing Eq. (24) and (25) we find that

$f_{i}(\bm{w}_{i})=\frac{1}{\pi^{2}|\bm{\Psi}_{i}|}e^{-\bm{w}_{i}^{\dagger}\bm{\Psi}_{i}^{-1}\bm{w}_{i}}\,.$
(26)

We now insert Eq. (21) into (23) and obtain,

$f_{i}(\bm{w}_{i})=\frac{1}{\pi^{2}|\bm{\Psi}_{i}|}e^{-\mathcal{H}_{i}}\int_{\mathbb{C}}\prod_{i\sim
j}\left(e^{-\mathcal{H}_{ij}-\mathcal{H}_{ji}}f_{j}^{(i)}d\bm{w}_{j}\right)\,,$
(27)

where we write $f_{j}^{(i)}$ for the cavity marginals,

$f_{j}^{(i)}(\bm{w}_{j})=\frac{1}{\mathcal{Z}^{(i)}}\int_{\mathbb{C}}e^{-\mathcal{H}^{(i)}}\prod_{k\neq
i,j}d\bm{w}_{k}\,.$ (28)

In essence, the above discussion amounts to organising the $2N$ integrals in
(20) in a convenient way, with the advantage of providing a simple intuition
for the role of the underlying network. The superscript $(i)$ is used to
indicate that the quantity corresponds to the cavity network where node $i$
has been removed.

Next we perform the integration in Eq. (27) and compare to the form in Eq.
(26). We thus obtain a recursion formula for the covariance matrix
$\bm{\Psi}_{i}$ and the cavity covariance matrices $\bm{\Psi}_{jl}^{(i)}$,

$\bm{\Psi}_{i}=\left(\bm{\chi}_{i}-\sum_{\begin{subarray}{c}i\sim j\\\ i\sim
l\end{subarray}}\bm{\chi}_{ij}\bm{\Psi}_{jl}^{(i)}\bm{\chi}_{li}\right)^{-1}\,,$
(29)

where the notation $i\sim j$ indicates that we sum over nodes $j$ connected to
node $i$. Unless there is some specific structure underlying the network, we
assume that most real world cases have a “tree-like” structure from the local
view point of a single node $i$. Hence, it is highly unlikely that the nodes
$j$ and $l$ are nearby in the cavity network where node $i$ is removed, and
thus $\bm{\Psi}_{jl}^{(i)}$ only gives non-zero contributions if $j=l$. We
therefore reduce Eq. (29) and obtain for the covariance matrix,

$\bm{\Psi}_{i}=\left(\bm{\chi}_{i}-\sum_{i\sim
j}\bm{\chi}_{ij}\bm{\Psi}_{j}^{(i)}\bm{\chi}_{ji}\right)^{-1}.$ (30)

Similarly, the cavity covariance matrix obeys the equation,

$\bm{\Psi}_{j}^{(i)}=\left(\bm{\chi}_{j}-\sum_{j\sim k,k\neq
i}\bm{\chi}_{jk}\bm{\Psi}_{k}^{(j)}\bm{\chi}_{kj}\right)^{-1}.$ (31)

Here we use that $\bm{\Psi}^{(i,j)}=\bm{\Psi}^{(j)}$ when the nodes $i$ and
$k$ are not connected. In other words, removing node $j$ from the cavity
network where node $i$ is missing, has the same effect as removing it from the
full network. The system (30) describes a collection of nonlinear matrix
equations that must be solved self-consistently.

For networks with high enough connectivity (and to good approximation even
with modest connectivity), the removal of a single node does not affect the
rest of the network. Hence the system (30) can be reduced to a smaller set of
equations approximately satisfied by the matrices $\bm{\Psi}_{i}$:

$\bm{\Psi}_{i}\approx\left(\bm{\chi}_{i}-\sum_{i\sim
j}\bm{\chi}_{ij}\bm{\Psi}_{j}\bm{\chi}_{ji}\right)^{-1}.$ (32)

The power spectral density $\phi_{i}$ can be obtained as the top left entry of
$\bm{\Psi}_{i}$.

In order to progress further, we now consider specific approximations that
help us compute the power spectral density. First we take a mean-field
approach in order to obtain the mean power spectral density for all nodes part
of the network; we then use the result for the mean field in order to compute
a close approximation to the local power spectral density of a single node.
Later, we adapt the method to partitioned networks where nodes belong to
different types of connected groups.

#### Mean Field

For the following we assume that all agents in the system behave the same on
average. In practice, the terms governed by self-interactions $A_{ii}$ are
drawn from the same distribution for all agents. Similarly, the terms
including $B_{ii}$ are governed by one distribution. Interaction strengths and
connections with other nodes in the network are also sampled equally for all
agents (we have explored a large Lotka-Volterra ecosystem as an example of
such a network). In the mean-field (MF) formulation we assume that the mean
degree and excess degree are approximately equal, and replace all quantities
in Eq. (30) and (31) with their average.
$\bm{\Psi}_{i}=\bm{\Psi}^{\mathrm{MF}}\,\forall i$. We then obtain the
following recursion equation,

$\bm{\Psi}^{\mathrm{MF}}=\left[\mathbb{E}[\bm{\chi}_{i}]-\mathbb{E}\left(\sum_{i\sim
j}\bm{\chi}_{ij}\bm{\Psi}^{\mathrm{MF}}\bm{\chi}_{ji}\right)\right]^{-1}.$
(33)

In order to solve this equation, we parameterise,

$\bm{\Psi}^{\mathrm{MF}}=\left(\begin{matrix}\phi&r\\\
-\bar{r}&0\end{matrix}\right),$ (34)

where the top left entry $\phi$ corresponds to the mean power spectral
density, and we introduce $r$ as the mean diagonal element of the resolvent
matrix $\bm{\mathcal{R}}$. Finally by inserting the ansatz (34) into equation
(33) we obtain,

$\begin{split}&\left(\begin{matrix}\phi&r\\\
-\bar{r}&0\end{matrix}\right)^{-1}=\left(\begin{matrix}0&\mathbb{E}[A_{ii}]+i\omega\\\
-\mathbb{E}[A_{ii}]+i\omega&\mathbb{E}[B_{ii}]\end{matrix}\right)\\\
&+c\left(\begin{matrix}0&\bar{r}\mathbb{E}[A_{ij}A_{ji}]\\\
-r\mathbb{E}[A_{ij}A_{ji}]&\phi\mathbb{E}[A_{ij}^{2}]+(r+\bar{r})\mathbb{E}[A_{ij}B_{ij}]\end{matrix}\right),\end{split}$
(35)

where $c$ is the average degree (i.e. number of connections) per node.
Moreover, the expectations in the second term are to be taken over connected
nodes $i\sim j$ (i.e. non-zero matrix entries).

From Eq. (35) above, we obtain the equations,

$\begin{split}\frac{\phi}{|r|^{2}}&=\mathbb{E}[B_{ii}]+c\left(\phi\mathbb{E}[A_{ij}^{2}]+2\mathrm{Re}(r)\mathbb{E}[A_{ij}B_{ij}]\right),\\\
\frac{\bar{r}}{|r|^{2}}&=-\mathbb{E}[A_{ii}]+i\omega-
cr\mathbb{E}[A_{ij}A_{ji}].\end{split}$ (36)

We solve the second equation in (36) for $r$ and write the mean power spectral
density in terms of $r$,

$\begin{split}\phi=&|r|^{2}\frac{\mathbb{E}[B_{ii}]+2c\mathrm{Re}(r)\mathbb{E}[A_{ij}B_{ij}]}{1-c|r|^{2}\mathbb{E}[A_{ij}^{2}]},\\\
r=&\frac{1}{2c\mathbb{E}[A_{ij}A_{ji}]}\left[-\mathbb{E}[A_{ii}]+i\omega\right.\\\
&\left.-\sqrt{(-\mathbb{E}[A_{ii}]+i\omega)^{2}-4c\mathbb{E}[A_{ij}A_{ji}]}\right]\end{split}$
(37)

This equation informs the first part of the results presented in the main
text.

#### Single Defect Approximation

The Single Defect Approximation (SDA) makes use of the mean-field
approximation for the cavity fields, but retains local information about
individual nodes. We parameterise similarly to Eq. (34) for a single
individual. Moreover, we replace all other quantities with the respective
mean-field approximation. Specifically, we obtain

$\begin{split}&\left(\begin{matrix}\phi_{i}^{\text{SDA}}&r_{i}^{\text{SDA}}\\\
-\bar{r}_{i}^{\text{SDA}}&0\end{matrix}\right)^{-1}=\left(\begin{matrix}0&A_{ii}+i\omega\\\
-A_{ii}+i\omega&B_{ii}\end{matrix}\right)\\\ &+\sum_{i\sim
j}\left(\begin{matrix}0&\bar{r}^{\text{MF}}A_{ij}A_{ji}\\\
-r^{\text{MF}}A_{ij}A_{ji}&\phi^{\text{MF}}A_{ij}^{2}+(r^{\text{MF}}+\bar{r}^{\text{MF}})A_{ij}B_{ij}\end{matrix}\right)\,.\end{split}$
(38)

We solve this equation for $\phi_{i}^{\text{SDA}},r_{i}^{\text{SDA}}$, which
delivers

$\begin{split}\frac{\phi_{i}^{\text{SDA}}}{|r_{i}^{\text{SDA}}|^{2}}&=\phi^{\text{MF}}\sum_{i\sim
j}A_{ij}^{2}+2\mathrm{Re}(r^{\text{MF}})\sum_{i\sim
j}A_{ij}B_{ij}+B_{ii}\,,\\\
r_{i}^{\text{SDA}}&=\left(A_{ii}+i\omega+\bar{r}^{\text{MF}}\sum_{i\sim
j}A_{ij}A_{ji}\right)^{-1}\,.\end{split}$ (39)

#### Partitioned Network

Previously we assumed that all nodes in a network are interchangeable in
distribution. However, many real-world applications feature agents with
different properties, imposing a high-level structure on the network. We
realise this by partitioning nodes into distinct groups that interact with
each other (see section Trophic Structure for a simple example).

In order to handle different connected groups we make use of the cavity method
as in Eq. (30) and (31). In particular, we split the sum in the second term on
the right-hand side of these equations into contributions from each group in
the partitioned network. Let $M$ denote the number of subgroups $V_{m}$ in a
partitioned network then we write,

$\begin{split}\bm{\Psi}_{i}&=\left(\bm{\chi}_{i}-\sum_{m}^{M}\sum_{\begin{subarray}{c}i\sim
j\\\ j\in
V_{m}\end{subarray}}\bm{\chi}_{ij}\bm{\Psi}_{j}^{(i)}\bm{\chi}_{ji}\right)^{-1}\,,\\\
\bm{\Psi}_{j}^{(i)}&=\left(\bm{\chi}_{j}-\sum_{m}^{M}\sum_{\begin{subarray}{c}j\sim
k\\\ k\in
V_{m}\end{subarray}}\bm{\chi}_{jk}\bm{\Psi}_{k}^{(j)}\bm{\chi}_{kj}\right)^{-1}\,.\end{split}$
(40)

Similar to the previous sections we replace all quantities with a mean-field
average $\bm{\Psi}_{m}^{\mathrm{MF}}$, but for each group separately. Hence we
obtain $M$ equations of the form

$\bm{\Psi}_{i}^{\mathrm{MF}}=\left[\mathbb{E}[\bm{\chi}_{i}]-\mathbb{E}\left(\sum_{m}^{M}\sum_{\begin{subarray}{c}i\sim
j\\\ j\in
V_{m}\end{subarray}}\bm{\chi}_{ij}\bm{\Psi}_{m}^{\mathrm{MF}}\bm{\chi}_{ji}\right)\right]^{-1}.$
(41)

In order to compute the mean power spectral density for different groups
separately, we use a parameterisation as in Eq. (34) for each group. Therefore
we have,

$\bm{\Psi}_{m}^{\mathrm{MF}}=\left(\begin{matrix}\phi_{m}&r_{m}\\\
-\bar{r}_{m}&0\end{matrix}\right),$ (42)

for all $m=1,\dots,M$. This delivers $2M$ equations to solve for all $r_{m}$
and $\phi_{m}$. Numerically this is straight forward, although algebraically
long-winded for the general case. However, the equations simplify for special
cases. In section Trophic Structure we demonstrate this method for a bipartite
network where a lack of intra-group interactions simplifies the analysis.

###  IV.2 Large Lotka-Volterra Ecosystem

#### Model Description

First, we define the framework for a general Lotka-Volterra ecosystem with $N$
species and a large but finite system size $V\gg 1$. Let $X_{i}$ denote the
number of individuals and $x_{i}=X_{i}/V$ the density of species
$i=1,\dots,N$. We start from the following set of reactions that define the
underlying stochastic dynamics of the system:

$\begin{split}&X_{i}\overset{b_{i}}{\longrightarrow}2X_{i}\quad\text{(birth)}\\\
&2X_{i}\overset{R_{ii}}{\longrightarrow}X_{i}\quad\text{(death)}\\\
&X_{i}+X_{j}\overset{R_{ij}}{\longrightarrow}\begin{cases}2X_{i}+X_{j}&\mathrm{(mutualism)},\\\
X_{i}&\mathrm{(competition)},\\\
2X_{i}&\mathrm{(predation)}.\end{cases}\end{split}$ (43)

The self-interactions are governed by the birth rate ${b_{i}>0}$ and density-
dependent mortality rate ${R_{ii}>0}$. Furthermore, we define three
interaction types between species $i$ and $j$, namely mutualism, competition
and predation. In the case of mutualistic interactions, both species benefit
from each other, whereas competition means that both species have a higher
mortality rate, depending on the density of the other species. For predator-
prey pairs, one predator species benefits from the death of a prey species.
The predator and prey species are chosen randomly, such that species $i$ is
equally likely to be a predator or prey of species $j$.

With probability $P_{c}$ we assign an interaction rate $R_{ij}>0$ to the
species pair $(i,j)$, and with probability $1-P_{c}$ there is no interaction
between species $i$ and $j$ (i.e. $R_{ij}=0$). In other words, each species
has on average $c=NP_{c}$ interaction partners. The reaction rates are
considered to be i.i.d. random variables drawn from a half-normal distribution
$|\mathcal{N}(0,\sigma^{2})|$, where we write for the mean reaction rate
$\mu=\mathbb{E}[R_{ij}]=\sigma\sqrt{2/\pi}$ and raw second moment
$\sigma^{2}=\mathbb{E}[R_{ij}^{2}]$. For each interaction pair, the
interaction type is chosen such that the proportion of predator-prey pairs is
$p\in[0,1]$, and all non-predator-prey interactions are equally distributed
between mutualistic and competitive interactions (i.e. the overall proportion
of mutualistic/competitive interactions is $1/2(1-p)$). Lastly, we define the
symmetry parameter $\gamma=1-2p$, where $\gamma=-1$ if all interactions are of
predator-prey type ($p=1$), and similarly $\gamma=+1$ if there are no
predator-prey interactions ($p=0$). In a mixed case where predator-prey and
mutualistic/competitive interactions have equal proportion ($p=1/2$), we have
$\gamma=0$. Later we will see that $\gamma$ is equivalent to the correlation
of signed interaction strengths.

In the limit $V\rightarrow\infty$, the dynamics of the species density $x_{i}$
obey the ordinary differential equations,

$\frac{dx_{i}}{dt}=x_{i}\left(b_{i}+\sum_{j}^{N}\alpha_{ij}x_{j}\right),$ (44)

where $\alpha_{ij}$ are the interaction coefficients with
${|\alpha_{ij}|=|\alpha_{ji}|=R_{ij}}$. The signs of the interaction
coefficients are determined by the type of interaction between species $i$ and
$j$. For mutualistic interactions we have ${\alpha_{ij}=\alpha_{ji}>0}$, and
${\alpha_{ij}=\alpha_{ji}<0}$ for competitive interactions. In the case of
predator-prey interactions the coefficients have opposite sign
${\alpha_{ij}=-\alpha_{ji}}$. Hence the symmetry parameter as described above
is given by the correlation of interaction coefficients
$\gamma=\mathbb{E}[\alpha_{ij}\alpha_{ji}]$. Furthermore, in order to ensure
bounded species densities, we require negative self-interactions
$\alpha_{ii}=-R_{ii}<0$.

If species live in isolation (i.e. when ${\alpha_{ij}=0\,\forall i\neq j}$),
we see that the densities approach the “effective” carrying capacity
$K_{i}=-b_{i}/\alpha_{ii}$. For the following computations we consider a large
Lotka-Volterra system. Since we are only interested in the effects of
interactions between species, we assume that all self-interactions are
approximately equal. Thus we write for the birth rate $b_{i}=b$ and mortality
rate $\alpha_{ii}=-b$. This gives the effective carrying capacity $K=1$ for
all species.

The fixed point $\bm{x}^{*}$ at the deterministic equilibrium state is given
by,

$x_{i}^{*}=1+\sum_{j\neq i}\alpha_{ij}x_{j}^{*}.$ (45)

We assume a random mixture of mutualistic and competitive interactions with
equal proportions, and therefore the interaction coefficients $\alpha_{ij}$
have zero mean ($\forall i\neq j$). Furthermore, we postulate that for large
ecosystems where $N\rightarrow\infty$, the equilibrium state
$x_{i}^{*}=\mathbb{E}[x_{i}^{*}]\equiv x^{*}$. Hence we obtain the expected
equilibrium density $x^{*}=1$ for all species $i$ (note that the following
computations are valid for any known fixed point $\bm{x}^{*}$, and our
assumptions are for mathematical simplification only). This assumption allows
us to write the Jacobian matrix for a linearisation around the equilibrium
state, with elements,

$\begin{split}J_{ii}|_{\bm{x}=\bm{x}^{*}}&=\alpha_{ii}=-b,\\\
J_{ij}|_{\bm{x}=\bm{x}^{*}}&=\alpha_{ij}.\end{split}$ (46)

In other words, the community matrix of a large Lotka-Volterra system as
described above has the same form as the interaction matrix, i.e.
$A_{ij}=\alpha_{ij}$. The local stability of such community matrix $A$ is
given by the elliptic law allesina_stability_2012 ;
allesina_stabilitycomplexity_2015 . It states that with high probability all
eigenvalues of the random matrix $A$ are distributed on an ellipse in the
complex plane, centered at $(-b,0)$ on the real axis. Thus for a stable matrix
we require all eigenvalues to be negative, and hence the horizontal semi-axis
of the ellipse determines the allowed range for the centre. It follows the
stability criterion,

$\sqrt{c\sigma^{2}}(1+\gamma)<b,$ (47)

with the average number of connections $c$ per species, and the correlation
$\gamma=\mathbb{E}[A_{ij}A_{ji}]$. For a random community matrix (i.e.
$\gamma=0$), we recover the stability criterion that has been proven by May
may_will_1972 , $\sqrt{c\sigma^{2}}<b$. If $\gamma<0$, where the proportion of
predator-prey type interactions is larger, the horizontal semi-axis of the
ellipse becomes smaller. In other words, the stability criterion relaxes for
predator-prey interactions. For $\gamma=-1$ (i.e. $A_{ij}=-A_{ji}\forall
i,j$), all interactions are of predator-prey type and all eigenvalues become
purely imaginary. Therefore the stability criterion becomes $0<b$, as the
ellipse stretches vertically into the imaginary plane. The opposite is true
for mutualistic/competitive interactions (i.e. $\gamma=+1$), where eigenvalues
are distributed on an ellipse with large horizontal radius along the real
axis. Thus it is more likely that some eigenvalues have positive real part and
the system destabilises. We choose the parameter $b$ for each case, such that
the stability criteria are fulfilled.

For a large but finite system size $V$, we write the stochastic differential
equations,

$\frac{dx_{i}}{dt}=x_{i}\left(b+\sum_{j}^{N}\alpha_{ij}x_{j}\right)+\frac{1}{\sqrt{V}}\eta_{i}(t),$
(48)

where $\eta_{i}(t)$ are Gaussian random variables with
$\langle\eta_{i}(t),\eta_{j}(t^{\prime})\rangle=\delta(t-t^{\prime})B_{ij}$.
The noise matrix $\bm{B}$ can be obtained from the reactions that determine
the process. The diagonal elements are given by the self-interactions and
total interaction from all other species, and the off-diagonal elements depend
on the type of interaction between species $i$ and $j$. Here we assume that
only predator-prey type interactions contribute to the covariance of species.
Therefore, we write

$\begin{split}B_{ii}&=x^{*}\left(b+\sum_{j=1}^{N}R_{ij}x^{*}\right)\approx
2b+c\mu,\\\ B_{ij}&=\begin{cases}-R_{ij}x^{*}x^{*}=-R_{ij}&\text{if
}\alpha_{ij}=-\alpha_{ji},\\\ 0&\text{else.}\end{cases}\end{split}$ (49)

Recall that $\mu$ is given as the mean reaction rate
$\mu=\mathbb{E}[R_{ij}]=\sigma\sqrt{2/\pi}$.

Finally, we linearise around the fixed point to obtain a new equation for the
fluctuations $\bm{\xi}=\sqrt{V}(\bm{x}-\bm{x}^{*})$,

$\frac{d\bm{\xi}}{dt}=\bm{A}\bm{\xi}+\bm{\zeta}(t),$ (50)

which has the form of an Ornstein-Uhlenbeck process as defined in Eq. (14).

#### Computing the Power Spectral Density

Let us now compute the mean power spectral density $\phi$ of the process
described above using Eq. (37) as starting point. We replace the necessary
quantities that we obtain from the community matrix $\bm{A}$ and noise matrix
$\bm{B}$ as defined in the previous section. In particular, we have the
expected diagonal elements of the community matrix $\mathbb{E}[A_{ii}]=-b$,
and the noise matrix $\mathbb{E}[B_{ii}]=2b+c\mu$. Moreover the raw second
moment of the non-zero interactions is given by
$\mathbb{E}[A_{ij}]=\sigma^{2}$ and the correlation
$\mathbb{E}[A_{ij}A_{ji}]=\gamma\sigma^{2}$. We use that
${\mathbb{E}[A_{ij}B_{ij}]=0\,\forall i,j}$ since the off-diagonal elements of
the noise matrix are only non-zero if there is a predator-prey interaction
between species $i$ and $j$. However, the elements of $A_{ij}$ have opposite
signs in the case of predator-prey pairs and thus sum to zero.

Plugging in these quantities into Eq. (37) we obtain,

$\begin{split}\phi=&|r|^{2}\frac{2b+c\mu}{1-|r|^{2}c\sigma^{2}}\,,\\\
r=&\frac{1}{2c\gamma\sigma^{2}}\left[b+i\omega-\sqrt{(b+i\omega)^{2}-4c\gamma\sigma^{2}}\right]\,.\end{split}$
(51)

In the main text we explore the theoretical ecological consequences of this
result.

###  IV.3 Trophic Structure Model

Figure 5: Example of a bipartite predator-prey network. Here we have $n_{x}=4$
predators with $c_{x}=3$ prey each, and $n_{y}=6$ prey with $c_{y}=2$
predators each. The bold lines illustrate the connections to the focal node
before and after removing a node from the network. Predator nodes (black) only
have local contributions from prey nodes (white) and vice versa. In the mean-
field approximation for the power spectral density, these contributions are
replaced by the average of each group.

#### Model Description

In the following we define a model analogous to the one described in the
previous section. For a large but finite system size $V$ we write the model in
terms of a stochastic process. Previously we allowed for different types of
interactions, however in this model we only focus on predator-prey
interactions. More specifically, the interaction network is partitioned into
$N_{x}$ predator species and $N_{y}$ prey species, where $N=N_{x}+N_{y}$ is
the total number of species. We assume that predators only interact with prey
and vice versa (i.e. we assume no inter-species interactions within the groups
of predators or prey) as illustrated in Fig. 5. Moreover, each predator and
prey species interacts with themselves (density-dependent mortality).

In the previous model we assigned the same birth rate to all species in the
ecosystem. Here we assume that predators decline at rate $d$ in absence of
prey, and $b$ is the birth rate of prey species. For simplicity, we assume
that $d,b$ are fixed quantities, equal for all predators and prey
respectively. Furthermore $R_{ij}$ is the interaction rate between predator
$i$ and prey $j$. Each predator species has a fixed number of prey $c_{x}$ and
each prey species has a fixed number of predators $c_{y}$, such that
$N_{x}c_{x}=N_{y}c_{y}$. The parameters $c_{x},c_{y}$ can be interpreted as
outgoing degrees of predator and prey nodes respectively. Connections between
predators and prey are then wired randomly. The interaction strength is set to
$R_{ij}=\alpha$, and considered equal for all predator-prey interactions
(analogous to a mean reaction rate). Where there is no interaction between
species, the interaction rate is simply set to zero. Note that this means that
the total sum of interaction strength is constant $\alpha c_{x}$ and $\alpha
c_{y}$ for all predator and prey species respectively. In contrast to the
previous model, now only the network structure contributes to the randomness
of the system.

Let $x_{i}$ denote the density of predator species ${i=1,\dots,N_{x}}$, and
$y_{j}$ the density of prey species ${j=1,\dots,N_{y}}$. In the deterministic
limit where $V\rightarrow\infty$ we then write the following ODEs,

$\begin{split}\frac{dx_{i}}{dt}&=x_{i}\left(-d-x_{i}+\sum_{j=1}^{N_{y}}R_{ij}y_{j}\right),\\\
\frac{dy_{j}}{dt}&=y_{j}\left(b-y_{j}-\sum_{i=1}^{N_{x}}R_{ji}x_{i}\right).\end{split}$
(52)

Given the fixed number of connections $c_{x},c_{y}$ and interaction strength
$\alpha$, we can simplify the ODEs to two equations for the average predator
and prey densities,

$\begin{split}\frac{dx}{dt}&=x\left(-d-x+c_{x}\alpha y\right),\\\
\frac{dy}{dt}&=y\left(b-y-c_{y}\alpha x\right).\end{split}$ (53)

In the limit of large $N$, the equilibrium state of the system converges to
the average quantities obtained from this reduced form. The biologically
relevant equilibrium states for this system are given by the trivial fixed
points ${(x^{*},y^{*})=(0,0),(0,b)}$, and the non-trivial fixed point,

$\begin{split}x^{*}=\frac{c_{x}\alpha b-d}{c_{x}c_{y}\alpha^{2}+1},\\\
y^{*}=\frac{c_{y}\alpha d+b}{c_{x}c_{y}\alpha^{2}+1}.\end{split}$ (54)

Next, we write the Jacobian matrix for a linearisation around the non-trivial
fixed point. The community matrix takes the form,

$\bm{A}=\left(\begin{matrix}-x^{*}\bm{I}&x^{*}\bm{R}\\\
-y^{*}\bm{R}^{T}&-y^{*}\bm{I}\end{matrix}\right),$ (55)

where the first ${i=1,\dots,N_{x}}$ rows and columns represent the predator
species, and the remaining ${j=N_{x}+1,\dots,N_{x}+N_{y}}$ rows and columns
correspond to the prey species.

For a large but finite system size $V$ we write the corresponding stochastic
differential equations,

$\begin{split}\frac{dx_{i}}{dt}&=x_{i}\left(-d-x_{i}+\sum_{j}^{N_{y}}R_{ij}y_{j}\right)+\frac{1}{\sqrt{V}}\eta_{i}(t),\\\
\frac{dy_{j}}{dt}&=y_{j}\left(b-y_{j}-\sum_{i}^{N_{x}}R_{ji}x_{i}\right)+\frac{1}{\sqrt{V}}\eta_{j}(t),\\\
\end{split}$ (56)

where $\eta_{i,j}(t)$ are Gaussian noise with
$\langle\eta_{i}(t),\eta_{j}(t^{\prime})\rangle=\delta_{ij}(t-t^{\prime})B_{ij}$.
The noise matrix is given by the self- and total interactions on the diagonal,
and the interactions between predators and prey on the off-diagonal. We
therefore write,

$\bm{B}=\left(\begin{matrix}2x^{*}(x^{*}+d)\bm{I}&-x^{*}y^{*}\bm{R}\\\
-x^{*}y^{*}\bm{R}^{T}&2y^{*}b\bm{I}\end{matrix}\right).$ (57)

Again, this allows us to write the dynamics in form of an Ornstein-Uhlenbeck
process as defined in Eq. (14).

#### Computing the Power Spectral Density

In the following, we use features of the bipartite interaction network. For
instance, all nodes that are connected to e.g. node $x_{i}$, will be prey
nodes $y_{j}$, and thus are not connected with each other (see Fig. 5). This
allows us to write the following recursion formulas for the mean power
spectral densities according to Eq. (41),

$\begin{split}\bm{\Psi}_{x}^{-1}&=\mathbb{E}[\bm{\chi}_{i}]-\mathbb{E}\left[\sum_{i\sim
j}^{N_{y}}\bm{\chi}_{ij}\bm{\Psi}_{y}\bm{\chi}_{ji}\right]\,,\\\
\bm{\Psi}_{y}^{-1}&=\mathbb{E}[\bm{\chi}_{j}]-\mathbb{E}\left[\sum_{j\sim
i}^{N_{x}}\bm{\chi}_{ji}\bm{\Psi}_{x}\bm{\chi}_{ij}\right]\,.\end{split}$ (58)

Recall that the top left entries of $\bm{\Psi}_{x}$ and $\bm{\Psi}_{y}$
deliver the mean power spectral densities for predators $\phi_{x}$ and prey
$\phi_{y}$ respectively. For the bipartite model, the helping matrices
$\bm{\chi}_{i},\bm{\chi}_{ij}$ (as defined in Eq. (22)) are given by,

$\begin{split}\bm{\chi}_{x}&=\left(\begin{matrix}0&-x+i\omega\\\
x+i\omega&2x(x+d)\end{matrix}\right),\\\
\bm{\chi}_{y}&=\left(\begin{matrix}0&-y+i\omega\\\
y+i\omega&2yb\end{matrix}\right),\\\
\bm{\chi}_{xy}&=\left(\begin{matrix}0&-\alpha y\\\ -\alpha x&-\alpha
xy\end{matrix}\right),\quad\bm{\chi}_{yx}=\left(\begin{matrix}0&\alpha x\\\
\alpha y&-\alpha xy\end{matrix}\right).\end{split}$ (59)

Inserting and writing out Eq. (58) gives,

$\begin{split}\left(\begin{matrix}\phi_{x}&r_{x}\\\
-\bar{r}_{x}&0\end{matrix}\right)^{-1}=&\left(\begin{matrix}0&-x+i\omega\\\
x+i\omega&2x(x+d)\end{matrix}\right)\\\
&+\alpha^{2}c_{x}\left(\begin{matrix}0&-\bar{r}_{y}xy\\\
r_{y}xy&\phi_{y}x^{2}-(r_{y}+\bar{r}_{y})x^{2}y\end{matrix}\right)\,,\\\
\left(\begin{matrix}\phi_{y}&r_{y}\\\
-\bar{r}_{y}&0\end{matrix}\right)^{-1}=&\left(\begin{matrix}0&-y+i\omega\\\
y+i\omega&2yb\end{matrix}\right)\\\
&+\alpha^{2}c_{y}\left(\begin{matrix}0&-\bar{r}_{x}xy\\\
r_{x}xy&\phi_{x}y^{2}+(r_{x}+\bar{r}_{x})xy^{2}\end{matrix}\right)\,,\end{split}$
(60)

where $c_{x},c_{y}$ are the number of connections per predator and prey
species respectively. Analogous to Eq. (37) we now derive a system of
equations and solve for $r_{x},r_{y}$ and $\phi_{x},\phi_{y}$. In the main
text we describe the features of the power spectral density deduced from this
system of equations.

###  IV.4 Simulations

#### Figure 1, left panels

A large random Lotka-Volterra ecosystem of the type described above was
generated. Parameters used were:
$N=1000,c=50,\gamma=-1,\sigma^{2}=1/4/c,b=0.05$. The solid line is the result
of Eq. (6), noting that $\mu=\sqrt{2\sigma^{2}/\pi}$. For the empirical power
spectrum, we used an Euler-Maryama time-stepping method to simulate a time
series of length $t_{\max}=2^{10}$ and time step $h=2^{-7}$. The power
spectrum for each species was calculated with a Fast Fourier Transform, and
the result averaged over all species. The top panel shows part of the time
series generated for the first species.

#### Figure 1, right panel

A two-trophic level model ecosystem was generated as described above.
Parameters in this case were:
$Nx=200,N_{y}=800,c_{x}=20,c_{y}=5,\alpha=10,b=1,d=1$. Time series and power
spectra were computed similarly to the left panels.

#### Figure 2

For the left panel, we generated Lotka-Volterra ecosystems with parameters
$N=1000,c=50,\sigma^{2}=0.5$, using $b=2+(1+\gamma)\sqrt{(}c*\sigma^{2})$ for
the simulations with $\gamma=0$ and $\gamma=1$. Time series and spectra were
computed similarly to Figure 1. For the right panels more care is needed.
Finite random matrices typically have a small number of eigenvalues that are
order $1/\sqrt{N}$ larger than predicted by the stability boundary in the
limit $N\to\infty$. To achieve the near-instability results in this figure, we
first generated the off-diagonal entries of the community matrices, then chose
the birth rate $b$ to put the rightmost eigenvalue of $\bm{A}$ exactly at
zero.

#### Figure 3

Parameters here are: $N=500,c=20,\gamma=-1,\sigma^{2}=1/4/c,b=0.2$. For the
‘direct’ results we numerically computed the power spectral density according
to the matrix formula (16). This was preferable to simulations of the time
series, as an extremely long time horizon is required to achieve good
resolution of the individual contributions to the power spectral density.

#### Figure 4

Parameters here are: $N_{x}=100,N_{y}=200,c_{x}=20,c_{y}=10,\alpha=5,b=1,d=1$.

## References

  * (1) Robert M. May.  Will a Large Complex System be Stable?  Nature, 238(5364):413, August 1972.
  * (2) Mark R. Gardner and W. Ross Ashby.  Connectance of Large Dynamic (Cybernetic) Systems: Critical Values for Stability.  Nature, 228(5273):784–784, November 1970.
  * (3) Jean Ginibre.  Statistical ensembles of complex, quaternion, and real matrices.  Journal of Mathematical Physics, 6(3):440–449, 1965.
  * (4) Alex James, Michael J. Plank, Axel G. Rossberg, Jonathan Beecham, Mark Emmerson, and Jonathan W. Pitchford.  Constructing Random Matrices to Represent Real Ecosystems.  The American Naturalist, 185(5):680–692, May 2015.  Publisher: The University of Chicago Press.
  * (5) Claire Jacquet, Charlotte Moritz, Lyne Morissette, Pierre Legagneux, François Massol, Philippe Archambault, and Dominique Gravel.  No complexity–stability relationship in empirical ecosystems.  Nature communications, 7:12573, 2016.
  * (6) Egbert G. Leigh Jr, James Rosindell, and Rampal S. Etienne.  Unified neutral theory of biodiversity and biogeography.  Scholarpedia, 5(11):8822, November 2010.
  * (7) Stefano Allesina and Si Tang.  Stability criteria for complex ecosystems.  Nature, 483(7388):205–208, March 2012.  Number: 7388 Publisher: Nature Publishing Group.
  * (8) Stefano Allesina and Si Tang.  The stability–complexity relationship at age 40: a random matrix perspective.  Population Ecology, 57(1):63–75, January 2015.  Publisher: John Wiley & Sons, Ltd.
  * (9) Stefano Allesina, Jacopo Grilli, György Barabás, Si Tang, Johnatan Aljadeff, and Amos Maritan.  Predicting the stability of large structured food webs.  Nature communications, 6(1):1–6, 2015.
  * (10) Jacopo Grilli, Tim Rogers, and Stefano Allesina.  Modularity and stability in ecological communities.  Nature communications, 7(1):1–10, 2016.
  * (11) Ian Donohue, Helmut Hillebrand, José M. Montoya, Owen L. Petchey, Stuart L. Pimm, Mike S. Fowler, Kevin Healy, Andrew L. Jackson, Miguel Lurgi, Deirdre McClean, Nessa E. O’Connor, Eoin J. O’Gorman, and Qiang Yang.  Navigating the complexity of ecological stability.  Ecology Letters, 19(9):1172–1185, 2016.
  * (12) V. Grimm and Christian Wissel.  Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion.  Oecologia, 109(3):323–334, February 1997.
  * (13) Richard Levins.  Coexistence in a Variable Environment.  The American Naturalist, 114(6):765–783, December 1979.  Publisher: The University of Chicago Press.
  * (14) A. R. Ives, K. Gross, and J. L. Klug.  Stability and Variability in Competitive Communities.  Science, 286(5439):542–544, October 1999.
  * (15) Clarence L. Lehman and David Tilman.  Biodiversity, Stability, and Productivity in Competitive Communities.  The American Naturalist, 156(5):534–552, November 2000.
  * (16) David Tilman, Peter B. Reich, and Johannes M. H. Knops.  Biodiversity and ecosystem stability in a decade-long grassland experiment.  Nature, 441(7093):629–632, June 2006.
  * (17) Lin Jiang and Zhichao Pu.  Different Effects of Species Diversity on Temporal Stability in Single‐Trophic and Multitrophic Communities.  The American Naturalist, 174(5):651–659, November 2009.
  * (18) Michel Loreau and Claire de Mazancourt.  Species Synchrony and Its Drivers: Neutral and Nonneutral Community Dynamics in Fluctuating Environments.  The American Naturalist, 172(2):E48–E66, August 2008.
  * (19) Veronik Campbell, Grace Murphy, and Tamara N. Romanuk.  Experimental design and the outcome and interpretation of diversity-stability relations.  Oikos, 120(3):399–408, March 2011.
  * (20) Ian Donohue, Owen L. Petchey, José M. Montoya, Andrew L. Jackson, Luke McNally, Mafalda Viana, Kevin Healy, Miguel Lurgi, Nessa E. O’Connor, and Mark C. Emmerson.  On the dimensionality of ecological stability.  Ecology Letters, 16(4):421–429, 2013.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/ele.12086.
  * (21) J-F. Arnoldi, M. Loreau, and B. Haegeman.  Resilience, reactivity and variability: A mathematical comparison of ecological stability measures.  Journal of Theoretical Biology, 389:47–59, January 2016.
  * (22) Kurt Wiesenfeld.  Noisy precursors of nonlinear instabilities.  Journal of Statistical Physics, 38(5-6):1071–1097, 1985.
  * (23) Rep Kubo.  The fluctuation-dissipation theorem.  Reports on progress in physics, 29(1):255, 1996.
  * (24) David Alonso, Alan J McKane, and Mercedes Pascual.  Stochastic amplification in epidemics.  Journal of the Royal Society Interface, 4(14):575–582, 2007.
  * (25) Tobias Galla.  Intrinsic noise in game dynamical learning.  Physical review letters, 103(19):198702, 2009.
  * (26) A. J. McKane and T. J. Newman.  Predator-Prey Cycles from Resonant Amplification of Demographic Stochasticity.  Physical Review Letters, 94(21):218102, June 2005.
  * (27) Eugene P Wigner.  On the distribution of the roots of certain symmetric matrices.  Annals of Mathematics, pages 325–327, 1958.
  * (28) Maria Dornelas et al.  Biotime: A database of biodiversity time series for the anthropocene.  Global Ecology and Biogeography, 27(7):760–786, 2018.
  * (29) Terence Tao and Van Vu.  Random matrices: the circular law.  Communications in Contemporary Mathematics, 10(02):261–307, 2008\.
  * (30) Terence Tao, Van Vu, Manjunath Krishnapur, et al.  Random matrices: Universality of esds and the circular law.  The Annals of Probability, 38(5):2023–2065, 2010.
  * (31) Lewi Stone.  The feasibility and stability of large complex biological networks: a random matrix approach.  Scientific Reports, 8(1):8246, May 2018.  Number: 1 Publisher: Nature Publishing Group.
  * (32) Theo Gibbs, Jacopo Grilli, Tim Rogers, and Stefano Allesina.  Effect of population abundances on the stability of large random ecosystems.  Physical Review E, 98(2):022410, August 2018.
  * (33) Crispin Gardiner.  Stochastic Methods: A Handbook for the Natural and Social Sciences.  Springer Series in Synergetics. Springer-Verlag, Berlin Heidelberg, 4 edition, 2009.
  * (34) Tobias Galla.  Dynamically evolved community size and stability of random lotka-volterra ecosystems (a).  EPL (Europhysics Letters), 123(4):48004, 2018.
  * (35) Theo Gibbs, Jacopo Grilli, Tim Rogers, and Stefano Allesina.  Effect of population abundances on the stability of large random ecosystems.  Physical Review E, 98(2):022410, 2018.
  * (36) Alan J McKane, Tommaso Biancalani, and Tim Rogers.  Stochastic pattern formation and spontaneous polarisation: the linear noise approximation and beyond.  Bulletin of mathematical biology, 76(4):895–921, 2014.
  * (37) Jeffrey Pennington and Pratik Worah.  Nonlinear random matrix theory for deep learning.  In Advances in Neural Information Processing Systems, pages 2637–2646, 2017.
  * (38) José Moran and Jean-Philippe Bouchaud.  Will a large economy be stable.  Available at SSRN, 2019.
  * (39) Feng Luo, Yunfeng Yang, Jianxin Zhong, Haichun Gao, Latifur Khan, Dorothea K Thompson, and Jizhong Zhou.  Constructing gene co-expression networks and predicting functions of unknown genes by random matrix theory.  BMC bioinformatics, 8(1):299, 2007.
  * (40) Assaf Almog, M Renate Buijink, Ori Roethler, Stephan Michel, Johanna H Meijer, Jos HT Rohling, and Diego Garlaschelli.  Uncovering functional signature in neural systems via random matrix theory.  PLoS computational biology, 15(5):e1006934, 2019.
  * (41) R. L. Stratonovich.  On a Method of Calculating Quantum Distribution Functions.  Soviet Physics Doklady, 2:416, July 1957.
  * (42) J. Hubbard.  Calculation of partition functions.  Phys. Rev. Lett., 3:77–78, Jul 1959.
  * (43) M Mézard, G Parisi, and M. A Virasoro.  SK Model: The Replica Solution without Replicas.  Europhysics Letters (EPL), 1(2):77–82, January 1986.
  * (44) Marc Mezard, Giorgio Parisi, and Miguel Angel Virasoro.  Spin glass theory and beyond.  Number v. 9 in World Scientific lecture notes in physics. World Scientific, Singapore ; New Jersey, 1987.  OCLC: ocm14929802.
  * (45) Tim Rogers, Isaac Pérez Castillo, Reimer Kühn, and Koujin Takeda.  Cavity approach to the spectral density of sparse symmetric random matrices.  Physical Review E, 78(3):031116, September 2008.
  * (46) Tim Rogers and Isaac Pérez Castillo.  Cavity approach to the spectral density of non-Hermitian sparse matrices.  Physical Review E, 79(1):012101, January 2009.
  * (47) Fernando Lucas Metz, Izaak Neri, and Tim Rogers.  Spectral theory of sparse non-hermitian random matrices.  Journal of Physics A: Mathematical and Theoretical, 52(43):434003, 2019.