Dispersal-induced instability in complex ecosystems

In his seminal work in the 1970s, Robert May suggested that there is an upper limit to the number of species that can be sustained in stable equilibrium by an ecosystem. This deduction was at odds with both intuition and the observed complexity of many natural ecosystems. The so-called stability-diversity debate ensued, and the discussion about the factors contributing to ecosystem stability or instability continues to this day. We show in this work that dispersal can be a destabilising influence. To do this, we combine ideas from Alan Turing’s work on pattern formation with May’s random-matrix approach. We demonstrate how a stable equilibrium in a complex ecosystem with trophic structure can become unstable with the introduction of dispersal in space, and we discuss the factors which contribute to this effect. Our work highlights that adding more details to the model of May can give rise to more ways for an ecosystem to become unstable. Making May’s simple model more realistic is therefore unlikely to entirely remove the upper bound on complexity.

In his seminal work in the 1970s Robert May suggested that there was an upper limit to the number of species that could be sustained in stable equilibrium by an ecosystem. This deduction was at odds with both intuition and the observed complexity of many natural ecosystems. The socalled stability-diversity debate ensued, and the discussion about the factors making an ecosystem stable or unstable continues to this day. We show in this work that dispersal can be a destabilising influence. To do this, we combine ideas from Alan Turing's work on pattern formation with May's random-matrix approach. We demonstrate how a stable equilibrium in a complex ecosystem with two trophic levels can become unstable with the introduction of dispersal in space. Conversely, we show that Turing instabilities can occur more easily in complex ecosystems with many species than in the case of only a few species. Our work shows that adding more details to the model of May gives rise to more ways in which an equilibrium can become unstable. Making May's simple model more realistic is therefore unlikely to remove the upper bound on complexity.
It was once a widely-held belief in ecology that greater ecosystem complexity promoted stability [1][2][3][4][5][6]. A priori, this is an intuitive proposition. Surely the more species and connections there are in a food web, the less sensitive the network ought to be to perturbation.
Providing a firm counterpoint to this view, Robert May used a simple statistical model [7,8] to argue that increasing the number of species in an ecosystem could in fact reduce stability. By analysing the eigenvalues of a randomly constructed community matrix, May deduced the following criterion for stability [7] σ 2 N C < 1. (1) In this criterion σ 2 is the variance in the inter-species interactions, N is the number of species and C is the connectance (the probability that any given pair of species interact with one another). May's result shows that stability is decided by the product c = σ 2 N C. We will call this quantity the 'complexity' of the ecosystem. May's model suggests that more complex ecosystems can become unstable. For a fixed variance of interactions, there is an upper bound to the number of species and food web connections that the ecosystem can sustain. This idea quickly became controversial, and May's work sparked the so-called complexity-stability (or diversitystability) debate, which continues to this day [5,6,9].
Since the 1970s, the discussion around stability has been made more precise and subtle. It is now understood that there are a number of senses in which an ecosystem can be unstable, and indeed there are a number of ways one can define diversity [5,10]. An ecosystem can be unstable with respect to, for example, the introduction * joseph.baron@manchester.ac.uk † tobias.galla@manchester.ac.uk of new species, the extinction of existing species, or environmental changes. May's work addresses stability with respect to fluctuations in species abundance. In order to understand the influence of the many aspects of real ecosystems on stability, May's fairly austere model has since been augmented and improved upon. Features not captured by May's initial model include food-web structure (e.g., trophic levels, modularity, and nestedness) [8,[11][12][13], the feasibility of the equilibrium [14,15], variation in interaction strength [16,17], and variability of the environment and of species susceptibility to environmental change [18][19][20][21]. In many models of complex ecosystems only the total abundance of each species is considered, without appreciating how the members of that species are distributed in space [7,8,[11][12][13][14][15][16][17][18][19][20][21]. In such models there is no notion of space, and hence no dispersal. In this work we explicitly include the effects of diffusive dispersal in space, and study the effects on stability.
While dispersal may intuitively be expected to be a homogenising and stabilising influence, Turing showed in his seminal work [22] that it can in fact destabilise a dynamical system. Such an instability has been studied in meta-population predator-prey models with small numbers of species [23][24][25][26]. We combine Turing's idea with May's random-matrix approach to show that a similar destabilising effect can be seen in models of complex ecosystems.
In order not to obscure the key effects at work, we opt to modify May's paradigmatic model sparingly. This allows us to highlight the consequences of the inclusion of dispersal. We suppose that the abundances of the species rest in a steady, homogeneous equilibrium. In order to study stability, we examine the Jacobian matrix governing perturbations about this equilibrium. Like May, we ask what statistical properties are required of the Jacobian matrix in order for the ecosystem to return to equilibrium when perturbed.
Unlike May however, we allow for different trophic The matrix is composed of three parts: a diffusion matrix D, a selfinteraction matrix d, and an interaction matrix A with entries from a Gaussian distribution. Each matrix is split into blocks due to the separation into two trophic levels -we use the subscript u to denote prey species and v for predator species. The approach can be generalised for greater numbers of trophic levels. The figure illustrates the content and structure of the blocks. The stability of the non-spatial ecosystem is described by the matrix for q = 0, or equivalently by setting Du = Dv = 0 (see Methods, and Section S1 in the Supplementary Information). levels, with statistical differences between them. It is the combination of dispersal and trophic structure which gives rise to the Turing-type instability. For the sake of mathematical simplicity, we confine our model to only two broad trophic classifications: predator and prey species. Our approach can be generalised to more complicated food web structures. We now postulate the form of the Jacobian matrix central to our problem, M q . The elements of this matrix describe how spatial disturbances of wavelength λ = 2π/q in the abundance of one species affect the abundances of the other species; q is known as the wavenumber [27].
Because of the separation of the population into trophic levels, the matrix M q has a block structure where each block has different statistics. The matrix is comprised of three terms: a diffusion term, an intra-species interaction term and an inter-species interaction term. We write The diffusion coefficients for prey and predator species are D u and D v respectively (Fig. 1). The interaction matrix A is modelled as having elements drawn from a Gaussian ensemble. Further details on the structure of M q and how one arrives at this form are given in Fig. 1 and in the Methods section. The problem of analysing stability reduces to finding the eigenvalue spectrum of the matrix M q . If the eigenvalues of this matrix all have negative real parts, then the equilibrium is stable with respect to disturbances of wavenumber q. Else, it is unstable. In order for the equilibrium to be stable on the whole, therefore, all eigenvalues of M q must have negative real parts for all values of q.
If the number of species in the ecosystem is large, then the eigenvalue spectrum is dependent only on the statistics of M q and not on its specific entries. Using randommatrix theory and ideas from statistical physics, we are able to deduce a mathematical expression for the support of the eigenvalue spectrum of M q . That is, we can find the region in the complex plane in which the eigenvalues sit and, most importantly, whether or not they have positive real parts. Examples are shown in Fig. 2.
With this analytical approach we can calculate what properties of M q make the equilibrium unstable. Thus, we can deduce how May's upper bound on ecosystem complexity is modified by the inclusion of dispersal and trophic levels. Most crucially, we show that equilibria which would be stable without spatial effects can be destabilised by dispersal. We find that this dispersalinduced instability is possible not only in a linear model, but also in a non-linear system where the equilibrium is arrived at dynamically, and hence feasible by construction.

Results
Eigenvalue spectra. We show some example eigenvalue spectra of the matrix M q in Fig. 2. The vast majority of the eigenvalues group into a 'bulk' region, with the exception of a few outliers. These outliers cannot be ignored -the excursion of even one eigenvalue across the imaginary axis to the positive real side makes the equilibrium unstable.
Using the statistical properties of M q we are able to calculate mathematically the bulk regions to which most of the eigenvalues are confined, and the locations of any outliers. In Figs. 2 (a)-(c) we show that these calcu- Eigenvalue spectra of the stability matrices in Fig. 1. In May's original model, the eigenvalues all lay uniformly within a circle in the complex plane. In our model the circle is warped and split into more complicated shapes. A small number of outlier eigenvalues now stray from the bulk. Eigenvalues of computer-generated random matrices Mq are shown as blue crosses, and compared to theoretical predictions for the boundary of the bulk region (red solid line). Predictions for the outliers are show as open red circles. For q = 0 all of the eigenvalues have negative real part (panel (a)). The equilibrium is stable in the non-spatial ecosystem. For disturbances with larger wavenumber q, one of the eigenvalues crosses the real axis (panel (b)). The equilibrium is unstable against such perturbations. If the wavenumber q is larger still the most unstable eigenvalue returns to the negative half-plane (panel (c)). This is characteristic of a Turing instability [27] -the equilibrium is unstable with respect to disturbances of a finite range of wavelengths. This is shown in panel (d), where we plot the real part of the most unstable eigenvalue as a function of q. The equilibrium is unstable against perturbations of wavenumber q whenever Re[ωmax] > 0. The green triangles mark the wavenumbers from panels (a)-(c).
lations agree very well with the spectra of computergenerated random matrices. We can therefore predict what community properties lead to stability or instability.
As Fig. 2 demonstrates, it is possible to find circumstances under which the model community is destabilised by the inclusion of dispersal. Fig. 2 (a) shows the eigenvalue spectrum for the model without spatial effects (q = 0, or D u = D v = 0, see Methods). All eigenvalues in panel (a) have negative real parts and we conclude that the equilibrium is stable for the non-spatial model ecosystem. In Fig. 2 (b) we take into account dispersal and show the eigenvalue spectrum for a non-zero wavenumber. All other parameters are the same as in panel (a). An outlier eigenvalue strays over the imaginary axis, demonstrating that the equilibrium is now unstable. For perturbations with higher wavenumbers the outlier then returns to the negative half-plane (panel (c)) -the equilibrium is stable with respect to perturbations of higher wavenumber. Fig. 2 (d) shows the real part of the most unstable eigenvalue as a function of the wavenumber q. We see that the equilibrium is unstable against perturbations in a finite band of wavenumbers. In Turing's original work on chemical reaction systems [22], this signalled the formation of stable periodic patterns. The exact shape of these patterns is usually determined by non-linearities in the differential equations describing the reactions. Our model is valid only in the vicinity of the supposed equilibrium, and similar to May [7] we have not specified the nature of any non-linearities. We therefore do not speculate for now about what might happen after the system has departed from the fixed point. We merely point out here the dispersal-induced instability of the equilibrium about which we have linearised.
Modifiying May's bound: stability with and without dispersal. In order to further appreciate the effect that the inclusion of dispersal has on stability, we first consider the conditions under which the non-spatial ecosystem becomes unstable (see Methods). This enables us to study how May's bound on the complexity c changes for our model, which has distinct predator and prey species. A stability plot is shown in Fig. 3 (a). The horizontal axis shows the average degree of predation p = CN µ vu (see Methods), the vertical axis is the complexity parameter c. The solid line indicates the upper bound on the complexity: below the line the equilibrium is stable, above this line it is unstable.
FIG. 3. Dispersal as a destabilising influence. Panel (a) shows the stability diagram for the non-spatial ecosystem model with two trophic levels. The solid line is the upper bound on complexity c for the equilibrium to be stable. Increasing complexity leads to instability. There is a lower bound on the amount of predation p required for the community to be stable (if p is too small, then the equilibrium is unstable even for small values of the complexity). Stability for the spatial ecosystem with dispersal is shown in panel (b). In the blue region the equilibrium is stable. In the yellow region the equilibrium would be stable in a non-spatial model, but dispersal in space induces instability. Crossing from the blue into the yellow region in panel (b) the ecosystem undergoes a Turing instability. In the red region in (b) the equilibrium is unstable both in the non-spatial and in the spatial ecosystem.
We see from Fig. 3 that greater predation p increases the amount of complexity c that can be sustained in stable equilibrium by the ecosystem. Notably, in order to have stability at all there is a lower bound on the predation parameter p.
If we now include dispersal, the stability diagram changes ( Fig. 3 (b)). In particular the upper bound on the complexity can become lower than in the nonspatial system. This is because a new type of instability is now possible -the Turing instability. Thus, there are instances in which the model is stable without dispersal, but unstable when dispersal is introduced (yellow area in Fig. 3 There are no situations in which an unstable equilibrium is stabilised by dispersal (see Section S5 in the Supplementary Information). Dispersal is a purely destabilising influence in our model.
How does complexity affect the Turing instability? So far we have concerned ourselves with the effects of dispersal on complex ecosystems. We now ask the reverse question: Spatial instability and pattern formation have been found in 'simple' models of ecosystems with a small number of species [28][29][30][31]. What are the effects of complexity on this Turing instability?
Turing instabilities in simple systems typically occur when the diffusion coefficients of the activator and inhibitor components are quite disparate -this is known as the 'fine-tuning issue' with the Turing mechanism [32,33]. The activating components in our system are the prey species, predators play the role of the inhibitors.
We would like to know what the effects of complexity are on the threshold ratio D v /D u at which the Turing instability occurs. To answer this question, we compute this threshold for different values of the complexity parameter c.
The case c = 0 can be achieved by setting the width σ of the Gaussian distribution for the elements of the matrix A to zero. All matrix elements within each block are then identical to each other. This means that all predator species become identical to one another, and similarly there is no distinction between different prey species. The model reduces to a simple two-species predator-prey system. Increasing c from zero introduces heterogeneity between the species.
We find that complexity can lower the ratio of diffusion coefficients required for Turing instability (Fig. 4). For more complex ecosystems the Turing instability therefore sets in more easily. So, not only can complexity decrease the stability of non-spatial model ecosystems, it can also reduce stability in spatial models. Conversely, increasing the ratio of diffusion coefficients D v /D u reduces the complexity that can be sustained in stable equilibrium. As can be seen in Fig. 4, the upper bound on c is lower when disparity of dispersal rates is large.
Spatial instability in a non-linear model with complexity. The linear analysis we have focused on so far, although informative, has drawbacks. It only deals with the dynamics in the vicinity of a homogeneous equilibrium, and it tells us nothing about how the ecosystem behaves in the long-run if the equilibrium is unstable. FIG. 4. Complexity reduces the ratio of diffusion coefficients required for Turing instability. In the yellow region the community is unstable due to Turing instability. A high value for the ratio of diffusion coefficients Dv/Du reduces the complexity that can be sustained in stable equilibrium. For sufficiently high complexity the non-spatial community becomes unstable (red area to the right of the vertical line). The equilibrium of the spatial ecosystem is then also unstable, irrespective of the diffusion coefficients.
Further, one could object that the linear model is somewhat contrived and that it does not capture how a 'real' ecosystem constructs itself in equilibrium.
We now present simulation data of a complex ecosystem obeying non-linear Levin-Segel-type dynamics [26] (see Methods). In Fig. 5 we demonstrate that a dispersal-induced instability can occur in this model as well. Panel (a) shows a realisation of the dynamics without dispersal. The model ecosystem converges to an equilibrum composed only of surviving species. By definition this is a feasible equilibrium [34]. Fig. 5 (b) shows the same model (with the same interaction matrix A), but now with dispersal. The abundances do not settle in the long run. Instead, they display quite erratic behaviour. We stress that this is different to the typical behaviour seen in two-species systems with a Turing instability. In such simple systems, the abundance of each species typically settles down eventually. The fixed-point value varies across space, generating a periodic pattern. The complex nature of the interactions in Fig. 5 (b) leads to more complicated dynamics known as diffusion-induced chaos [35,36].
That being said, at any point in time one can take a snapshot of the spatial profile of the species abundances. An example is shown in Fig. 6. One finds some spatial structure for the abundances of prey species (green lines), but it differs from classic Turing patterns, which are typically more periodic and regular. We note that the quickly diffusing predator species (red lines) have more smoothly undulating spatial profiles than the slowly diffusing prey species in this example.

Discussion
The random-matrix approach to modelling ecosystems has developed substantially since May's original work. We contribute to this development with a model of a complex ecosystem which includes different trophic levels and dispersal in space. These additional features change May's bound on complexity.
How the bound changes tells us about the influence of the new model components on stability. For example, we find that predation acts as a stabilising influence (Fig. 3, and Section S5 of the Supplementary Information).
Inspired by Turing's mechanism for pattern formation [22], we show how dispersal can be a destabilising factor in a complex ecosystem. An equilibrium which would be stable in a non-spatial model can become unstable through dispersal. Conversely, we observe that increased complexity can lower the threshold for the diffusion coefficients required for a Turing instability. So, not only does complexity reduce stability in a non-spatial model as was May's conclusion [7], it also destabilises spatial models with dispersal. This is interesting especially in light of the so-called 'fine-tuning' issue with the Turing mechanism. The ratio of diffusion coefficients required for a Turing instability is usually large, making it hard to find experimental examples of Turing pattern formation [27,33,37,38]. Based on our findings, we speculate that Turing patterns may be easier to observe in complex dynamical systems.
To demonstrate that the Turing instability can also be seen in model systems where the equilibrium is arrived at more naturally, we performed simulations of a complex ecosystem with Levin-Segel dynamics (Fig. 5). We found that this model also has a dispersal-induced instability. Because of the complex nature of the interactions the system does not reach a stable patterned state, normally characteristic of Turing instabilities in models with fewer species. Instead, one sees persistently volatile dynamics known as diffusion-induced chaos [35,36].
Our notion of dispersal as a destabilising influence is in contrast to a recent study by Gravel et al. [39] who assert that dispersal can in fact stabilise a complex metaecosystem. One key difference between this work and ours is the way in which dispersal is implemented. In Ref. [39], the ecosystem is organised into patches, and dispersal occurs from any patch to any other patch with equal probability. In other words there is all-to-all connection between the patches, and therefore no sense of spatial arrangement of these patches. In our model, dispersal is local and mediated via a standard diffusion term. We regard local dispersal as a more realistic way to model species movement. The fact that our findings are different from those in [39] indicates that the way one models dispersal can have quite an impact on its effects on stability. We note that the theory we have developed can also be used to predict stability in meta-ecosystems with local dispersal between neighbouring patches (see Section S1 of the Supplementary Information).
One criticism levelled at May's model is that it is too simple and that perhaps through the inclusion of further aspects of natural ecosystems, the upper bound on the complexity would be alleviated. Our results do not sup-port this hypothesis. Like May, we also find that there is always an upper limit on the complexity c = σ 2 N C that an ecosystem can stably sustain. Other recent studies using random-matrix approaches arrive at similar conclusions. For example, Allesina and Tang take into account more realistic food-web structures and still find upper limits on the number of interconnected species [12,13]. May's result therefore generalises to models capturing more aspects of 'real' communities in ecology. This prediction is supported by the observation of 'diversity regulation' in some ecosystems [40][41][42][43][44].
A final observation that we wish to convey is that making models for complex ecosystems more detailed introduces the opportunity for new types of instability. In May's original model, for example, the mean of the community matrix elements was zero. As a consequence, any one species is equally like to benefit or suffer from the presence of another species. Mathematically, all eigenvalues then reside within one bulk region, and it is this bulk region that determines stability. Mutualism can be introduced through interaction coefficients which are positive on average, and competition through a negative average interaction [12]. This leads to additional outlier eigenvalues, which can make an equilibrium unstable even though it would otherwise be stable. Introducing trophic levels can generate complex-conjugate pairs of outliers (Fig. 2  a)), allowing for more opportunity for instability (Section S5 of the Supplementary Information). Dispersal, finally, leads to the possibility of a Turing-type instability. Overall, adding more details to the model of May tends to give rise to more ways in which an equilibrium can become unstable. Making May's simple model more realistic is therefore unlikely to remove the upper bound on complexity.
-Methods -Linear model. We imagine that we find the ecosystem at an homogenenous equilibrium. Our model is concerned with the dynamics of small perturbations of the species abundances about this fixed point. The stability of the homogeneous fixed point is determined by whether or not these perturbations decay or increase with time. We write u i (x, t) and v j (x, t) for the perturbations of the prey and predator species abundances respectively at position x and time t. These are the deviations away from the fixed point. There are N u prey species and N v predator species with N = N u + N v species in total. We define the constants γ u = N u /N and We assume that prey species diffuse at rate D u and predator species at rate D v in a spatially homogeneous environment. Similar to May [7], we also imagine that all species have negative self-interaction. This is a necessary condition for stability and reflects intra-species competition. The probability that a particular pair of species interact with one another is C. This parameter is known as the 'connectance' [7]. The effect of a change in the abundance of species j on species i, where j belongs to trophic level β and i belongs to α, is A αβ ij [α, β ∈ {u, v}]. The linearised reaction-diffusion equations thus take the following form If species i and j are non-interacting, A αβ ij = A βα ji = 0. If the two species interact, then A αβ ij and A βα ji are drawn from a joint Gaussian distribution with means A αβ ij = µ αβ and A βα ji = µ βα . The elements in the random matrix have variance σ 2 and they are correlated according to All other correlations are set to zero. Each of the model parameters can be interpreted ecologically. We define and interpret the complexity c = CN σ 2 in the main text. The interaction mean µ uu indicates the degree to which different prey species cooperate (if µ uu > 0) or compete (if µ uu < 0). The coefficient µ vv has a similar interpretation for predators. The means of the off-diagonal blocks µ uv < 0 and µ vu > 0 indicate the degree to which (on average) prey species suffer and predator species gain from predator-prey interactions. The parameters Γ αβα β describe the correlations between interaction coefficients. The only non-zero entries are taken to be Γ u ≡ Γ uuuu , Γ v ≡ Γ vvvv and Γ uv ≡ Γ uvvu . That is, only elements which are diagonally opposite one another in A are correlated. A positive value of Γ αβα β indicates that if one species benefits more than average from an interaction, the other species involved does so as well. The opposite is true if Γ αβα β is negative.
Taking the Fourier transform with respect to the spatial coordinate x of Eq. (3), we arrive at dynamical equations (see Supplementary Information Section S1) for disturbances of wavenumber q in the abundances of the various species (the wavenumber is related by q = 2π/λ to the wavelength λ). We denote the combined vector of the Fourier transforms of species abundances byX q = (· · · ,ũ i (q, t), · · · ,ṽ j (q, t), · · · ), and arrive at the more compact matrix equationẊ The matrix entry (M q ) αβ ij tells us what the effect of a disturbance of wavelength q in species j (belonging to trophic level β) is on species i (belonging to trophic level α).
The vector X q has dimension N = N x + N y and is arranged such that the first N x elements are the Fouriertransformed abundancesũ i (q, t) and the last N y elements are theṽ j (q, t). The matrix M q is depicted in Fig. 1. It has a block structure due to the separation into trophic levels. Its three contributions are a diagonal diffusion matrix, a diagonal self-interaction matrix and a random interaction matrix, whose variance and correlations are given in Eqs. (4) and (S2). The indices α and β correspond to the different blocks and i and j correspond to the position within the block. If the matrix M q has eigenvalues with positive real parts for any value of q ≥ 0, the disturbances u i and v i will grow with time, indicating an unstable equilibrium.
By setting q = 0 in Eq. (6), one recovers Eqs. (3) with the diffusion term removed. Focusing on q = 0 in our model (or equivalently setting D u = D v = 0) thus allows one to study the stability of a non-spatial model ecosystem with two trophic levels. See also Sections S1 and S5 in the Supplementary Information.
We note that the same form for the matrix M q can be found for a model of a meta-ecosystem (with dispersal between nearby patches instead of in continuous space), see Section S1 of the Supplementary Information.
The values of the model parameters used in the figures are given in full in Section S6 of the Supplementary Information.
Calculation of the boundary surrounding the bulk of the eigenvalues. The vast majority of the eigenvalues of the random matrices M q reside within one or two 'bulk' regions of the complex plane. To determine stability we need to know if there are bulk eigenvalues with positive real parts. Identifying the boundaries of the regions containing the eigenvalues is sufficient for this purpose. Our calculation uses methods originally developed in condensed matter physics, and follows lines similar to those of [45,46].
We write ω = ω x + iω y for the eigenvalues of M q . The general expression for the boundary surrounding the bulk of the eigenvalues (ω y as a function of ω x ) is given by the simultaneous solution of the following equations (see We note the free index α in the second of these equations (α ∈ {u, v}). This is therefore a system of three coupled equations. One first eliminates the auxiliary variables χ α , and then expresses ω y in terms of ω x . This results in the red curves in Fig. 2.
The solution simplifies in several special cases, which we exploit to provide explicit stability criteria analogous to May's bound (Section S5 of the Supplementary Information).
Calculation of outliers. In addition to the bulk eigenvalues the stability matrix can have isolated outlier eigenvalues. If any of these outliers have positive real part, the equilibrium is unstable. Their position in the complex plane is calculated along the lines of [47]. Details can be found in the Supplementary Information. The outlier eigenvalues are given by the complex values ω satisfying the following equation (see Supplementary Information Eq. (S82)) The auxiliary quantities χ α (ω) in this relation satisfy subject to the condition Eqs. (S20) and (S21) need to be solved simultaneously, subject to Eq. (S22). If there are no solutions then there are no outliers. In special cases the above expressions can be simplified, and explicit stability criteria can be found (Section S5 in the Supplementary Information).
Finding the threshold for instability. Instability can occur in one of several different ways: (1) The bulk region of eigenvalues for q = 0 can cross into the positive half-plane; (2) One of the outlier eigenvalues for q = 0 can cross the imaginary axis; (3) An outlier eigenvalue for q = 0 can stray into the positive half-plane. We have not observed any circumstances under which the bulk crosses the imaginary axis for non-zero q where it does not for q = 0. The eigenvalues must have negative real parts for all q in order for the equilibrium to be stable in the spatial system. This includes q = 0. Cases (1) and (2) therefore indicate instabilities occurring both in the non-spatial and the spatial ecosystem. In case (3) the spatial system is unstable, but the non-spatial system remains stable. In each of these cases, the threshold for instability is found by identifying sets of parameters for which either the boundary for the bulk eigenvalues touches the imaginary axis (case 1) or where the outlier eigenvalues touch the imaginary axis (cases 2 and 3).
This can be done using our analytical results for the spectrum of eigenvalues (Section S5 of the Supplementary Information), leading to the results shown in Figs. 3

and 4.
Simulating the non-linear model. Results in Figs. 5 and 6 are from a numerical integration of the Levin-Segel-type model [26], where a > 0 is a constant. To integrate these equations numerically, the diffusion terms are discretised. The integration is then carried out using the Runge-Kutta (RK4) method [48].

A. General setup: trophic levels, species interactions and diffusive dispersal
We suppose that the total population consists of a sub-population of predator species and a sub-population of prey species, which are distributed in space. We assume that there are N u = γ u N different prey species and We imagine that the system reaches an equilibrium and that the equations governing the dynamics of the species abundances, linearised about this equilibrium, take the following form The u i (x, t), i = 1, . . . , N u , describe perturbations of the abundances of the prey species about the equilibrium at position x and time t, and similarly v j (x, t), j = 1, . . . , N v , are perturbations of the predator abundances.
We presume that the interactions within a species are competitive (such that d u > 0 and d v > 0), and that each species disperses via diffusion. The corresponding dispersal coefficients are D u for prey, and D v for predators. Crucially, the prey and predator species can have different dispersal coefficients D u = D v . In the model there are inter-population and intra-population interactions. These are described by the interaction coefficients A αβ ij , where α and β can takes values α, β ∈ {u, v}, where u stands for prey and v for predator. For example, A uu ij describes interactions between different prey species, and A uv ij describes the effect that predator species j has on prey species i. The interaction coefficients A αβ ij are random numbers drawn from a Gaussian distribution with All other correlations are zero. We assume that the only non-zero Γ αβα β are Γ uuuu ≡ Γ u , Γ uvvu ≡ Γ uv and Γ vvvv ≡ Γ v . With a probability 1 − C, species i and j are chosen not to interact with one another and thus A αβ ij and A βα ji are set to zero. The parameter C therefore describes the 'connectance' (the probability that any given pair of species interact with one another).
Each of the model parameters can be interpreted ecologically. We define and interpret the complexity c = CN σ 2 in the main text. The interaction mean µ uu indicates the degree to which two different prey species cooperate on average (if µ uu is positive) or compete (if µ uu is negative). The coefficient µ vv has a similar interpretation for predators. The means of the off-diagonal blocks µ uv < 0 and µ vu > 0 indicate the degree to which (on average) prey species suffer and predator species gain from predator-prey interactions. The coefficients Γ αβα β indicate statistical correlations between interaction coefficients. A positive value of Γ αβα β indicates a statistical tendency for one species to benefit more (than average) from an interaction if the other species involved does so as well. The opposite is true if Γ αβα β is negative. The only non-zero entries are taken to be Γ u ≡ Γ uuuu , Γ v ≡ Γ vvvv and Γ uv ≡ Γ uvvu . That is, only elements which are diagonally opposite one another in A are correlated [49,50].

B. Fourier transform and stability matrix
The perturbations u i (x, t) and v j (x, t) in Eq. (S1) are functions of position (and time). Taking the Fourier transform with respect to the spatial coordinate x, we writẽ and similar forṽ j (q). Theũ i (q),ṽ j (q) describe perturbations of wavenumber q, or equivalently, of wavelength λ = 2π/q. Using this in Eqs. (S1) one obtainṡ We note that by setting q = 0 in Eqs. (S4) one recovers Eqs. (S1) with the diffusion term removed. From the definition of the Fourier transform Eq. (S3), we see these are then equations for the total u i and v j summed throughout space. This removes any sense of spatial distribution of species. Setting q = 0 therefore recovers the naïve 'non-spatial' model where species dispersal is not taken into account. Eqs. (S4) can be written in a more compact matrix forṁ where the vector X is of length N , and with The matrix M q is of dimension N × N . It is the sum of three terms: (i) The first term is −q 2 D, where the matrix D is of size N × N , diagonal and of the form where 1 1 Nu and 1 1 Nv are identity matrices of dimensions N u × N u and N v × N v respectively.
(ii) The second term is of a similar form and describes intraspecific interaction.
(iii) The third term in M q is a random matrix with the Gaussian statistics given in Eq. (S2). We note that the matrix M q has a block structure; it is comprised of submatrices of dimensions N u × N u , N u × N v , N v × N u and N v × N v . This is illustrated in Fig. 1 of the main text.
The eigenvalues of the matrix M q tell us about the stability of the hypothetical fixed point about which we have linearised. More precisely, the fixed point is stable against perturbations of wavelength λ = 2π/q if and only if all eigenvalues of M q have negative real part. The fixed point is stable against perturbations of any wavelength only if this is the case for all q. That is to say, in order for the equilibrium to be stable, the matrix M q must not have an eigenvalue with positive real part for any q.
We therefore wish to identify the eigenvalue spectrum of M q . We can then use this to deduce under what conditions the equilibrium is stable. The matrix M q has N u + N v eigenvalues in total. The eigenvalue spectrum consists of a 'bulk' of eigenvalues with some possible outliers (see Fig. 2 in the main text). The bulk contains the vast majority of eigenvalues. They are confined to a compact area in the complex plane. Both the boundary of the bulk region and the location of the outliers can be calculated based only on the statistical properties of the matrix M q (as opposed to the exact entries). Further, the boundary of the bulk region can be calculated independently of the outliers. We therefore split our calculation into two parts, and present the calculation of bulk eigenvalues in Section S2, and that for the outliers in Section S3.

C. Correspondence with a model of a meta-ecosystem with patches
We now briefly consider a model in which species hop between discrete patches, as opposed to diffusion in a continuous space. We label patches by an integer index , and assume that patch is located at position ∆x in space. That is to say, the spacing between two adjacent patches is ∆x. The continuous diffusion term in Eqs. (S1) is then replaced to yield (We have suppressed the dependence on time to keep the notation compact). The quantity φ( , ) is a 'hopping kernel', describing the probability that an individual hops from patch to patch , given that a hop occurs. We suppose hopping is only between neighbouring patches, and occurs with equal probability. We therefore choose The constantsD u andD v in Eqs. (S9) describe the overall hopping rates (for prey and predators respectively) to any of the two neighbouring patches. If the migration between patches is diffusive (or just by dimensional arguments) these rates must be proportional to (∆x) −2 , so we Next, we take the discrete Fourier transform with respect to the spatial variable in Eqs. (S9) , We obtain withφ (q) = cos (q∆x) . (S13) Provided the arrangement of patches is 'dense', i.e., that there are many patches per unit distance, we can assume ∆x 1, and expandφ(x) = cos (q∆x) in powers of ∆x. To lowest non-trivial order we find Eqs. (S12) then become and we thus arrive at Eqs. (S1). This shows that any stability criteria derived for the ecosystem with diffusion in continuous space also applies to a model with discrete patches and hopping between adjacent patches. The approximation is valid provided the distance between two patches in space is small (∆x 1), and we have assumed the standard scalingD u = D u /(∆x) 2 ,D v = D v /(∆x) 2 . It is also possible to analyse the model with discrete patches, and to find the eigenvalues as a function of q without performing this approximation. One would still find the possibility for Turing instability (as demonstrated for a different system in [33]).

S2. BOUNDARY OF THE BULK OF THE EIGENVALUE DISTRIBUTION
In this section, we broadly follow the structure of similar calculations in [45,46,51].

A. Preliminaries -setting up the calculation
For an ensemble of N × N random matrices m we consider the so-called resolvent, defined by The overbar represents an average over the ensemble of matrices. The expression 1 1 N denotes the identity matrix of size N × N . We note the following identities from general linear algebra (for an N × N matrix B): (i) Tr P −1 B P = Tr [B], for any invertible N × N matrix P; (ii) If B is invertible, then the eigenvalues of B −1 are reciprocals of those of B; (iii) The trace is the sum of the eigenvalues; (iv) Adding the identity matrix to a matrix increases each eigenvalue of this matrix by 1; (v) Multiplying a matrix by a constant scales its eigenvalues by that same constant.
Using these, one can show that where the λ are the eigenvalues of m. In the limit N → ∞, this sum can be written as an integral over the complex plane where ρ(λ) is the density of eigenvalues. Bearing in mind that G(ω) is a complex quantity, we perform a contour integral around a closed path C, which we presume contains none of the poles of G, In this expression S denotes the region bounded by C, and we have used the residue theorem to evaluate the contour integral. We next note the complex version of Gauss' law (letting ω = x + iy) Defining the real-valued potential Φ(ω) via and using Eqs. (S4) and (S5), along with the fact that the region S is arbitrary, we obtain Poisson's equation In analogy with electrostatics, the resolvent G(ω) can be thought of as being related to the 'electric field' with components E x = 2 Re G and E y = −2 Im G [46]. Therefore we can obtain the distribution of eigenvalues, if we can find the potential Φ. The potential Φ can be related directly to the matrix m via This can be seen to agree with Eqs. (S6) and (S2) once one realises that det (1 1 N ω − m) = λ (ω − λ).
We note the following about this potential and its relationship with the eigenvalue density. By definition, the resolvent (or average Green's function) is given by Indeed, we also have The resolvent is related to the eigenvalue density by The left-hand side of this equation is familiar from the Cauchy-Riemann equations for analytic functions. From this we deduce that if G(x, y) is an analytic function in a particular region (i.e. we can write it in terms of ω only), the eigenvalue density must be zero in this region.

Ensemble of random matrices
Adding a low-rank perturbation to a large random matrix is known not to change the bulk of the eigenvalue spectrum [12,47,52]. This means that we can assume that the mean of the distribution from which the matrix elements are drawn is zero. The mean can have an effect on outlier eigenvalues, and we will therefore restore it for the calculation of the outliers in Section S3.
For the time being we consider a Gaussian random matrix with statistics as follows The matrix m is similar to M q (the matrix we focus on in the main paper), but it is different in a number of ways: First, we note that we have re-scaled the matrix elements relative to the matrix M q in the main paper -the variance and co-variances are now proportional to 1/N . This approach is standard in the theory of random matrices, its objective is of a technical nature. The re-scaling allows us to carry out the calculation in the limit N → ∞. In practice the results derived in this way are often accurate also for finite N . The scaling is undone by replacing σ 2 → N σ 2 in the result. For a more detailed discussion, see Section S4.
Further, we are assuming all-to-all interaction (C = 1) for now between the species. We discuss how choosing a general connectance 0 < C < 1 affects the result in Section S4.
Finally, the wavenumber q is set to zero here. The result for non-zero q can obtained by making the replacement d α → d α + q 2 D α , see also Section S5 B.

Replica method and Hubbard-Stratononich transformation
In order to calculate the spectrum of the ensemble of random matrices we need to compute ln det [(ω * 1 1 N − m T )(ω1 1 N − m)]. We recall that the overbar stands for the average over the randomness in m. In order to calculate averages of logarithms, one can use the so-called replica method [46,[52][53][54], based on the identity ln y = lim n→0 y n −1 n . This allows one to calculate y n instead of the more intricate average ln y. The object y n in turn represents an n-fold 'replicated' system in statistical physics [46,[52][53][54]. In the context of the present problem one finds that the n replicas decouple in the limit of large N , similar to what was observed in [45,46]. This means that (S13) Using a Gaussian integral we can write the determinant in Eq. (S13) in the form where is a positive infinitesimal quantity introduced so as to avoid singularities which occur when ω is equal to one of the eigenvalues of m. To carry out the average over the ensemble of random matrices we need to evaluate We next perform a Hubbard-Stratonovich transformation [55] in order to linearise the terms involving m αβ ij . This yields Average over the random matrix and introduction of order parameters Carrying out the average over the random matrix elements in Eq. (S16) then gives [recalling the definition of the potential Φ(ω) in Eq. (S8)] (S17) Following [45,52] we disregard contributions containing N −1 and N −1 i z α i y α i . These do not contribute to leading order in N . We now introduce the order parameters which we impose in the integral Eq. (S17) using Dirac delta functions in their complex exponential representation. For example where we note that w α is a complex quantity. We can thus rewrite Eq. (S17) as where where D [· · · ] denotes integration over all of the order parameters and their conjugate ('hatted') variables, and where We recall that α can take the values α = u (prey) and α = v (predators), and that γ u and γ v are the fraction of prey and predator species respectively. We note that the integrals over y α i and z α i are uncoupled for different values of i as a result of introducing the order parameters in Eq. (S18). Carrying out the integrals over the variables y α and z α in the expression for Ω one obtains Ω = α γ α ln(ŵ αŵ α −ûv). (S22)

Saddle-point integration and evalulation of order parameters
We now take the limit N → ∞, and carry out the integral in Eq. (S20) in the saddle-point approximation. To do this, we extremise the expression Ψ + Θ + Ω. Extremising with respect to the conjugate variablesû,v,ŵ α andŵ α , we find One makes the following observation Further, we find that which means that we can now expressŵ αŵ α in terms of w α w α andûv. We also find that which means that we can in turn expressûv in terms of w α w α and uv. Defining r = uv, we can writeûv andŵ αŵ α in terms of w α w α and r. This means that we can replaceŵ αŵ α −ûv in the expression for Ω. So, lettingq α = (ŵ αŵ α −ûv) −1 , we obtain an implicit equation forq α in terms of r and w α w α : We now extremise Ψ + Θ + Ω in Eq. (S20) with respect to u, v, w α and w α . We find Combining these with Eqs. (S23), one sees that The first equation in this set has two solutions for → 0 + .

First solution:
One solution is v = 0, implying r = 0, andq α = −w α w α . The last two of the above equations yield From this, one can solve for the resolvent G(ω) = α γ α w α , which we see is an analytic function of ω. This means that the eigenvalue density vanishes in regions of the complex plane for which r = 0 is the only valid solution.
Second solution: The other solution to the first of Eq. (S29) is This can be solved simultaneously with the following equations in order to find r, w α and w α as functions of ω = x+iy

Significance of the order parameters
It is important at this point to note that the order parameters defined in Eq. (S18) have some broader significance. The quantities w α and w α , when evaluated at the saddle point, are related to the resolvent: if we had redone the calculation up until this point but had replaced ω1 1 N in the definition of the resolvent Eq. (S1) with a diagonal matrix (ω ω ω) αβ ij = ω α δ αβ δ ij ]. Then, instead of arriving at Eq. (S33), we would have been able to define The result of this line of reasoning is that, at the saddle point, we deduce iw α = (iw α ) . Such a deduction is still valid once we set ω α = ω. This observation will prove important for ruling out invalid solutions when we solve the saddle-point equations.
Although the fact that at the saddle-point iw α = (iw α ) may at first seem troubling, one should note the following. Let w α = q α + is α and w α = q α − is α . The (real) variables over which we integrate in Eq. (S20), in order to impose the delta function in Eq. (S19), are q α and s α . Because the argument of the exponential in Eq. (S20) is complex, we have to use analytical continuation to deform the integrals, which are initially over the real axis, to contours in the complex plane. These contours include the saddle point. As such, the ostensibly 'real' variables q α and s α take on complex values at the saddle point. Hence, w α and w α are no longer necessarily complex conjugates when evaluated at the saddle point.

Calculation of boundary of bulk spectrum
Inside the support of the eigenvalue spectrum (where the eigenvalue density is non-zero), we have r > 0 and the condition Eq. (S31). Assuming no discontinuities in the value of r as a function of ω, the boundary of the support of the eigenvalue spectrum is defined by points ω in the complex plane where r approaches zero, but where Eq. (S31) is still satisfied. This yields the following set of equations which one can solve for the curve defining the boundary ω y (ω x ) (where ω = ω x + iω y ): That is, values of ω which lie on the boundary of the support of the eigenvalue spectrum must satisfy the above simultaneous equations.

A. Setup and general relations
We now wish to find the outliers that arise when we perturb the random matrix m by introducing a finite mean. That is, we are now interested in the matrix m ≡ m + 1 N µ, where the N ×N matrix µ accounts for the non-zero mean. Specifically m has elements in particular independently of i and j. The elements of m have zero mean as before. So, the elements within each of the 4 blocks [(α, β) = (u, u), (u, v), (v, u), (v, v)] of the matrix m have the same mean, µ αβ /N . Following the conventions of statistical physics, the mean is chosen to be proportional to 1/N for technical convenience. In the main paper, the mean is independent of N . We explain the mapping between the two setups in Section S4, where we also discuss general values of connectance 0 < C < 1.
The problem amounts to finding the eigenvalues of the matrix m -i.e. m with a rank-2 perturbation µ. We follow a procedure similar to Ref. [47]. The bulk of the eigenvalue spectrum is unchanged by the low-rank perturbation. Therefore we look for additional eigenvalues that lie outside of the bulk of eigenvalues of m. Such an eigenvalue λ must satisfy Any outlier λ, is by definition, not an eigenvalue of the unperturbed matrix m. Therefore the matrix 1 1 N λ − m is invertible. Hence,

B. Simplification to an effective 2-species problem
Now we define the resolvent matrix G = (λ1 1 N − m) −1 .
We suppose for now that there are the same numbers of prey species as predator species, i.e. N u = N v . We will relax this assumption once we have outlined the proof for this simpler case. Suppose we rearrange the rows and columns of the matrices so that the first species is a predator species, the second a prey species, the third again a predator species, and so on. The matrix µ then consists of N u × N u block of size 2 × 2 of the form µ µ µ (2) ≡ µ uu µ uv µ vu µ vv . Then, the matrix µ can be written in the form where u and v are vectors of length N u . Each entry of v is given by the 2 × 2 matrix µ µ µ (2) , v T = µ µ µ (2) , µ µ µ (2) , · · · , The vector u is of the form We next use Sylvester's determinant identity valid for combinations of m × k matrices A and k × m matrices B. We note that the matrices on the left-hand side of Eq. (S8) are of size m × m, and those on the right-hand-side of size k × k.
Using this identity, we find that det(1 where General case (N u not necessarily equal to N v ): We now wish to show that Eq. (S9) also applies in the case N u = N v . This is done using the following identity for block matrices: If A, B, C and D are matrices of sizes m × m, m × k, k × m and k × k respectively, then we have (S10) First we note that the determinant det 1 is unchanged by adding extra rows and columns to the matrices in the following way: We add rows and columns µ to produce a square matrix µ N of dimension N = N + |N u − N v | which can be decomposed as in Eq. (S5). We then augment G to obtain Then, using the theorems in Eq. (S10), one can show Hence, we can apply Sylvester's determinant identity as before and arrive at Eq. (S9).

C. Proof that the resolvent matrix is diagonal
We have shown that Eq. (S9) holds in all cases. It remains now to find the averaged Green's functions . For this, we follow [51] and use the so-called 'locator expansion' to show that the resolvent matrix is diagonal in the limit N → ∞, i.e. that G αβ from which we obtain a Dyson equation G = f + Gmf . This allows us to construct a series for the resolvent where we have omitted sums over repeated indices for brevity. We find that we have to evaluate terms such as k,γ m αγ ik f γ k m γβ kj ≈ kγ f γ k m αγ ik m γβ kj . This involves averages over products of Gaussian variables. Such terms can be simplified using standard Gaussian combinatorics which is valid for even n and is zero otherwise, and where stands for the sum of the products of all possible one can see that (since the only free indices are i, j, α and β in Eq. (S13)) that G αβ ij ∝ δ ij δ αβ . We therefore find that for large N .

D. Final expression for outliers
Given that the resolvent matrix is approximately diagonal for large N , we only need to calculate the sum of diagonal elements 1 N i G αα ii in Eq. (S9). This can be found from the quantities involved in the bulk calculation as follows.
We know that at the saddle point . Now, differentiating with respect to ω α and using Jacobi's formula for the derivative of a determinant, we have Now setting ω α = ω, we see from Eq. (S9) that the outliers λ obey This can be written as The quantities χ α (λ) satisfy [c.f. Eq. (S36)] Note that not all solutions of this set of equations are valid. At the saddle point, the order parameter u = 1 N iα z α i z α i must satisfy u ≥ 0, and we have u = (1 + σ 2 u) α γ α |χ α | 2 . This means that we must have in order for the set (λ, {χ α (λ)}) to be a valid solution to the combination of Equations (S20) and (S21). Eqs. (S20) and (S21) can be solved numerically for λ, subject to the constraint in (S22). Note that when the left-hand side is equal to the right-hand side in the inequality (S22) the point λ in question is inside the bulk region.  N ; (b) Not all species interact with all others in the main text (i.e. the connectance can take values smaller than one, 0 < C < 1). We now describe how to relate the results of the previous sections to the setup in the main text.

A. Rescaling with N
In this Supplement, we used the statistics given in Eqs. (S12) and (S1) for the entries of the random matrix m. These are The statistics used in Eq. (4) and (5) of the the main paper (and also given in Eq. (S2) of this Supplement) are different. They can be obtained from the above by making the replacements µ αβ → N µ αβ and σ 2 → N σ 2 . The results in the main paper have been obtained from those in Sections S2, S3 and S5 by making this replacement.
The theoretical predictions for stability or instability are formally derived in the limit N → ∞, but in practice they are accurate for finite values of N , provided N is not too small. Sizes of N ≈ 50 are often sufficient to obtain satisfactory agreement (see Fig. S1 below for an example). The spectra in Fig. 2 of the main paper are obtained numerically from matrices of size 750 × 750, and can be seen to agree very well with the results of the mathematical theory.
B. Taking into account connectance C A connectance C < 1 between species can be accounted for by making a similar replacement of σ 2 and µ αβ to the one described above. We now describe how to see this.
In order to find the potential Φ(ω) (which contains all the information about the eigenvalue spectrum and could in principle also be used to find the outliers [52]) defined in Eq. (S8), one needs to carry out an average over the ensemble of random matrices, see e.g. the step from Eq. (S16) to Eq. (S17). We now wish to take into account that any one pair of species only interacts with probability C. We then obtain [c.f. Eq. (S15)] Upon evaluating the integral, one obtains If we set C = 1, we recover the term we would normally obtain with all species connected For individual sets of i, j, α and β, and assuming N 1, the expression in the curly brackets in Eq. (S3) can be expanded in orders of 1/N to obtain Comparing Eq. (S3) [with Eq. (S5) substituted in] to Eq. (S4), we see that the consequence of introducing a connectance 0 < C < 1 is to scale the system parameters in the following way (S6)

C. Combined effect
In summary, the scaling that one needs to perform in order to obtain the results in the main text from those in Sections S2 and S3 is given by All simulation results in the main text were performed using the model as defined in Section S1 (which is the model defined in the main text). We note that Eqs. (S36) and (S20-S22) become respectively and where c = N Cσ 2 .

S5. STABILITY CRITERIA
We class the instabilities in our system into two types: (1) The bulk of the eigenvalue distribution may cross the imaginary axis. This is the type of instability May identified [7].
(2) One of the outlier eigenvalues may cross the imaginary axis, with the bulk staying in the negative half-plane. Instabilities of this type are discussed for non-spatial systems for example in [13].
Below, we provide analytical criteria for each of these instabilities. Crucially, we note that an instability of the second kind can be introduced to an otherwise stable system by the inclusion of diffusion. That is, for a finite band of wavenumbers, the spatial system is unstable where the non-spatial system would be stable. First, we provide criteria for the stability of the non-spatial system, and then we characterise the spatial instability. In the following the model parameters σ 2 and µ αβ are scaled in the same way as in the main text, and we allow for general connectance 0 < C < 1.
A. Instabilities of the non-spatial system Instability due to bulk eigenvalues, and relation to May's bound Mathematical condition for instability: The transition occurs when the boundary of the bulk eigenvalue spectrum first hits the imaginary axis. We generally observe that the bulk of the eigenvalue spectrum is a convex set. Since the matrix we are dealing with is real, we also know that the spectrum is symmetric with respect to the real axis. As a consequence we find that point of first contact is at ω = 0. We note that, from the correspondence with the resolvent shown in the previous section, this means that the variables χ α (0) must be real. In order for the system to be stable, λ = 0 must necessarily be outside the bulk (otherwise there are positive real eigenvalues). Hence, a necessary condition for stability is We have used Eq. (S22) with the rescaling µ αβ → Cµ αβ , and σ 2 → Cσ 2 , and we recall that the complexity parameter is given by c = N Cσ 2 . The quantities χ α satisfy The equations above cannot be solved analytically in general, but a numerical solution can be obtained.

Special cases:
We can make the connection with May's original criterion [7] in several special cases. First, consider the case where γ v = 0 -that is, we do not distinguish between predator and prey species, and consider a general population of N = N u species (as May did). In this case, Eqs. (S1) and (S2) yield We recover May's criterion for Γ uu = 0 and d u = 1. Secondly, the criterion for stability also simplifies greatly when the prey and predator populations have uncorrelated internal interactions so that Γ u = Γ v = 0, the same number of species (i.e., γ u = γ v = 1 2 ), and equal death rates d u = d v = d. In this case, Eqs. (S1) and (S2) give Similarly, if Γ uv = 0 and Γ u = Γ v = Γ, we obtain A further special case is a situation in which all interactions are uncorrelated (Γ αβ = 0 for all α and β). Here we find From these examples, we get a general idea of how each of the model parameters affect the onset of instability. The intra-species interaction I α = −d α is required to be negative in order to have stability. The more negative I α is, the more stable the equilibrium is. The equilibrium can also be stabilised by asymmetry in the interaction coefficients. That is, the lower the parameters −1 ≤ Γ αβ ≤ 1 are, the more stable the equilibrium. Finally and crucially, as per May's conclusion, a lower complexity c = N Cσ 2 contributes to making the equilibrium more stable in all the special cases outlined above.

Instability due to outlier eigenvalues
Mathematical condition for instability: If we suppose that the bulk eigenvalues remain to the left of the imaginary axis, an instability may still occur if one or more of the outlier eigenvalues take on a positive real part. This may occur in two possible ways: (I) A single outlier with no imaginary part may cross at λ = 0; or (II) A pair of complex conjugate outliers may cross the imaginary axis. The point at which an instability of the first kind occurs can be obtained by substituting λ = 0 in Eqs. (S20) and (S21). Re-scaling σ 2 and µ αβ to match the conventions in the main paper, we obtain The model parameters at which an instability of type (II) above occurs are found as follows: We evaluate the outlier eigenvalues by solving Eqs. (S20) and (S21) subject to Eq. (S22), and we carry out the re-scaling µ αβ → N Cµ αβ , σ 2 → N Cσ 2 . We do this for varying model parameters until one of the outliers gains a positive real part. This way, we ensure that we spot any complex conjugate pairs crossing the imaginary axis. The sloped lines in Fig. 3 in the main text, for example, were generated in this way.
Example: Fig. 2 in the main text shows an example in which a single outlier crosses the imaginary axis. In Fig. S1 we show a case in which a complex-conjugate pair of outliers has crossed the axis. The model parameters are the same as in Fig. 2 of the main text, with the exception of the µ αβ . This shows that altering the means of the matrix elements can turn an instability of type (I) above into one of type (II). We note that the bulk of the eigenvalue spectrum is the same in Fig. 2a and in Fig. S1. This confirms that the bulk spectrum is not affected by changes in the means µ αβ .

Special cases:
We can make analytical progress in the same set of special cases as in the previous section. First, consider the case γ u = 0 (we could similarly use γ v = 0). Then, the outlier eigenvalue [from Eq. (S21)] is real and given by For C = 1 this reduces further to results mentioned in [47,52]. A necessary condition for stability is then In the case where the prey and predator populations have uncorrelated internal interactions so that Γ u = Γ v = 0, the same numbers of species γ u = γ v = 1 2 , and the same death rates d u = d v = d and µ uu = −µ vv , the outlier eigenvalues are given by If µ 2 uu + µ uv µ vu < 0, there is no real solution for λ and there is guaranteed to be no instability, so long as d > 0. We note that complex conjugate pairs of outliers are made possible only by having different trophic levels.
Similarly, in the case Γ uv = 0 and Γ u = Γ v = Γ, the eigenvalues are given by λ = −d ± (N C) 2 (µ 2 uu + µ uv µ vu ) + 2Γc 2N C µ 2 uu + µ uv µ vu . (S11) In both special cases, the more negative the quantity µ uv µ vu is, the more stable the equilibrium. We note that −µ uv µ vu is a measure for the 'strength' of predation (taking into account both the effects on prey and predators). So, predation is a stabilising influence. Further, asymmetry in the interaction coefficients has a stabilising effect. A further special case is Γ αβ = 0 for all α and β. Here, we obtain Once again, the more negative the quantity µ uv µ vu is, the more stable the equilibrium.

B. Dispersal-induced instability
So far we have considered only a non-spatial ecosystem in this Supplement, i.e. we have set q = 0 in Eq. (S4), or equivalently D u = D v = 0. We now address instabilities induced by dispersal. The additional terms arising from a non-zero dispersal term can be absorbed into a re-definition of d α , (S13) This can be seen directly from Eq. (S6). With the replacement in (S13), we can then use the stability criteria derived in Sec. S5 to obtain criteria for dispersal-induced instability.
The instability due to bulk eigenvectors is unaffected by dispersal With the inclusion of dispersal Eqs. (S1) and (S2) become α γ α χ 2 α < 1 c , This is the condition for the stability of the Fourier mode with wavenumber q.
Although it is difficult to prove analytically, we find in general that the inclusion of dispersal serves only to shift the bulk of the eigenvalues in the negative real direction, never in the positive. This is seen to be the case in Fig.  2 in the main text. This can also be seen to be true in the special cases we have discussed already -increasing the magnitude of d α (which is effective what one does by introducing non-zero q) in Eqs. (S3)-(S6) serves only to expand the range of parameters for which the equilibrium is stable.

Turing-type instability due to outlier eigenvalues
The Turing instability is a 'stationary' instability [56] which exhibits no oscillations. This means that the eigenvalue crossing the imaginary axis has no imaginary part (unlike the example given in Fig. S1). The point in parameter space at which an instability of this type occurs in the non-spatial case (q = 0) is given by the solution of Eq. (S7). Using this, and the replacement d α → d α + q 2 D α to account for q = 0, we arrive at the condition for stability with respect to disturbances of non-zero wavenumber q: In order for the ecosystem to show a Turing instability this condition must be satisfied for q = 0 (i.e., the system is stable with respect to perturbations with q = 0), but violated for some non-zero q (there are unstable modes with q = 0). That this can indeed be possible is exemplified by the special case where Γ αβ = 0 for all α and β. In this case, the stability criterion for the Fourier mode q becomes The expression for P 0 (q 2 ) is quadratic in q 2 . It has a minimum at If P 0 (0) > 0 and P 0 (q 2 min ) < 0, then there is a Turing instability. This occurs when where N Cγ u µ uu = N Cγ u µ uu − d u and N Cγ v µ vv = N Cγ v µ vv − d v . This is very reminiscent of the criterion for Turing instability in more conventional reaction-diffusion systems [27].
In the more general case, we need to find the minimum possible value of the expression on the left-hand side of the first relation in Eq. (S15). If the minimum value of P (q 2 ) is negative (i.e. the condition in the first line of Eq. (S15) is violated), but the value of P (0) is positive, then we say that there is a Turing (dispersal-induced) instability for this parameter set.
Again, we first find the minimum value of P (q 2 ) for a given set of model parameters. To do this, we differentiate P (q 2 ) with respect to q 2 , noting that the three relations in Eq. (S15) are coupled. This leads to the following set of simultaneous equations to find P (q 2 min ) where χ α = ∂χα ∂q 2 . One solves the above equations simultaneously to find χ u , χ u , χ v , χ v and q 2 min . One then evaluates So, one can find the minimum value of P (q 2 ) for a particular parameter set and thus deduce whether or not there is dispersal-induced instability (i.e. whether the condition in Eq. (S15) is satisfied). For the purpose of constructing stability diagrams, one can instead set P (q 2 min ) = 0, equivalent to assuming that the outlier eigenvalue is just about to cross the imaginary axis (λ = 0). One then solves for the critical model parameters instead which give P (q 2 min ) = 0. This can be used to produce solid sloped lines in Figs. 3 (b) and 4 in the main text .
Although the above set of equations appears cumbersome, we note that they are at most linear in χ u and χ v , so these can be solved for relatively easily. In practice, Eqs. (S19) reduce to three simultaneous equations for q 2 min , χ u and χ v . C. Dispersal cannot make an unstable equilibrium stable In our model there are no instances in which an unstable equilibrium of the non-spatial ecosystem becomes stable when dispersal is introduced. This is because a stable equilibrium in the spatial model must be stable with respect to disturbances for values of q, including q = 0. This means that such an equilibrium is also stable in the non-spatial system.