Spatial heterogeneity enhance robustness of large multi-species ecosystems.

Understanding ecosystem stability and functioning is a long-standing goal in theoretical ecology, with one of the main tools being dynamical modelling of species abundances. With the help of spatially unresolved (well-mixed) population models and equilibrium dynamics, limits to stability and regions of various ecosystem robustness have been extensively mapped in terms of diversity (number of species), types of interactions, interaction strengths, varying interaction networks (for example plant-pollinator, food-web) and varying structures of these networks. Although many insights have been gained, the impact of spatial extension is not included in this body of knowledge. Recent studies of spatially explicit modelling on the other hand have shown that stability limits can be crossed and diversity increased for systems with spatial heterogeneity in species interactions and/or chaotic dynamics. Here we show that such crossing and diversity increase can appear under less strict conditions. We find that the mere possibility of varying species abundances at different spatial locations make possible the preservation or increase in diversity across previous boundaries thought to mark catastrophic transitions. In addition, we introduce and make explicit a multitude of different dynamics a spatially extended complex system can use to stabilise. This expanded stabilising repertoire of dynamics is largest at intermediate levels of dispersal. Thus we find that spatially extended systems with intermediate dispersal are more robust, in general have higher diversity and can stabilise beyond previous stability boundaries, in contrast to well-mixed systems.

Introduction

There are many ways to think about and represent ecosystem functioning and stability, one prominent direction in mathematical ecology is the search for limits of stability in dynamical species population models in terms of some specified parameters or characteristics of the system. To this end Generalised-Lotka-Volterra (GLV) dynamics and modifications, with for example higher order interaction [1] (e.g. third species modifying how one species interacts with another) as well as many examples of more advanced interaction functions between species have been used. With a proliferation of stability results both confusing and enlightening, such as the role of diversity for stabilising ecosystems [2-4] or to what extent properties such as modularity [5] or nestedness [6] in interactions act to stabilise or destabilise or are aspects of the same property [7].

Specifically the question whether diversity (used as synonymous with species richness) and/or strength of interaction between species act to destabilise or not, has drawn a lot of attention since the first mathematical formulations of large ecosystem stability with May’s paper in the 1970s [2]. The conclusions are not unanimous but a majority of studies find that at a certain level of diversity (or interaction strengths) the system will no longer have a stable equilibrium solution [8]. This means that ecosystems cannot sustain too high diversity, variation or mean of interaction strengths and can possibly collapse when interactions are increased if for example an ecosystem’s available space is reduced. Studies have also located limits to feasibility (abundance >0 for all species) where higher interaction strengths lead to extinctions (loss of feasibility) [9, 10], and systems vulnerable to perturbations in structural features (such as growth rates, or interaction strengths) [11, 12]. With even higher interaction strengths, abundances can start oscillating or the system collapses to a significantly smaller subset of species [9]. This latter limit is the classical limit recognised as a drastic change in system characteristics and many times referred to as collapse.

A prominent model indicating diversity as destabilising is the original Generalised-Lotka-Volterra model. The GLV equations have again become widely used for studies of large ecological communities to investigate generic properties of complex ecosystems using relatively few parameters [13, 14]. Although many insights have been gained by this approach, one of the many simplifications of the GLV model is the lack of space. This simplification leans on the assumption that a spatial average of both the interactions among species and the species abundances is sufficient to capture the dynamics and stability aspects of an ecosystem.

There are on the other hand studies specifically taking space into account, studying for example pattern formation where the diversity is small [15-17], space influence on interaction structure [18] and the spread of specific species [19]. Another type of spatial models are meta-communities. A meta-community is a collection of ecological communities in space connected by dispersal. Most meta-community studies find that intermediate dispersal is favourable for high diversity and robustness [20, 21], although systems under investigation comprised of very few species. It is also known that high dispersal can synchronise fluctuations in population abundances in different habitats, increasing the likelihood of extinctions [20-23]. The mechanism behind both these results is the presence or absence respectively of a buffering effect between the local communities. Other meta-community studies have found meta-communities exhibiting macro-ecological patterns at the edge of collapse [24] and investigated local stability in the feasible domain [25]. A feature that still holds in the latter spatial stability analysis is that there is a boundary of radical shift and loss of stability.

Recent studies using the GLV equations with dispersal have found that spatial extension enables crossing of stability boundaries and higher diversity in comparison to spatially unresolved systems. These features rely on chaotic dynamics [26, 27] and spatial heterogeneity in interactions [26]. Using a spatial gradient representative of abiotic factors, [28] has shown that the level of dispersal influences the abruptness in species composition between neighbouring communities.

Motivated by the large literature on stability limits, structural stability (robustness) and recent interest in spatial extension, in this paper we study a GLV model with diffusion on a lattice with no differentiation between spatial locations such as spatial heterogeneity in interactions, carrying capacities, or intrinsic growth rates. In contrast to spatially unresolved and spatially differentiated models, this formulation lifts the assumption of average species abundances and allows us to investigate the bare consequences of spatial extension. We re-investigate stability limits, robustness and diversity for varying amounts of dispersal.

Theory and methods

The version of the classical GLV model we will use as a base is stated below where ϕ, r and K are species abundances, intrinsic growth rates and carrying capacities for species i respectively. The web of interactions between species is represented by A, a N × N matrix with a density of interactions, c = 0.5, and whose non-zero entries are drawn from a normal distribution with negative mean, μ = −0.5 and variance one. This means all types of interactions (trophic, mutualistic, competitive etc.) are included although there is a larger percentage of competitive (−, −) and amensialistic (−, 0) interactions. The parameter σ becomes the standard deviation (s.d.) of the interaction strengths and is often used as a tuning parameter and proxy for complexity. An increase in interaction strengths means an increase in complexity.

For small σ in systems governed by Eq 1 with diversity N the system will settle in a unique stable fixed point with all N species present. Under increases in σ (increasing complexity) at a certain threshold feasibility is lost. With continued increase beyond this threshold to stay in a stable fixed point species will successively go extinct. Further increase in σ eventually pushes the system across the final stability boundary and the system transitions to either oscillations, chaos or a fixed point with a substantial loss of species. This latter limit has many times been referred to as collapse. The region between the two boundaries is structurally unstable meaning a small perturbation in parameters (r, σ, c etc.) lead to qualitative change, in effect species extinctions [9]. This region can also have multiple stable fixed points with differing patterns of extant (non-extinct) species [29].

We introduce a spatial dimension into the GLV model by setting up a grid (or line) where all grid-points have the same interaction matrix. This in effect, sets the same maximum amount of species with the same interspecies interactions at every grid-point but allows for the possibility of different dynamics and local population abundances. The grid-points are connected by diffusion, representing the dispersal of species in space. This can be seen as for example animals colonising areas where there is a larger abundance of prey per predator or plants dispersing to areas where there are less competitors (more available resources).

To avoid effects from the boundary of the grid we use periodic boundary conditions (shape the grid as a torus in two-dimensional space or ring in one-dimensional space) as shown in Fig 1.

Fig 1

Schematic grid and line.

The figure shows schematic pictures of how the grid is formed like a torus (to the right) and how the line is formed like a ring (to the left). This is done to simulate an ecosystem situated in a large connected space and to minimise boundary effects.

The GLV equations in continuous space with diffusion are where ϕ(x, t) is the relative abundance of species i which now depends on both spatial location x and time t. The diffusion rates for species i are D, which are different for all species but same for all grid-points for the same species. The D are random variables drawn from uniform distributions with means μ = 10 with a, a parameter we varied from -5 to 1 in steps of 1. The different D for different species is to represent that species have varying dispersal rates, although all the D are kept at the same order of magnitude for each run by a standard deviation proportional to the mean σ = bμ, with 0.3 < b < 0.4, such that the diffusion intervals for different means are non-overlapping. Since we wish to use a grid we discretise the spatial dimension of Eq 2 using the a discrete Laplace operator in two or one dimension given as where α and β are grid indices and the denominator h is the “distance” between patches set to 1. The resulting dynamical equation for species i in grid-point (α, β) in a two dimensional grid is thus

The dispersal in our model is thus only between nearest neighbour grid-points.

Real ecosystems are often situated in either two- or three-dimensional space. In most of our simulations on the other hand we used one-dimensional space to reduce the computational load and effect of boundaries. This was acceptable since we noticed no qualitative difference in results between two- and one-dimensional simulations.

To analyse the dynamics of the model, we screen each grid-point for fixed points and oscillations in species abundances. If there are oscillations we measure amplitudes and the synchronisation of local species abundance oscillations between grid-points. The synchronisation of dynamics we measure with the maximum phase-shift in oscillations between the same species in different grid-points, in effect the maximum difference between arguments of the Fourier transform at the dominant oscillation frequency. When the mean maximum phase shift of all extant species is zero the grid is synchronised, when π radians (180°) completely unsynchronised.

Results

With the addition of a connected space we add the possibility of different dynamics or solutions, and timing of dynamics in different grid-points. Spatial heterogeneity (local species abundance differences) because of different dynamics is only possible when σ is large enough (larger than the feasibility boundary) so that different solutions of the equations are available in the grid-points. In this case the system can end up with multiple fixed points, oscillatory patterns, or a combination of fixed points and oscillations. The same oscillatory dynamics can also give rise to spatial heterogeneity by phase-shifts in local species abundances between grid-points. The amount of spatial heterogeneity in dynamics and abundances depends on the magnitude of the diffusion constants D as shown in Fig 2. If diffusion constants are zero this represents a number of identical disconnected systems with GLV dynamics, while very large constants lead to complete synchronisation of the grid, smoothing out all local differences in the abundances of species. All these local variations at lower diffusion magnitudes lead to a multitude of system characteristics in terms of combinations of grid-point dynamics and phase-shifts.

Fig 2

Spatial structure for varying diffusion rates.

The figure shows three columns of local species abundances at a specific time for each grid-point from three example runs. Diffusion rates (uniformly distributed with μ ∝ σ) increase from top to bottom. The trend toward synchronisation of the grid with increasing diffusion magnitudes is clearly seen, although the three systems are seen to synchronise at different magnitudes.

Abundance oscillations and crossing stability limits

Species abundance oscillations (global or local) can appear in the entire non-feasible region. We find that in many systems for the spatial GLV, increasing σ can make a fixed point transition to a oscillatory pattern with the same diversity. This behaviour is not as prevalent in the non-spatial GLV, which instead commonly switches to another fixed point with a lower species richness, in effect species go extinct if the system is perturbed or pressured. The connected space enables the system to keep its species richness unchanged if perturbed with either global species abundances (species abundances for the meta-populations) oscillating if the grid-points are synchronised or “stay” in a global fixed point by local unsynchronised abundance oscillations. The spatial systems are therefore more robust in terms of extinctions and in the latter case also in global abundances.

If diffusion constants are small so that synchronisation does not take place, and systems are at interaction strengths nearing previous collapse boundaries this means averaging over local species abundances results in global average abundances in accordance with fixed points that are unstable in the non-spatial GLV. Thus, this robustness of the spatial GLV makes it possible for ecosystems with internal heterogeneity to stay stable in parameter regions beyond stability limits of the non-spatial GLV.

In Fig 3 is an example of a system of N = 20 maximum biodiversity for a range of σ larger than the feasibility limit, and diffusion magnitudes leading to unsynchronised oscillations (μ = 10−2, σ = 0.4μ). We see that the unsynchronised local species abundance oscillations lead to stable global species abundances at an interaction strength where there is no fixed point in the non-spatial GLV. Thus the previous stability boundary is crossed by the spatial GLV without chaotic dynamics and almost without any change of the global species abundances.

Fig 3

Unsynchronised local oscillations stabilise.

The figure shows an example of a system of N = 20 maximum number of species with out-of-phase oscillations when σ is larger than the feasibility limit and diffusion rates are small (μ = 10−2, σ = 0.4μ). On the left side on top we show global average fixed point or oscillation species abundances for a system with small diffusion. On the bottom we show fixed point abundances for a system with the same interaction matrix but with no spatial dimension (non-spatial GLV). Both the global average abundances for the spatial system and the non-spatial are shown for σ ranging from zero to collapse values with an enlargement of the latter part in the red box. The Green shading in the plots indicate a region where the only stable fixed point for the non-spatial GLV is with 3 species going extinct. On the other hand the global average abundances show no change at all (green area in top left plots). With higher σ in the orange region of the non-spatial GLV (bottom left plot) the system is seen to be structurally unstable, while the spatial system on top shows little if any structural instability. These systems behave almost the same with and without a connected space until σ is large enough approaching collapse values. To the right on the green background are shown example dynamics for σ in the green marked area in the left plots, in different grid-points for the spatial system (top) and for a oscillatory solution for the non-spatial system (bottom). We see in the spatial system oscillatory dynamics in each grid-point example with the same Fourier spectra, but differing phases (the panels to the right). Together the different phases and amplitudes but same frequencies of the local abundance oscillations average to the values corresponding to an unstable fixed point of the non-spatial GLV. For the non-spatial system there is a oscillatory pattern, note however the increase in sharpness in both frequency and amplitude as well as some species going extinct and reappearing, which is not a biologically realistic or stable solution for a ecological system.

Lower variability in abundances

There are different ways the inclusion of a connected space acts to stabilise global and local species abundance oscillations making them less variable (lower amplitudes and/or lower frequencies). The oscillations in the spatial GLV usually reflect a pattern found in the non-spatial GLV, but with the addition of diffusion local abundance oscillations become less violent with lower amplitudes.

Other systems with oscillations might have different “preferred” oscillations for the non-spatial and spatial GLV respectively. In these cases the non-spatial oscillation patterns are higher in amplitudes and sometimes with higher frequencies than in the spatial GLV, see for example the non-spatial GLV oscillations in Fig 3 for one version of such an oscillation pattern. Yet another possibility is that diffusion makes possible a different lower amplitude oscillation pattern not present in the non-spatial GLV, an example of this is shown in Fig 4.

Fig 4

Local and global species abundances in spatial extended system.

The figure shows an example of a system with varying rates of diffusion from top to bottom μ = 2.5σ = 10−5, 10−3, 10−1, 1 (the green shading increases with increasing diffusion rates). The left column shows the system’s dynamics at arbitrary grid-points from high amplitude with low diffusion as the amplitudes are dampened with increasing diffusion and then back to high amplitudes as the diffusion rate is high enough for synchronisation over the grid. If there are two panels in the left column this shows two unsynchronised grid-point dynamics. If there is only one panel the system is synchronised. The right column shows the spatial mean for all grid-points (global average abundances). Note that the unsynchronised oscillations lead to almost constant global species abundances while the synchronised system on the bottom shows violent global oscillations.

Many of these violent oscillations in the non-spatial system lead to periodic species extinctions, an unrealistic scenario. The spatial GLV on the other hand can avoid global extinctions both by less violent local oscillations leading to less local extinctions as well as recolonisation-extinction dynamics between local areas. Significant lowering of the variability in local species abundances we find when systems are unsynchronised at low to medium diffusion rates. Although such systems commonly have a large spread in local oscillation amplitudes, all grid-points have dampened oscillation amplitudes compared to the non-spatial system and the global abundances are almost constant. When the diffusion rates in systems like these are increased, engendering synchronised oscillations, they usually reflect oscillation patterns present in the non-spatial GLV but again with reduced amplitudes.

In the limit of large diffusion (μ > 1 and σ = bμ, 0.3 < b < 0.4) the ecosystems again retrieve the high amplitudes of the non-spatial GLV oscillations synchronised in all grid-points. Thus when diffusion is increased first oscillations start to synchronise leading to global abundance oscillations and with continued increase in diffusion the species oscillations increase in amplitude until they mimic that of the non-spatial GLV. It is therefore in the middle range of diffusion rates we find the most stabilising effects as seen in Fig 4.

Parameter-space

We have shown examples of dynamics of the spatial GLV and also pointed out the multitude of ecosystem dynamics in a connected space. We can bring some order to this multitude and substantiate our claim that moderate diffusion rates lead to increased stability by gathering statistics on measures that contribute to stability, while varying σ and the diffusion rate. The measures we chose for each run are A) Number of grid-points with oscillatory patterns (potentially different patterns), B) Average maximum phase shift for a species between grid-points if oscillations (degree of synchronisation), C) Average amplitude if oscillations, and D) Average diversity (n). The measures are both in combination and one-by-one the basis for different aspects of stability. Statistics over these measures in the non-feasible region for varying magnitude of diffusion rates are shown in Fig 5.

Fig 5

Parameter-space—Diffusion magnitude vs. s.d. of interaction strength σ.

The figure shows statistics for 70 systems in diffusion magnitude vs. standard deviation parameter-space for panel A) Number of grid-points with oscillatory patterns, panel B) Average maximum phase shift for a species between grid-points if oscillations (degree of synchronisation), panel C) Average amplitude if oscillations, and panel D) Average diversity (n). The lowest diffusion in the diagrams is zero, corresponding to completely disconnected space (non-spatial GLV) and the largest D ∼ 1 (μ = 1, σ = bμ, 0.3 < b < 0.4). Worthy of noting is the lower degree of oscillations and the lower diversity in the non-spatial systems. It is also clear that oscillations are present in almost the entire non-feasible region.

An example of the measures contributing to stability is, the oscillations appearing in almost the entire non-feasible region, although most prominently for higher interaction strengths (panel A), while the average diversity for systems with diffusion is higher than without diffusion (panel D). This shows the increase in robustness across the entire span of interaction strengths, σ. Systems can react with local oscillations instead of extinctions if pressured or perturbed. Another example, the lesser fraction of stable fixed points for large σ together with less oscillations at D = 0 for all i (bottom right panel A), show where systems without a connected space start collapsing. Adding moderate rates of diffusion so that the system has a connected space gives rise to unsynchronised oscillations and the ability for systems to stay stable across previous stability limits. Panel C shows the decrease in local species abundance variability (oscillation amplitudes) for the middle range of diffusion magnitudes, where both disconnected and high diffusion lead to violent oscillations. This feature is connected to the degree of synchronisation shown in panel B, where we can see the increase in synchronisation as the diffusion constants are increased in magnitude.

In Fig 6 we summarise the results from the parameter-space in terms of behaviour and stability. The table shows that the middle range of diffusion, in effect systems which allow for a certain amount of dispersal of species between local areas, are the more robust systems over a larger range of interaction strengths.

Fig 6

Stability of spatial GLV.

The table shows in text the type of dynamics present for low, medium and high diffusion rates and s.d of interaction strength. Green indicates stable robust systems and the darker the colour the more robust. Red on the other hand are collapsed or unstable volatile systems. Note that the middle range is green for all s.d of interaction strengths.

The diffusion magnitude, as we have shown, we found to have a large impact on dynamics, the variation of diffusion rates between species however had little impact. Different random draws of species diffusion constants with the same magnitude could phase-shift oscillations between runs and introduced a very slight nudge towards asynchronisation, but otherwise introduced no qualitative change in oscillatory patterns or fixed point solutions in the grid-points. The slight nudging towards unsynchronised dynamics could point towards a stabilising effect from large variation in dispersal rates between species, but this was not further investigated.

Discussion

Ever since May’s mathematical argument for the instability of large complex systems [2] challenging the previous view of ‘complexity begets stability’ advocated by MacArthur [30], it has been an open question if the biodiversity, multitude of interactions between species and the amount of interactions (in effect complexity) of an ecosystem makes it more or less stable. Since then a limit to the complexity an ecosystem can sustain, derived in May’s paper and later extended, has persisted. Ecosystems at the edge of this limit or externally pushed so that their complexity exceeds it are predicted to radically change, lose species or collapse. Crucial aspects of such models are that the ecosystems are modelled as isolated homogeneous systems. A resolution to the long standing question of the limit of complexity might not be primarily the structure of interactions but that real ecosystems are both internally heterogeneous and externally connected. Connected ecosystems might gain their capacity of sustaining previously thought to be unstable patterns of cohabitation using the flexibility of local variations and buffer of surrounding ecosystems.

Our model represents an ecosystem or high diversity meta-community for which the abiotic conditions are constant, reflected by our choice of identical interaction strengths, carrying capacities and intrinsic growth rates throughout the connected space. Large ecosystems might include many different types of habitats connected by dispersal, which is mostly found to be positive for biodiversity and robustness [25, 26, 31]. But even in ecosystems with only one type of habitat, the area could be large enough such that external perturbations are not equal, there might be boundaries artificial or natural, or aggregation of species due to conspecific attraction [32, 33], leading to patches with equal species interactions connected by dispersal. As we have shown, in such ecosystems the mere possibility of differential species abundances at different locations enhances stability and diversity.

The classical stability boundaries have long been known to be a transition from a stable equilibrium fixed point to some qualitatively different dynamics in terms of abundance oscillations, chaos or drop to substantially smaller community [2, 9, 34]. Thus oscillations have been known to be a possible dynamics resulting from crossing the stability boundary. In a connected space the amplitudes of such oscillations are dampened. In addition, if the dispersal rates between patches are low enough, oscillations in local species abundances can be out of phase leading to global abundances remaining largely unchanged. This result is a spatial extension of the time-average effect, which is that time-averages of fluctuating species abundances in an ecosystem are nearly constant [35, 36]. We can then speculate that the “Portfolio effect”, keeping ecosystem properties by unsynchronised oscillations of functionally equal species, is an additional aspect of robustness gained with spatial extension.

The effect of local stabilising oscillations is as we have shown not only present when nearing collapse, but can also prevent extinctions when the system is perturbed. The non-spatial GLV is structurally unstable in the non-feasible region, meaning small structural changes (small increase in σ) can lead to species extinctions. The connected space, on the other hand can stabilise by local oscillations, thus avoiding the irreversible event of an extinction. With real ecosystems more alike the spatial GLV, we can expect ecosystem to be more robust than many models might suggest [37, 38].

One of the reported mechanisms for high diversity in spatial GLV models with spatial heterogeneity in interactions is that some local habitats function as sources thereby preventing some species from extinction [26]. In our study, this is not a possible mechanism for high diversity, instead we find that spatial connectivity both promotes oscillatory patterns with higher diversity than fixed points as well as allowing for a larger spectrum of the possible dynamics some of which suffer less extinction and thus contribute to higher global diversity. A similar mechanism is found in chaotic dynamics, where it is reported that booms and bust (asynchronous abundances) at different identical locations nurture persistent high diversity dynamics [27].

Adding a connected space increases the dynamics available for a system. This is apparent in our study as well as in [28], studying the abruptness between such different dynamical realisations along a spatial gradient. We find that not only is there room for different combinations of dynamics in the grid but the connected space can also facilitate dynamics not found or unstable in spatially unresolved systems. This is a source of stability, but we also argue that the increased dynamical possibilities is itself a mechanism of robustness for the system. Importantly we also find that this repertoire of dynamics is largest at intermediate levels of dispersal, an observation in line with previous studies [21, 25, 26, 39].

It is quite remarkable that the inclusion of space even with homogeneous patches and without external forcing is sufficient for dynamic spatial heterogeneity to emerge in the GLV, including regimes of stable global abundances and persistent local abundance fluctuations. High diversity meta-communities have been shown to also display local species turn-over in absence of external influence [40]. The fact that undisturbed high diversity meta-communities have such a wealth of dynamical behaviour can have a large impact on biological monitoring. To correctly asses and manage ecosystems it is vital that we can predict what behaviour can be expected from intrinsic system dynamics so that we do not wrongfully assign such dynamics to environmental change or anthropogenic pressures and invest our efforts where it is uncalled for.

With the help of spatially extended models many novel mechanisms promoting diversity and stability have been found and intrinsic dynamics uncovered, but there are avenues relevant to real ecosystems still unexplored. In this study we used different diffusion constants for different species but all within the same order of magnitude, with no definite qualitative change in stability properties compared to constant diffusion rates. It is still unclear how a larger spread in species mobility, larger spread in diffusion constants, would effect stability. This, both because of the very slight nudge towards unsynchronised dynamics and since spatial patterning in reaction diffusion systems tend to appear in systems where the components have a large spread in diffusion constants. We can speculate that spatial species patterning might appear in such systems. Varying the diffusion constants between grid-points is also a direction left unexplored. Varying constants can represent varying distance, some natural obstacle, a road or systems with internal dynamics connected to surrounding ecosystems. With varying diffusion constants we can therefore both study “unsymmetrical” habitats where interactions between local areas differ as well as how local disturbances might affect an ecosystem and connected ecosystems.

Left for future work is also varying the standard deviation of interaction strength σ over the grid. The inverse interaction strength can be interpreted as a proxy for available habitat, since a reduction of available space will force the existing species closer together and to interact more which is an increase in σ. Thus, with different σ we can model a grid with different local habitat sizes as well as a system’s response to local habitat losses if perturbing one or some of the σ at a time.

Conclusion

A large ecosystem can have internal spatial heterogeneity in local species abundances. Although, many insights into ecosystems’ functioning and stability have been gained through studies where such heterogeneity is averaged out, in this study we find that adding this possibility promotes ecosystem robustness and diversity. We add the possibility of abundance heterogeneity by modelling an ecosystem as collection of communities connected by dispersal. In this spatially resolved system a limit of complexity for stability is no longer a limit to qualitative change or collapse. With the connected space acting as a buffer harbouring unsynchronised local abundances, global abundances are kept constant beyond stability limits of spatially homogeneous systems. We also found that spatially extended systems have lower variability in abundance fluctuations and the ability to avoid extinctions by local species oscillations, thus promoting high diversity. The increase in dynamical possibilities and combinations is a source of robustness for the global system. This repertoire of dynamics is maximal at intermediate dispersal which we therefore find to promote the most robust ecosystems.

1 Jun 2021

Dear Miss Pettersson,

Thank you very much for submitting your manuscript "Spatial heterogeneity enhance robustness of large multi-species ecosystems" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Thank-you for your submission and apologies for the delay in providing a recommendation.

Both reviewers see the question of how spatial mechanisms influence both stability and robustness as an extremely important area of investigation. However, both reviewers also question whether there is a sufficient amount of novel insight (relative to existing literature) in the current manuscript to merit publication. Reviewer 1 highlights specific recent papers combining pairwise species interactions and spatial processes---I agree with the reviewer that it is important to highlight what precisely are the differences and novel findings here so that readers will understand what is being added to our understanding. Reviewer 2 raises more general concerns about the broader context of diversity-stability relationships. Like the reviewer, I am not sure that the goals here are identical to those papers (and perhaps they speak more to the question of robustness in the current ms than to local stability) but it is fair to say that readers may find this context extremely helpful.

In seeing these concerns, I believe that revisions may be difficult. I do think that it would be possible to send the paper out to review again if the authors were to choose to resubmit to PLOS Comp Biol, but the authors would need to highlight more clearly the comparison with existing results and broader concepts.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Natalia L. Komarova

Deputy Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************

Thank-you for your submission and apologies for the delay in providing a recommendation.

Both reviewers see the question of how spatial mechanisms influence both stability and robustness as an extremely important area of investigation. However, both reviewers also question whether there is a sufficient amount of novel insight (relative to existing literature) in the current manuscript to merit publication. Reviewer 1 highlights specific recent papers combining pairwise species interactions and spatial processes---I agree with the reviewer that it is important to highlight what precisely are the differences and novel findings here so that readers will understand what is being added to our understanding. Reviewer 2 raises more general concerns about the broader context of diversity-stability relationships. Like the reviewer, I am not sure that the goals here are identical to those papers (and perhaps they speak more to the question of robustness in the current ms than to local stability) but it is fair to say that readers may find this context extremely helpful.

In seeing these concerns, I am suggesting a rejection. I do think that it would be possible to send the paper out to review again if the authors were to choose to resubmit to PLOS Comp Biol, but the authors would need to highlight more clearly the comparison with existing results and broader concepts.

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors look at dynamics of interacting species (via the GLV equations), and add connected space to the system. The species interact in each location and are able to migrate between the different locations. They find that stability limits on GLV system when considered without space, are no longer restrictive due to the possibility of spatial heterogeneity, and the related appearance of unsynchronized oscillations of the abundances.

The manuscript looks well-written and draws clear and potentially interesting conclusions. But at the moment it is hard for me to judge the novelty of the work because it seems, as far as I can see, very closely related to the following publications:

* Roy, Felix, Matthieu Barbier, Giulio Biroli, and Guy Bunin. “Complex Interactions Can Create Persistent Fluctuations in High-Diversity Ecosystems.” PLOS Computational Biology 16, no. 5 (2020): e1007827.

* Pearce, Michael T., Atish Agarwala, and Daniel S. Fisher. “Stabilization of Extensive Fine-Scale Diversity by Ecologically Driven Spatiotemporal Chaos.” Proceedings of the National Academy of Sciences 117, no. 25 (2020): 14572–83.

In both these references, different locations in space with GLV (or replicator equations) are coupled via migration. Dynamics of the abundances appear and are argued to allow for stability bounds to be crossed. The theme of de-synchronization plays a role in these references as well, as does the non-trivial dependence on migration rates.

One difference between the present manuscript and these references, is that the locations in space are here placed on a grid versus all-to-all connections. But in the present manuscript the topology of the connections is claimed not to be very important to the qualitative results (when using 1d versus 2d or 3d grid).

It may very well be that the present manuscript contains novel results beyond these publications, but in order to make a judgement, I'd first like to see a comparison by the authors, between their work and the above references.

Reviewer #2: Please find my comments attached

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: No: I see no working links to data or code repositories

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: Reply to Spatial heterogeneity enhance robustness of large multi-species ecosystems.odt

Click here for additional data file.

4 Aug 2021

Submitted filename: ResponseToReviewers.docx

Click here for additional data file.

25 Aug 2021

Dear Miss Pettersson,

Thank you very much for submitting your manuscript "Spatial heterogeneity enhance robustness of large multi-species ecosystems" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Thanks for your work on addressing the reviewer comments. I would be happy to consider a revised manuscript that addresses the remaining (relatively minor) comments of reviewer 2

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

James O'Dwyer

Deputy Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Thanks for your work on addressing the reviewer comments. I would be happy to consider a revised manuscript that addresses the remaining (relatively minor) comments of reviewer 2

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have taken the comments of both myself and the other reviewer into account. The new version of the manuscript better connects with the existing literature on spatial effects, and in particular metacommunities.

The manuscript contains a number of interesting observations, specifically on the conditions under which space can stabilize communities at higher species richness, and the dynamics which allow this to happen.

Reviewer #2: Review uploaded as attachment

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: None

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Jacob O'Sullivan

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

Submitted filename: Reply to Spatial heterogeneity enhance robustness of large multi-species ecosystems resubmission.odt

Click here for additional data file.

1 Oct 2021

Submitted filename: ReplyToReviewers.docx

Click here for additional data file.

7 Oct 2021

Dear Miss Pettersson,

We are pleased to inform you that your manuscript 'Spatial heterogeneity enhance robustness of large multi-species ecosystems' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

James O'Dwyer

Deputy Editor

PLOS Computational Biology

Natalia Komarova

Deputy Editor

PLOS Computational Biology

***********************************************************

21 Oct 2021

PCOMPBIOL-D-21-00445R2

Spatial heterogeneity enhance robustness of large multi-species ecosystems

Dear Dr Pettersson,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Olena Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol