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

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 ({mathbb{E}}[{boldsymbol{zeta }}(t){boldsymbol{zeta }}{(t^{prime} )}^{T}]=delta (t-t^{prime} ){boldsymbol{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 ({mathbb{E}}[{boldsymbol{xi }}(t){boldsymbol{xi }}{(t+tau )}^{T}]) at equilibrium,

From ref. ³³ 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 ({boldsymbol{{mathcal{R}}}}={({rm{i}}omega -{boldsymbol{A}})}^{-1}=langle {boldsymbol{u}}{{boldsymbol{v}}}^{dagger }rangle). Thus we may write,

where ({mathcal{Z}}=| {boldsymbol{A}}-iomega {boldsymbol{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^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 sites i, ({{mathcal{H}}}_{i}), and contributions from interactions between i and j, ({{mathcal{H}}}_{ij}),

These terms can be decomposed as ({{mathcal{H}}}_{i}={{boldsymbol{w}}}_{i}^{dagger }{{boldsymbol{chi }}}_{i}{{boldsymbol{w}}}_{i}) and ({{mathcal{H}}}_{ij}={{boldsymbol{w}}}_{i}^{dagger }{{boldsymbol{chi }}}_{ij}{{boldsymbol{w}}}_{j}), where we introduce the compound spins ({{boldsymbol{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 ({{mathbf{Psi }}}_{i}=langle {{boldsymbol{w}}}_{i}{{boldsymbol{w}}}_{i}^{dagger }rangle). 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}_{j}^{(i)}) 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’ ({{mathbf{Psi }}}_{jl}^{(i)}) 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 ({{mathbf{Psi }}}_{jl}^{(i)}),

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 ({{mathbf{Psi }}}_{jl}^{(i)}) only gives non-zero 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.

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 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 ({boldsymbol{{mathcal{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.

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 for ({phi }_{i}^{{rm{SDA}}},{r}_{i}^{{rm{SDA}}}), which delivers

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 ({{mathbf{Psi }}}_{m}^{{rm{MF}}}), 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 straightforward, 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 intra-group interactions simplifies the analysis.

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 ≫ 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 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 ∈ [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 γ = 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 (gamma ={mathbb{E}}[{alpha }_{ij}{alpha }_{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 ({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 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 (gamma ={mathbb{E}}[{A}_{ij}{A}_{ji}]). For a random community matrix (i.e. γ = 0), we recover the stability criterion that has been proven by May¹, (sqrt{c{sigma }^{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 (langle {eta }_{i}(t),{eta }_{j}(t^{prime} )rangle =delta (t-t^{prime} ){B}_{ij}). 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 off-diagonal 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, ({boldsymbol{xi }}=sqrt{V}({boldsymbol{x}}-{{boldsymbol{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 (over non-zero entries) reaction rate (mu ={mathbb{E}}[{R}_{ij}]=sigma sqrt{2/pi }).

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 ({mathbb{E}}[{A}_{ii}]=-b), and the noise matrix ({mathbb{E}}[{B}_{ii}]=2b+cmu). 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. (38) we obtain,

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

Trophic structure model

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_xc_x = N_yc_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 corresponding stochastic differential equations,

where η_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,

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

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,

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 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’²⁰. 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²². 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⁵³;

The stationary distribution of ξ is then ({P}_{{rm{st}}}({boldsymbol{xi }})={mathcal{N}}({boldsymbol{0}},{mathbf{Sigma }})). 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.3 in the power spectrum in the upper right hand panel indicates that the trajectories in the upper left panel exhibit quasi-cycles (i.e. have a typical frequency, see inset), while the peak at ω = 0 at following decay of the power spectrum in the middle right hand panel indicates that the trajectories in the middle left panel do not exhibit 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).

Figure details

Figure 1, left panels

A large random Lotka-Volterra ecosystem of the type described above was generated. Parameters used were: N = 1000, c = 50, γ = −1, σ² = 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⁻⁷. 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, α = 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, σ² = 0.5, using (b=2+(1+gamma )sqrt{left(right.}c* {sigma }^{2}left.right)) for the simulations with γ = 0 and γ = 1. Time series and spectra were computed similarly to Fig. 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 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.

Figure 3

Parameters here are: N = 500, c = 20, γ = −1, σ² = 1/4/c, b = 0.2. For the ‘direct’ results we numerically computed the power spectral density according to the matrix formula in Eq. (17). This was preferable to simulations of the time series, as a 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, α = 5, b = 1, d = 1.

Figure 5

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 Eq. (38), with parameters b, cσ² (a composite parameter), γ, and an additional scale parameter for overall noise strength.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com