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

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. We show that different network structures leave distinct signatures in the spectrum of fluctuations, and demonstrate the application of our theory to the analysis of ecological time-series data of plankton abundances.


INTRODUCTION
"Will a large complex system be stable?" asks the title of Robert May's seminal 1972 paper [1] 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 [2], with results of Ginibre [3] 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 [4,5] -with some arguing that these serve best as a null model for ecosystem structure [6]-this growing field has continued to provide insights into the mechanisms that promote ecosystem stability. For instance, Allesina and Tang [7,8] 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 [9] or community [10] 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 [11]. In order to make clear the differences between these interpretations, Grimm and Wissel [12] 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 [13][14][15][16][17][18][19][20]. In [21], Suweis et al. propose to study the attenuation of perturbations as they propagate through ecological networks, introducing measures of reactivity and localization. Taking a different approach, recently Arnoldi et al. [22] 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 can drive a system away from equilibrium, and thus are important precursors to linear and nonlinear instabilities [23].
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). The relationship between such power spectra and temporal notions of ecological stability is multifaceted. Particular points of interest are the height of spectrum, which gives information about the magnitude of stochastic fluctuations, the locations of non-zero peaks corresponding to quasi-cyclic signals, or a peak at zero indicating baseline wander. Moreover, the Fourier transform of the spectrum yields the autocorrelation structure of the stochastic trajectories. We provide a brief guide to these concepts in the Methods section Interpreting the Power Spectral Density in the Context of Temporal Stability, though in the main text we will refrain from ascribing overly simplified interpretations to power spectra. Beyond providing a more detailed view of temporal stability, an investigation of power spectra yields a number of further advantages. For instance, power spectra are readily computed from empirical data and provide detailed information about intrinsic fluctuations and (via the fluctuation dissipation theorem [24]) 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 [25][26][27]. 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, while the applicability of the theory to larger systems is limited by comparatively slow numerical schemes, and difficulty parameterising large models. Here, by applying techniques from the statistical physics of complex systems, we demonstrate the possibility of deriving exact analytic formulae for the power spectra of large random and noisy dynamical systems.
We apply our method to characterise the stochas-tic 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 predatorprey 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 [28] and points to the wide applicability of the theory we develop. Importantly, the universal character of the power spectrum we derive is independent of the choice of random variables in the model, and only depends on the aggregate properties we identify. 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.
The paper is structured as follows. First we demonstrate how to compute the mean power spectral density of a large random Lotka-Volterra system in section Interaction types determine fluctuation spectra, showing how different dominant interaction types result in distinct fluctuation power spectra. We then show how to compute the spectrum for an individual species within the large random ecosystem system in section Species fluctuations exhibit strong heterogeneity, and further generalise the method in section Trophic structure induces fluctuation frequency gap to consider bipartite interaction networks, showing how a two-level trophic system can leave a distinct fingerprint in the power spectrum of an ecosystem. Readers interested in the potential of our results as a tool for analysis of real time series data may wish to jump to section Confronting RMT in theoretical ecology with time series data which provides a proof of concept in this direction. Here we explore an ecological time series dataset of plankton abundances, showing how our results provide a technique to infer the structural details of real ecosystems. Full derivations of our analytic results are provided in the Methods section, along with detailed descriptions of the models we use for demonstration throughout this paper.

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) : Here, the coefficients α ij for i = j describe the interaction between species i and j, and the η i (t) are Gaussian noise term with correlations . 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 ≡ b and density dependent mortality rate α ii ≡ −b. The other interaction coefficients α ij are chosen at random so that (i) each species interacts with an average of c others (for each possible interaction we include it with probability c/N , independent of all others), (ii) interactions have mean strength E|α ij | = µ and second moment Eα 2 ij = σ 2 , (iii) the correlation is controlled by the symmetry parameter 1]. Crucially, the full details of the distribution of the parameters α ij are not required, thanks to the universality of property of large random matrices [29,30].
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 γ. At γ = −1 all interactions are of the predator-prey type (α ij = −α ji ), at γ = 1 only purely mutualistic (α ij = α ji > 0) or competitive (α ij = α 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 [8], is stable provided b > √ cσ 2 (1 + γ) (we refer to [31,32] for stability and feasibility of equilibrium states with heterogeneous species abundance distributions). Around this fixed point, species in the stochastic system in Eq. (1) will exhibit approximately linear fluctuations ξ i , described by an Ornstein-Uhlenbeck process of the form Here A is the Jacobian of Eq. (1), known as the community matrix in the context of theoretical ecology, and ζ is an N -vector of Gaussian white noise with correlation matrix B = B(x * ). We assume that the equilibrium point at x = x * is linearly asymptotically stable (i.e. stable in the mathematical sense described by May [1]) and now proceed to investigate its temporal stability as characterised by the power spectra (see Methods section Interpreting the Power Spectral Density in the Context of Temporal Stability). The power spectral density of fluctuations Φ(ω) is defined as the Fourier transform of the covariance E[ξ(t)ξ(t + τ ) T ]. For multivariate Ornstein-Uhlenbeck processes one can show (see e.g. [33]) that 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 A and B. We derive an expression for the mean-field power spectral density φ(ω) = E[Φ ii ] in terms of the resolvent function r. Specifically, where expectation is taken only over the non-zero entries of A and B, and r ∈ C solves the self-consistent equation This result holds for general random matrix models in which interaction parameters are drawn from the same distribution for all species pairs (and we later show how the method can be extend for other model types with species-specific parameters using a single-defect approximation or partitioned networks). For the present case of our Lotka-Volterra model, the community matrix coincides with the interaction matrix (that is A ij = α ij ). In the methods we derive the rules E[B ii ] = 2b + cµ, E[A ij B ij ] = 0 for the statistics of the noise correlation matrix. To get a sense for the information contained in Eq. (4), we explore the result for several cases with varying interaction structures.
First consider a weak interaction limit where the difference between species is rather small, so that σ 2 1. In this case we find a simple Lorentzian spectral density: Fluctuations of this type are indicative of a highly stable system in which the balance of interaction types γ 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 γ = −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 ω 2 < 4cσ 2 then with an order 1/ω tail outside this range. Note that Eq. (7) 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 [28]. The result is illustrated in Fig. 1, left panel.
For a random mixture of interaction types with γ = 0 no approximations are necessary as Eq. (4) simplifies to The stability boundary here is given by b 2 = cσ 2 . The above result therefore implies the emergence of a 1/ω 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. γ = +1), the full solution to Eq. (7) in this case is complicated, but for stable systems appears qualitatively similar to the result Eq. (8) above. Near the stability boundary, however, we find another behaviour.
(9) In contrast to the previous case, this power spectrum exhibits a pole of order 1/ √ ω at low-frequency, followed by a 1/ω 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).

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 Eq. (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 Eq. (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 φ SDA i to the spectral density of fluctuations for species i, given by the expression where φ MF and r MF are the mean-field power spectrum and resolvent obeying Eq. (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. We immediately 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 quasiperiodic fluctuations in these populations 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.

Trophic structure induces fluctuation frequency 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. Here we focus on a bipartite predator-prey network as an example. 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 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 (x * , y * ), around which linear-order fluctuations will occur, analogously to Eq. (2) above. We compute a community matrix of the form where the first i = 1, . . . , N x rows and columns represent the predator species, and the remaining j = N x + 1, . . . , N x + N y rows and columns correspond to the prey species. The noise matrix is derived from the underlying individual-based model (see Methods) and given by 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 (61) to be solved for the mean contributions to the resolvent r x , r y , and to the power spectrum, φ x , φ 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 gap in excited frequencies around zero. Examination of the system in Eq. (61) allows us to determine the window of excited frequencies to be bounded by the critical frequencies 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 Eq. (14). 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 (see section  17), and the corresponding mean power spectral density (solid line) are obtained by solving Eq. (61) (inset shows the prey species). Excited frequencies are confined to a band (shaded) between critical frequences given in Eq. (14). The order 1/N peak at ω = 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).
Interpreting the Power Spectral Density in the Context of Temporal Stability).
Comparisons between simulations and our analytical results shows another interesting feature: an order 1/N disagreement at frequency ω = 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 [27]. 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.

Confronting RMT in theoretical ecology with time series data
Although hugely influential in the field of theoretical ecology over the last 50 years, traditional work on RMT has so far led to rather limited empirically testable insights. The central issue is that while many ecological considerations can be incorporated in a random matrix model, each leads to a binary outcome; the system is either stable or unstable to small perturbations. Thus testing the predictions of these models demands the time-intensive task of measuring real species interaction networks (which are assumed to be stable) and asking whether they indeed tend to be weakly connected (as suggested by May [1]), have a dominance of predator-prey interactions (as suggested by Allesina and Tang [7]), or satisfy some other prediction of the theory. In contrast, the approach presented in this paper offers the tantalizing prospect of directly linking the ecological RMT framework with comparatively easy-to-obtain time series data.
To trial the use of our methods in the analysis of real ecological data, we have investigated a high-resolution time series dataset for the abundance of coastal plankton species, taken over a period of 88 consecutive days [34]. In Figure 5 we show the estimated empirical mean power spectrum from the data (circles), compared to that of the best fit Lotka-Volterra random ecosystem model according to our theory. Full details of the data analysis and fitting are given in the Methods. Examination of this fit reveals several qualitative features of the implied ecological interactions.
First, we note that the best fit value for the interaction symmetry parameter is γ = 0.81, implying an ecosystem in which predator-prey interactions are scarce, and is more likely dominated by competition. Trust in this finding is strengthened by the fact that the fit is quite sensitive to this parameter; the dashed line in Figure 5 gives the best fit under the constraint γ < 0, which performs poorly, especially for low frequency.
Interestingly, when viewed in logarithmic axes ( Fig. 5  main panel), the plankton abundance power spectrum appears to exhibit a similar change of scaling between high and low frequency ranges to that seen in Fig. 2 for the case of symmetric interactions near instability. We can assess the closeness to instability by considering the spectrum inferred from the best fit model, as shown in the lower right panel of Figure 5. The rightmost edge is very close to zero, implying ecosystem dynamics which are close to instability. This feature corresponds to the large peak at zero in the power spectral density, which suggests that low frequency perturbations to the overall species abundances are very slow to relax.
One feature of the spectrum not reproduced by the simple models considered thus far is the smaller additional peak around ω ≈ 1.5. This peak has a few possible explanations: a external effect of some sort; possible secondary structure in the ecological interaction network, which could manifest on a system wide scale such as the trophic structure analysed in the previous section; or a feature isolated to a smaller number of more dominant species. A further limitation of the model used here is the assumption of uniform species abundance; in reality, species abundances tend to be distributed log-normally, with few species contributing to the majority of ecosystem biomass. Incorporating such model refinements are well-within the bounds analytical tractability for our approach (see, for instance, [32] and [31]); we hope and expect the theoretical groundwork we have developed here will pave the way for the investigation of such features in future studies. We plot the power spectral density of data taken from [34] (blue circles) alongside that of a fitted Lotka-Volterra random ecosystem model (solid line). Fitted parameters are b = 0.6643, cσ 2 = 0.1316, γ = 0.8078, data averaged over n = 3 samples per day. Also shown is the best fit under the restriction γ < 0 (dashed line). The lower right panel shows the spectral boundary inferred from the fit (black ellipse), along with the eigenvalues of a sample random community matrix for illustration.
In the above, we have illustrated the use of our methods to infer details of the structure and stability of real ecosystems from time series data, as well as to identify departures from the unstructured assumptions of standard RMT models. Indeed, such departures are present in many real world ecosystems, with important consequences for the validity of any predictions made within the standard RMT framework [4]. In contrast, RMT has recently found renewed attention in the field of microbiome research, where it is believed that the key conceit of standard RMT models (that communities are unstructured) holds [35]. However in this field, the spectre of model parameterization again raises its head.
In [35], a species-interaction network presented in [36] was used to parameterize an RMT model and show that ecological interactions in the microbiome tended to be weak and non-cooperative. The species interaction network determined in [36] was itself the result of fitting mouse intestinal microbiome abundances to a deterministic generalized Lotka-Volterra model. However, fully fitting this model required disturbing the mouse microbiota away from its equilibrium state using antibiotics (S fixed point species abundances are insufficient to parameterize an S × S species interaction network, so data on non-equilibrium transient trajectories were required). While such experimental manipulation may be permissible for studying the microbiota of model organisms such as mice, the ethical issues of such experimentation in humans has raised questions about the informativeness of temporal data for understanding microbial communities such as the human gut microbiota [37].
In contrast to the approach taken in [35], our methodology requires no external perturbation to a host's microbiome, relying as it does solely on the natural demographic fluctuations present in any finite population. In addition, our approach allows the RMT model itself to be directly parameterized through data, rather than requiring the fitting of an intermediate model.

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, demonstrating the flexibility of the method and usage across a broader range of models. Finally we fitted our model to existing time series data sets, that suggest a majority of competitive or mutualistic interactions within plankton ecosystems. 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 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 γ, 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/ω 2 or 1/ √ ω 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 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 [38], heterogeneous turnover rates [39], or explicitly spatial models where spatio-temporal patterning may persist [40].
From an ecological perspective, it is 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. Beyond simply providing more detail as to the temporal dynamics of an ecosystem around an equlibrium point, our method has also opened up the exciting possibility of identifying the signature of a certain interaction structures in the power spectrum of oscillations in data gathered in the field. We fitted our model to a highly resolved time series data set on a plankton ecosystem. We found that the empirical data is indicative of an ecology dominated by competitive and mutalistic interactions, with far fewer predator-prey interactions. This insight in consistent with recent results that suggest that self-regulation (competition) and facilitation (mutualism) are widespread in phytoplankton communities [41].
In order to further realise this vision, 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. Our theoretical approach encourages further work focused on the application to data sets gathered in field studies with modifications more suitable for the method we presented.
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 [42], finance [43], biochemistry [44] and neuroscience [45]. 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 insight into the model, which relies only on statistical meta parameters.

Power Spectral Density for a General Ornstein-Uhlenbeck Process
In the following we develop a method to compute the power spectral density of N -dimensional Ornstein-Uhlenbeck processes, where ζ(t) is an N -vector of Gaussian white noise with correlations E[ζ(t)ζ(t ) T ] = δ(t − t )B. The matrix A determines the mean behaviour of ξ and is considered to be locally stable, i.e. all eigenvalues of A have negative real part. Using the matrices A and 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 A and B to deduce a self-consistent scheme of equations whose solution contains information on the power spectral density.
We start with the definition of the power spectral density Φ(ω) as the Fourier transform of the covariance E[ξ(t)ξ(t + τ ) T ] at equilibrium, From [33] on multivariate Ornstein-Uhlenbeck processes, we know that the power spectral density can also be written in the form of the matrix equation, 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 Eq. (17) as a complex Gaussian integral reminiscent of problems appearing in the statistical physics of disordered systems. Our approach in the following is to treat ω as a fixed parameter and drop the explicit dependence from our notation. We begin by writing Simplification of the integrand is achieved by unpicking the matrix inversion in the exponent via a Hubbard-Stratonovich transformation [46,47]. 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' Introducing a bracket operator we can obtain succinct expressions for the power spectral density Φ = uu † as well as the resolvent matrix R = (iω − A) −1 = uv † . Thus we may write, where Z = |A − iωI| 2 /π 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 [48,49]. Further applications of the method include the analysis of the eigenvalue distribution in sparse matrices [50][51][52]. 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. (19) into the sum of its local contributions at site i, H i , and contributions from interactions between i and j, H ij , These terms can be decomposed as where we introduce the compound spins w i = (u i , v i ) T and transfer matrices, Let us focus on the power spectral density of a particular variable ξ i , obtained from the diagonal element φ i = Φ ii . For this we compute the single-site marginal f i by integrating over all other variables, Alternatively, φ i can be obtained as the top left entry of the covariance matrix Ψ i = w i w † i . We write the covariance matrix as the integral, which could also be expressed in terms of a Gaussian integral, By comparing Eqs. (25) and (26) we find that We now insert Eq. (22) into Eq. (24) and obtain, where we write f (i) j for the 'cavity marginals', In essence, the above discussion amounts to organising the 2N integrals in Eq. (21) 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. We will further use this notation for the 'cavity covariance matrix' Ψ (i) jl introduced in the following. Next we perform the integration in Eq. (28) and compare to the form in Eq. (27). We thus obtain a recursion formula for the covariance matrix Ψ i and the cavity covariance matrices Ψ where the notation i ∼ 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 Ψ (i) jl only gives nonzero contributions if j = l. We therefore reduce Eq. (30) and obtain for the covariance matrix, Similarly, the cavity covariance matrix obeys the equation, Here we use that Ψ (i,j) = Ψ (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 in Eq. (31) 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, as its contribution is negligible compared to the full system. Hence the system in Eq. (31) can be reduced to a smaller set of equations approximately satisfied by the matrices Ψ i : The power spectral density φ i can be obtained as the top left entry of Ψ 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.
a. 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 meanfield (MF) formulation we assume that the mean degree and excess degree are approximately equal, and replace all quantities in Eqs. (31) and (32) with their average. Ψ i = Ψ MF ∀i. We then obtain the following recursion equation, In order to solve this equation, we parameterise, where the top left entry φ corresponds to the mean power spectral density, and we introduce r as the mean diagonal element of the resolvent matrix R. Finally by inserting the ansatz of Eq. (35) into Eq. (34) we obtain, 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 ∼ j (i.e. non-zero matrix entries).
From Eq. (36) above, we obtain the equations, We solve the second equation in Eq. (37) for r and write the mean power spectral density in terms of r, This equation informs the first part of the results presented in the main text.
b. 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. (35) for a single individual. Moreover, we replace all other quantities with the respective mean-field approximation. Specifically, we obtain We solve this equation c. 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 the section Trophic Structure Model for a simple example).
In order to handle different connected groups we make use of the cavity method as in Eqs. (31) and (32). 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, Similar to the previous sections we replace all quantities with a mean-field average Ψ MF m , but for each group separately. Hence we obtain M equations of the form In order to compute the mean power spectral density for different groups separately, we use a parameterisation as in Eq. (35) for each group. Therefore we have, for all m = 1, . . . , M . This delivers 2M equations to solve for all r m and φ m . Numerically this is straight forward, although algebraically long-winded for the general case. However, the equations simplify for special cases.
In the section Trophic Structure Model we demonstrate this method for a bipartite network where a lack of intragroup interactions simplifies the analysis.
Large Lotka-Volterra Ecosystem d. Model Description First, we define the framework for a general Lotka-Volterra ecosystem with N species and a large but finite system size V 1. Note that this parameter can be interpreted as a scaling factor for the fluctuation amplitude and thus, larger systems exhibit higher stability and quantitative reliability for our analytic results. Let X i denote the number of individuals and x i = X i /V the density of species i = 1, . . . , N . We start from the following set of reactions that define the underlying stochastic dynamics of the system: 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 predatorprey 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 = N P c interaction partners. The reaction rates are considered to be i.i.d. random variables drawn from a half-normal distribution |N (0, σ 2 )|, where we write for the mean reaction rate µ = E[R ij ] = σ 2/π and raw second moment ij ]. For each interaction pair, the interaction type is chosen such that the proportion of predator-prey pairs is p ∈ [0, 1], and all non-predatorprey 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 γ = 1 − 2p, where γ = −1 if all interactions are of predator-prey type (p = 1), and similarly γ = +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 γ = 0. Later we will see that γ is equivalent to the correlation of signed interaction strengths.
In the limit V → ∞, the dynamics of the species density x i obey the ordinary differential equations, where α ij are the interaction coefficients with |α ij | = |α 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 α ij = α ji > 0, and α ij = α ji < 0 for competitive interactions. In the case of predator-prey interactions the coefficients have opposite sign α ij = −α ji . Hence the symmetry parameter as described above is given by the correlation of interaction coefficients γ = E[α ij α ji ]. Furthermore, in order to ensure bounded species densities, we require negative self-interactions α ii = −R ii < 0. If species live in isolation (i.e. when α ij = 0 ∀i = j), we see that the densities approach the 'effective' carrying capacity K i = −b i /α 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 α ii = −b. This gives the effective carrying capacity K = 1 for all species.
The fixed point x * at the deterministic equilibrium state is given by, We assume a random mixture of mutualistic and competitive interactions with equal proportions, and therefore the interaction coefficients α ij have zero mean (∀i = j). Furthermore, we postulate that for large ecosystems where N → ∞, the equilibrium state 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 x * , and our assumptions are for mathematical simplification only. The results are independent of the particular equilibrium configuration, as long as a stable equilibrium can be measured and extracted from data (we discuss a few caveats where we apply our method to time series data from a plankton ecosystem). This assumption allows us to write the Jacobian matrix for a linearisation around the equilibrium state, with elements, 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 = α ij . The local stability of such community matrix A is given by the elliptic law [7,8]. 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, with the average number of connections c per species, and the correlation γ = E[A ij A ji ]. For a random community matrix (i.e. γ = 0), we recover the stability criterion that has been proven by May [1], √ cσ 2 < b. If γ < 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 γ = −1 (i.e. A ij = −A ji ∀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. γ = +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, where η i (t) are Gaussian random variables with The noise matrix 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 offdiagonal elements depend on the type of interaction between species i and j. We assume that only predator-prey type interactions contribute to the covariance of species fluctuations (i.e. that only predator-prey interactions involve the simultaneous change in abundance of a species pair). Therefore, we write We next linearise around the fixed point to obtain a new equation for the fluctuations, ξ = √ V (x−x * ), which has the form of an Ornstein-Uhlenbeck process as defined in Eq. (15). Recall that in our simplified model the equilibrium abundance x * = 1 (note however, that in general the entries of the noise matrix B depend on the particular fixed point of a given system). Therefore we write for the noise matrix evaluated at the fixed point where µ is given as the mean reaction rate µ = E[R ij ] = σ 2/π. e. Computing the Power Spectral Density Let us now compute the mean power spectral density φ of the process described above using Eq. (38) as starting point. We replace the necessary quantities that we obtain from the community matrix A and noise matrix B as defined in the previous section. In particular, we have the expected diagonal elements of the community matrix E[A ii ] = −b, and the noise matrix E[B ii ] = 2b + cµ. Moreover the raw second moment of the non-zero interactions is given by E[A ij ] = σ 2 and the correlation E[A ij A ji ] = γσ 2 . We use that E[A ij B ij ] = 0 ∀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. Here we have nx = 4 predators with cx = 3 prey each, and ny = 6 prey with cy = 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.
Plugging in these quantities into Eq. (38) we obtain, In the main text we explore the theoretical ecological consequences of this result.
Trophic Structure Model f. 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. 6. 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 = α, 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 αc x and α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, . . . , N x , and y j the density of prey species j = 1, . . . , N y . In the deterministic limit where V → ∞ we then write the following ODEs, Given the fixed number of connections c x , c y and interaction strength α, we can simplify the ODEs to two equations for the average predator and prey densities, 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, Next, we write the Jacobian matrix for a linearisation around the non-trivial fixed point. The community matrix takes the form, where the first i = 1, . . . , N x rows and columns represent the predator species, and the remaining j = N x + 1, . . . , N x + N y rows and columns correspond to the prey species. For a large but finite system size V we write the cor-responding stochastic differential equations, where η i,j (t) are Gaussian noise with η i (t), η j (t ) = δ ij (t − t )B ij . The noise matrix is given by the selfand total interactions on the diagonal, and the interactions between predators and prey on the off-diagonal. We therefore write, Again, this allows us to write the dynamics in form of an Ornstein-Uhlenbeck process as defined in Eq. (15). g. 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. 6). This allows us to write the following recursion formulas for the mean power spectral densities according to Eq. (42), Recall that the top left entries of Ψ x and Ψ y deliver the mean power spectral densities for predators φ x and prey φ y respectively. For the bipartite model, the helping matrices χ i , χ ij (as defined in Eq. (23)) are given by, , Inserting and writing out Eq. (59) gives, −r y xy r y xy φ y x 2 − (r y +r y )x 2 y , where c x , c y are the number of connections per predator and prey species respectively. Analogous to Eq. (38) we now derive a system of equations and solve for r x , r y and φ x , φ y . In the main text we describe the features of the power spectral density deduced from this system of equations.
Interpreting the Power Spectral Density in the Context of Temporal Stability For orientation, we here provide some interpretation of the power spectral density in the context of temporal stability. Essentially when we talk about temporal stability, we can can be referring to one of two measures. The first is how far stochastic trajectories tend to stray from their equilibrium value over long time horizons. We refer to this as 'variability' [20]. The second is how quickly population abundances tend to change over finite time horizons. We will characterise this by the 'temporal autocorrelation'.
The variability can be characterised by the variance in time-averaged trajectories around the mean [22]. For a system such as Eq. (2), which we recall can be a linear approximation for a nonlinear system such as Eq. (1), we find that ξ is normally distributed with zero mean and a covariance matrix, Σ, that solves the following Lyapunov equation [53]; The stationary distribution of ξ is then P st (ξ) = N (0, Σ). For instance, in the left panels of Fig. 7, we show stochastic trajectories for two different systems, with standard deviations marked by black dashed lines. Meanwhile the marginal normal distribution for these trajectories is plotted in the inset of the right panels of Fig. 7. A system can then be said to be 'less stable' (in a temporal sense) if it has a greater variability. A consideration of the solutions to Eq. (62) shows that this measure of temporal stability is highly correlated with asymptotic stability; less stable deterministic systems tend to have stochastic counterparts with higher variance around equilibrium states. Despite the fact that the trajectories in Fig. 7 have the same variance (see black dashed lines in left panels and inset plots in right panels) it is clear that they have very different temporal structure. While these differences are entirely masked by the measure of variability (which time-averages out the temporal structure), such differences are captured by the power-spectral density (see Fig. 7, right panels). For instance, the peak at ω ≈ 0. hibit quasi-cycles (i.e. do not have a typical frequency, see inset).
In the context of temporal stability, the relationship between the power spectra and the autocorrelation function ξ(t)ξ(t − τ ) is of particular importance. By the Wiener-Khinchin theorem, we know that the autocorrelation function is given by the Fourier transform of the power spectrum. This is shown in the bottom panel of Fig. 7. The autocorrelation of the trajectory in the upper panels decays rapidly with time. In contrast, the autocorrelation of the trajectory in the middle panel decays more slowly. This can be clearly seen in the inset trajectory plots (left hand panels, top and middle). Thus we see that a distinct measure of temporal stability exists that is more appropriate over shorter time horizons; a system can be said to be 'less stable' over finite times if it has a more rapidly decaying autocorrelation function. This measure of temporal stability is more weakly correlated with asymptotic stability than its counterpart, variability, as it is affected by the magnitude imaginary parts of the system's eigenvalues (rather than their real parts, as in asymptotic stability).   6), noting that µ = 2σ 2 /π. 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.  For the left panel, we generated Lotka-Volterra ecosystems with parameters N = 1000, c = 50, σ 2 = 0.5, using b = 2 + (1 + γ) (c * σ 2 ) for the simulations with γ = 0 and γ = 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/ √ N larger than predicted by the stability boundary in the limit N → ∞. 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 A exactly at zero.   The dataset 41467 2017 2571 MOESM6 ESM.xlsx was imported into Matlab and processed as follows: We took the average of the three reported daily measurements to construct an 88-day time series for each species. To limit boundary effects we discarded all species with at least one with zero measured abundance, in doing so retaining 100 species. The mean was subtracted and then the power spectrum fitted using the covariance method with 8th order autoregression. The model fitting was achieved with a non-linear least squares method applied to our equation (38), with parameters b, cσ 2 (a composite parameter), γ, and an additional scale parameter for overall noise strength.

Data Availability
Plankton abundance data used in Fig 5 are