Far-ranging generalist top predators enhance the stability of meta-foodwebs

TL;DR

This study demonstrates that far-ranging generalist top predators significantly enhance the stability of complex meta-foodwebs, emphasizing their crucial role in ecosystem conservation and stability.

Contribution

It introduces a novel analysis showing that far-ranging top predators increase stability in large, realistic meta-foodwebs using generalized modeling and stability functions.

Findings

Meta-foodwebs with far-ranging predators are more stable.

Adding far-ranging predators stabilizes large foodwebs.

Far-ranging predators have a net stabilizing effect despite increasing web size.

Abstract

Identifying stabilizing factors in foodwebs is a long standing challenge with wide implications for community ecology and conservation. Here, we investigate the stability of spatially resolved meta-foodwebs with far-ranging super-predators for whom the whole meta-foodwebs appears to be a single habitat. By using a combination of generalised modeling with a master stability function approach, we are able to efficiently explore the asymptotic stability of large classes of realistic many-patch meta-foodwebs. We show that meta-foodwebs with far-ranging top predators are more stable than those with localized top predators. Moreover, adding far-ranging generalist top predators to a system can have a net stabilizing effect, despite increasing the food web size. These results highlight the importance of top predator conservation.

Figures

Tables

Table 1. Table 1: Generalized parameters used to describe the meta-foodweb.

Parameter	Interpretation
Exponent
$\varphi$	Sensitivity of primary production to own population density
$\gamma$	Sensitivity of predation to total available prey density
$\lambda$	Exponent of prey switching
$\psi$	Sensitivity of predation to predator density
$\mu$	Exponent of closure
$\omega$	Sensitivity of dispersal to own the species population density in the output patch
$\tilde{\omega }$	Sensitivity of dispersal to own species population density in the target patch
$\kappa$	Sensitivity of dispersal to population density of some other species in the output patch
$\hat{\kappa }$	Sensitivity of dispersal the population density of some other species in the target patch
Scale
$\alpha$	Biomass flow
$\sigma$	Fraction of biomass loss due to predation
$\tilde{\sigma }$	Fraction of biomass loss due to respiration
$\xi$	Fraction of loss by predation due to global species
$\tilde{\xi }$	Fraction of loss by predation due to local species
$\beta$	Relative contribution to biomass loss due to predation by a certain predator
$\delta$	Fraction of local growth by predation
$\tilde{\delta }$	Fraction of local growth by primary production
$\chi$	Relative contribution of population as prey to a certain predator
$\nu$	Fraction of total biomass gain due to dispersal
$\tilde{\nu }$	Fraction of total biomass gain due to predation and primary production
$\rho$	Fraction of total biomass loss due to dispersal
$\tilde{\rho }$	Fraction of total biomass loss due to predation and respiration

Table 2. Table 2: Constant values and distribution intervals of the exponent parameters. In the case of the ’Adaptive’ model the parameters κ X subscript 𝜅 𝑋 \kappa_{X} and κ ^ X subscript ^ 𝜅 𝑋 \hat{\kappa}_{X} depend on the feeding links of the local foodweb. If species i 𝑖 i is the prey of species j 𝑗 j then κ X i ​ X j subscript 𝜅 subscript 𝑋 𝑖 subscript 𝑋 𝑗 \kappa_{X_{i}X_{j}} is randomly assigned a value inside the interval [ 0 , 1 ] 0 1 [0,1] and is otherwise set to 0. This means the dispersal of prey is stimulated if predators are present. Additionally the parameter κ ^ X j ​ X i subscript ^ 𝜅 subscript 𝑋 𝑗 subscript 𝑋 𝑖 \hat{\kappa}_{X_{j}X_{i}} is randomly assigned a value inside the interval [0,1] if species j 𝑗 j does feed on species i 𝑖 i and is set to 0 otherwise. This corresponds to predators that prefer to disperse into patches with higher prey abundance.

	${\varphi }_X$ / ${\varphi }_Y$	${\psi }_X$ / ${\psi }_Y$	${\mu }_X$ / ${\mu }_Y$	${\gamma }_X$ / ${\gamma }_Y$	${\lambda }_X$ / ${\lambda }_{YX}$	${\omega }_X$	${\hat{\omega }}_X$	${\kappa }_X$	${\hat{\kappa }}_X$
Holling T2	1	1	[1,2]	[0,1]	[1,2]	[0,1]	[0,1]	0	0
Standard	[0,1]	[0,1]	[1,2]	[0,2]	[1,2]	[-2,2]	[-2,2]	[-2,2]	[-2,2]
Adaptive	[0,1]	[0,1]	[1,2]	[0,2]	[1,2]	[0,1]	[0,1]	[0,1]	[0,1]

Equations

$${\dot{X}}_i^k$$

$${\dot{Y}}_i$$

$$\mathbf{J}=\left(\begin{array}{cc}{\mathbf{J}}_X& {\mathbf{J}}_{XY}\\ {\mathbf{J}}_{YX}& {\mathbf{J}}_Y\end{array}\right),$$

$${\mathbf{J}}_X=\mathbf{I}\otimes {\mathbf{P}}_X-\mathbf{L}\otimes {\mathbf{C}}_X$$

$$\displaystyle\dot{x}_{i}^{k}=\alpha_{X_{i}^{k}}\Bigg{[}$$

$$\displaystyle+\sum_{l}\nu_{X_{i}^{k}X_{i}^{l}}e_{i}^{kl}(\mathbf{x}^{k},\mathbf{x}^{l})-\sum_{l}\rho_{X_{i}^{l}X_{i}^{k}}e_{i}^{lk}(\mathbf{x}^{l},\mathbf{x}^{k})-\tilde{\rho}_{X_{i}^{k}}\sigma_{X_{i}^{k}}\xi_{X_{i}^{k}}\sum_{j}\beta_{Y_{j}X_{i}^{k}}d_{Y_{j}X_{i}^{k}}(\mathbf{x},\mathbf{y})\Bigg{]}\,,$$

$$\displaystyle\dot{y}_{i}=\alpha_{Y_{i}}\Bigg{[}\tilde{\delta}_{Y_{i}}g_{Y_{i}}(y_{i})+\delta_{Y_{i}}f_{Y_{i}}(t_{Y_{i}},y_{i})-\tilde{\sigma}_{Y_{i}}m_{Y_{i}}(y_{i})-\sigma_{Y_{i}}\sum_{j}\beta_{Y_{j}Y_{i}}d_{Y_{j}Y_{i}}(\mathbf{x},\mathbf{y})\Bigg{]}\,,$$

$${\varphi }_{X_i^k}$$

$${\mu }_{X_i^k}$$

$${\lambda }_{X_j^k{X}_i^k}$$

$${\lambda }_{Y_j{X}_i^k}$$

$${\gamma }_{X_i^k}$$

$${\psi }_{X_i^k}$$

$${\hat{\omega }}_{X_i^k{X}_i^l}$$

$${\hat{\kappa }}_{X_i^k{X}_j^l}$$

$$r=\frac{r^{\prime }}{\sqrt{N}},$$

$${\alpha }_i$$

$${\tilde{\rho }}_i{\alpha }_i={\tilde{\nu }}_i{\alpha }_i=a_i\cdot 10^{-2n_i}$$

$${\rho }_i^{kl}{\alpha }_i={\nu }_i^{kl}{\alpha }_i=(1-a_i)\cdot 10^{-2n_i}.$$

$${\beta }_{ij}$$

$${\xi }_i=\frac{s_{gi}}{{\sum }_n{s}_{ni}},$$

$${\mathbf{J}}_X=\left(\begin{array}{ccccc}\ddots & ⋮& & ⋮& \\ \cdots & {{\mathbf{O}}_X}^{ik}& \cdots & {{\mathbf{O}}_X}^{il}& \cdots \\ & ⋮& & ⋮& \\ \cdots & {{\mathbf{P}}_X}^k-{\sum }_m{{\mathbf{C}}_X}^{km}+{{\mathbf{O}}_X}^{kk}& \cdots & +{\hat{\mathbf{C}}}_X{}^{kl}+{{\mathbf{O}}_X}^{kl}& \cdots \\ & ⋮& \ddots & ⋮& \\ \cdots & {\hat{\mathbf{C}}}_X{}^{lk}+{{\mathbf{O}}_X}^{lk}& \cdots & {{\mathbf{P}}_X}^l-{\sum }_m{{\mathbf{C}}_X}^{lm}+{{\mathbf{O}}_X}^{ll}& \cdots \\ & ⋮& & ⋮& \\ \cdots & {{\mathbf{O}}_X}^{jk}& \cdots & {{\mathbf{O}}_X}^{jl}& \cdots \\ & ⋮& & ⋮& \ddots \end{array}\right).$$

$${\mathbf{J}}_X=\mathbf{I}\otimes {\mathbf{P}}_X^{\prime }-\mathbf{L}\otimes {\mathbf{C}}_X+\mathbf{N}\otimes {\mathbf{O}}_X,$$

$$L_{kj}$$

$$\mathbf{L}$$

$$({\mathbf{P}}_X^{\prime }-\kappa {\mathbf{C}}_X+\mu {\mathbf{O}}_X)\cdot \mathbf{q}={\lambda }_{\kappa \mu }\mathbf{q}.$$

$$v=\mathbf{p}\otimes \mathbf{q},$$

$${\mathbf{J}}_X\cdot \mathbf{v}$$

$$=\mathbf{I}\cdot \mathbf{p}\otimes {\mathbf{P}}_X^{\prime }\cdot \mathbf{q}-\mathbf{L}\cdot \mathbf{p}\otimes {\mathbf{C}}_X\cdot \mathbf{q}+\mathbf{N}\cdot \mathbf{p}\otimes {\mathbf{O}}_X\cdot \mathbf{q}$$

$$=\mathbf{p}\otimes {\mathbf{P}}_X^{\prime }\cdot \mathbf{q}-\kappa \mathbf{v}\otimes {\mathbf{C}}_X\cdot \mathbf{q}+\mu \mathbf{p}\otimes {\mathbf{O}}_X\cdot \mathbf{q}$$

$$=\mathbf{p}\otimes ({\mathbf{P}}_X^{\prime }-\kappa {\mathbf{C}}_X+\mu {\mathbf{O}}_X)\cdot \mathbf{q}$$

$$=\mathbf{p}\otimes {\lambda }_{\kappa \mu }\mathbf{q}={\lambda }_{\kappa \mu }(\mathbf{p}\otimes \mathbf{q})={\lambda }_{\kappa \mu }\mathbf{v}.$$

$$\left(\begin{array}{cc}{\mathbf{J}}_X& {\mathbf{J}}_{XY}\\ {\mathbf{J}}_{YX}& {\mathbf{J}}_Y\end{array}\right)\cdot \left(\begin{array}{c}\mathbf{v}\\ 0\end{array}\right)$$

$${\mathbf{J}}_{YX}\cdot \mathbf{v}\stackrel{!}{=}0.$$

$${\mathbf{J}}_{YX}=(1,\dots ,1)\otimes {\mathbf{C}}_{YX}.$$

$${\mathbf{J}}_{YX}\cdot \mathbf{v}$$

Full text

Introduction

A persistent theme in community ecology is the quest to understand the factors that stabilize foodwebs. Between early notions of “Complexity begets stability”, May’s mathematical argument that more complex foodwebs should be less stable[1], and a number of recent explorations using simulations and stability analysis, a plethora of different hypotheses on foodweb stability have been explored. Among the most prominent of these is the stabilizing effect of weak links [2, 3], although other evidence suggests a potentially destabilizing effect of these links [4, 5]. Further stabilizing effects are ascribed to allometric scaling [6, 7], nontrophic interactions [8], and certain foodweb motifs [9, 10, 11, 5]. The topic of foodweb stability has gained renewed urgency by recent papers raising concerns that past perturbations have not only led to the loss of species but left ecological systems in a significantly more fragile state [12, 13]. A likely candidate contributing to this fragility of previously disturbed systems is the loss of top predators. Habitat loss, overexploitation, and numerous other stress factors have caused global declines in top predators [14]. For example in the British isle all of the large top predators such as the wolf, bear, and lynx have gone extinct.

Among the species in a foodweb top predators appear to have a disproportionate impact on stability. Althoug they constitute typically only a tiny fraction of the biomass [15], top-predators have been found to have a strong impact on stability [16, 5, 17]. Top predators are often generalists which connect and balance biomass from different types of specialist predators on lower trophic levels [18]. Thereby they also act at the apex point of weakly-linked long loops in foodwebs which are thought to be a stabilizing motif [9].

A less studied but perhaps equally important factor contributing to the role of top predators is the spatial nature of ecological systems. The effect of space has been studied in ecology for a long time in the context of classical metapopulation models [19, 20], island biogeography [21], and movement ecology [22]. However, spatially explicit studies of the dynamics of large foodwebs (meta-foodwebs) have only begun to appear recently [23, 24, 25, 26, 27, 28, 29, 30]. Based on earlier studies [31, 32, 33, 34] and the recent work of Pillai et al. [23, 35] we can expect that complex spatial environments aid the persistence of diverse communities by increasing the opportunities for generalists and enabling quick recovery from localized perturbations via the rescue effect.

In addition to beneficial consequences, spatial environments can also create intrinsic instabilities leading to complex and potentially dangerous spatiotemporal dynamics [36, 37, 30]. In [30] we introduced an efficient method for identifying such instabilities in large spatially explicit models of foodwebs dynamics, by combining a generalized modelling approach [38, 39] with a Master-Stability function method [40, 41].

A feature that is not captured in previous models is that top predators perceive the habitat differently than other species. Top predators tend to be large-bodied highly mobile species that need to roam large areas to cover their energy needs. While a landscape may appear as a complex network of many habitat patches from the perspective of a mouse it only consists of a single patch from the perspective of an eagle circling overhead.

The top-predator’s different perspective on the landscape is interesting because it means that the top predator cannot only link different biomass streams and close long loops in foodwebs by feeding on different species, it can also do so by simultaneously preying on populations in different habitat patches.

In the present paper we investigate the stability of meta-foodwebs that appear as a complex network of habitat patches to all but the top predator population, while the top predator regards the whole system as a single global patch. We refer to such ultra-mobile predators as global predators from now on. To study this system we use a generalized modelling approach and extend the master-stability function approach to meta-foodwebs [30] to accommodate a global top-predator.

We find that systems with a global top predator tend to be more stable than systems with the same number of species where all species are local. This effect is particularly pronounced when the global predator feeds on several species, and thus acts as a generalist. In this case the stability of the system with a global top predator exceeds the stability of the smaller system where the global predator is absent. Our results highlight the importance of global generalist top predators in safeguarding the systems stability against disturbances.

Model

We consider a system of $S$ local species on $N$ patches that disperse between these patches according to a network of dispersal routes, and $S^{\prime}$ global species that can feed on all patches simultaneously. Below we particularly consider $S^{\prime}=1$ in examples, while in the mathematical work we leave $S^{\prime}$ undefined for generality. We denote the biomass density of local species $i$ in patch $k$ by $X_{i}^{k}$ and the biomass density of the global species $i$ by $Y_{i}$. The dynamical equations of the model have the following form:

[TABLE]

where $G$ is the rate of growth by primary production, $M$ the rate of the loss by respiration and mortality, $F$ the rate of growth due to feeding on other species, $D$ the rate of biomass loss due to being eaten by other species, and $E_{i}^{kl}$ the rate of dispersal from patch $l$ to $k$. Furthermore, $\epsilon$ denotes the conversion efficiency of prey biomass into predator biomass. The subscripts indicate the affected populations. If two populations are affected, the order of the subscripts is such that biomass flows from the second to the first population. For instance, the function $D_{Y_{j}X_{i}^{k}}$ describes the biomass loss of the local population of species $i$ on patch $k$ due to predation by the global species $j$. For this generalized equation system, and without restricting the functions in the equation to a specific functional form, one can then derive a Jacobian matrix that describes the response of the system to perturbations. Because of the uncertainties that exist, this Jacobian matrix still contains a number of unknown (but directly interpretrable) parameters that describe the topology of the network, the scale of biomass turnover (scale parameters), and the strength of nonlinearities (exponent parameters).

We use the niche model [42] to generate foodwebs with realistic topologies. An example of a five-species niche web is shown in figure 8. We identify the niche value with the logarithm of the body mass, which allows us to implement allometric scaling of the biomass turnover rates. We then generate the scale and exponent parameters, taking the food web structure into account. The parameters are drawn at random from suitably defined intervals in accordance with ecological reasoning (see Methods, Tab. 2) and Fig. 8. The resulting sets of the generalized parameters closely follow previous publications, such as [38, 39, 43, 37].

We use random geometric graphs (RGG)[44] as a model for the spatial topology of the network of habitat patches (see Methods). Unless mentioned otherwise, we consider homogeneous systems, where all patches are characterized by the same parameters and have the same biomass densities. This assumption enables the construction of master stability functions (explained below) that allow us to efficiently study systems with arbitrary number of geographical habitat patches.

To construct the master stability function we first note that the Jacobian has the form

[TABLE]

where the indices $X$ and $Y$ refer again to local and global species. Furthermore, the block $\mathbf{J}_{X}$ has the internal structure

[TABLE]

familiar from systems without global predator [30], provided that the feeding rates of the global predator depend linearly on prey biomass. Here, $\otimes$ is a Kronecker product, $\mathbf{I}$ is an identity matrix of suitable size, $\mathbf{P}_{X}$ is the Jacobian of an isolated patch, $\mathbf{L}$ is the Laplacian matrix describing diffusion on the patch network, and $\mathbf{C}_{X}$ is the linearization of the dispersal terms [30]. More details, and the generalization to global predators with nonlinear feeding rates, are given in the Methods section. By exploiting the structure of the Jacobian its eigenvalues can be calculated by first finding the eigenvalues of the Laplacian matrix of the spatial network, and then computing the eigenvalues of a greatly reduced matrix in which the eigenvalues of the Laplacian appear as a parameter.

While we confine the details of the eigenvalue construction to the methods section, we remark that the procedure follows the same idea as the master stability function approach used in the study of syncronization [40, 41] but is applied to steady states as in [30] and is further modified in this paper to accomodate the global predators.

The master stability function approach provides several advantages. First and foremost, it separates the contributions of the spatial structure to the dynamics from those of the local food web dynamics. This allows clearer ecological insight into the source of observed instabilities.

Building on the separation of local and geographic structure one formulate precise conditions that the spatial structure must meet for a given food web to be stable. The function that represents these conditions is the actual master stability function after which this approach is named.

A byproduct that we use here, is a drastic reduction in the numerical effort for finding the eigenvalues of the Jacobian. In the present paper this reduction enables us to obtain solid statistics in the numerical investigations of large many-patch many-species meta-foodwebs.

To verify that the results described below also hold for systems with patches that violate the homogeneity assumption, we ran additional numerical analyses of systems with heterogeneous patches. In this case the spectrum of the Jacobian was determined by direct numerical diagonalisation of the full meta-foodweb Jacobian matrix. However, due to the large size of these matrices, we used smaller systems than for the homogeneous case.

Results

In the following we compare the mean stability of systems with a global top predator to that of systems with a local top predator or none at all. Then we evaluate how this stability depends on the sensitivity of the top predator to prey biomass, and on the number of prey populations.

Proportion of stable webs

We quantify stability as the proportion of stable webs in an ensemble of systems that correspond to the same model. A web is stable if all eigenvalues of the Jacobian have a negative real part. An ensemble is characterized by the number of local species $S$, of global species $S^{\prime}$ (which is 1 or 0), of spatial patches $N$, and a set of intervals for the generalized parameters, see Table 2. Each system of the ensemble is generated by randomly creating a spatial network (RGG), a niche foodweb with $S+S^{\prime}$ species, and a set of generalized parameters from the given intervals. For the foodwebs with a global predator the species with the highest niche value is selected as the global predator.

Figure 1 shows a comparison of the proportion of stable webs of systems with the top predator added as a global species vs as a local species. The proportion of stable webs is in all examined cases larger for the systems with the global top predator than in the systems with the local one, with the increase ranging between 20% and 200%. For the unstable systems, we evaluated whether the leading unstable mode is homogeneous or inhomogeneous. Homogeneous instabilities make a minor contribution to the number unstable webs. This is because there are many more inhomogeneous than homogeneous eigenmodes. Figure 2 shows that the proportion of stable webs remains almost unchanged when the global top predator is removed from the system.

Among the unstable systems, the proportion of homogeneous and inhomogeneous leading modes is also nearly the same. This is because the underlying eigenvalue problem for the stability of homogeneous systems can be separated into two types of solutions. One type represents the inhomogeneous eigenmodes where the global predator has the value 0 in the corresponding eigenvector. The other part represents the coupling of the global species to the spatial part of the system, which corresponds to the spatially homogeneous eigenvectors. When the global predator is removed, the spatially inhomogenoues part of the eigenvalue problem, which gives most of the eigenmodes, is essentially unchanged. Therefore the removal of the global predator has no strong impact on the proportion of stable webs. For more details on the eigenvalue problem see the Methods section.

To explore in more detail the influence of the global predator on stability we evaluate the change of stability with a key parameter related to predation by the top species, for systems with and without a global predator. This parameter, $\gamma$, describes the sensitivity of the top predator’s food intake to the total available prey biomass. It is known that it has a positive correlation with stability[38]. The reason is that the predator is close to saturation for small $\gamma$ and therefore can hardly control the population of its prey species.

At first, we use two types of three-species foodwebs, one with a generalist top species and one with a specialist top species (see figure 3 b and d). We compare the situation that the top species is global with the situations where it is local or completely absent. In the latter case, we have a two-species system. The proportion of stable webs in dependence of the exponent parameter $\gamma$ of the top species for the two types of three-species systems is shown in figure 3.

As before, the stability of the systems with the global top predator is higher than that with the local top predator. Furthermore, the stability of the systems with the generalist global top predator is also higher than that of the two-species system when $\gamma$ is larger than 0.25. With the specialist top predator stability is smaller than for the two-species system but approaches it for larger $\gamma$.

In the case of the three-species systems the generalist global species has an stabilizing effect, in contrast to the specialist one. However, from studying this small system it is not clear whether the stabilizing effect is due to the larger number of prey species or to the fact that the global species is an omnivore that feeds on different trophic levels. In order to investigate this issue, we also studied systems with more species, and we explored how stability depends on the number and trophic range of the top predator’s prey species in these larger foodwebs. Figure 4 a shows the proportion of stable webs for system with $S=10$ local species and a global predator in dependence of the number of prey species of the global predator. A larger number of prey species is correlated with a higher proportion of stable webs. When the global predator has four prey species the stability reaches that of the system without global predator and exceeds it for larger numbers of prey species.

In figure 4 b the proportion of stable webs is shown in dependence of the niche interval size and the number of prey species of the global predator. The niche interval size is defined as the difference between the highest and the smallest niche values of the global predator’s prey. Larger niche interval size implies that the global predator feeds on more trophic levels on average. The figure shows that the global predator has the strongest stabilizing effect when it feeds on a large number of species with relatively similar niche values. This means that feeding on different trophic levels is not a stabilizing factor in our system.

Heterogeneous systems

The results presented so far were obtained using a spatially homogeneous model. Making the systems heterogeneous makes them more realistic as natural systems always have heterogeneities. Since we cannot use the Master Stability Function approach any more, the numerical effort becomes much larger, as we have to perform a direct numerical evaluation of the eigenvalues. Heterogeneity is introduced by assigning the exponent parameters independently for every patch from the same uniform distributions as before (table 2). We use the same foodweb structure in every patch. Hence we get the same scale parameters in every patch. Note that our way of constructing heterogeneous systems ignores correlations between the exponent parameters in different patches and is therefore a stronger type of heterogeneity than we would expect from nature.

Figure 5 shows the proportion of stable webs. It is analogous to figure 1, but for smaller, heterogeneous five-patch systems. For comparison the figure also includes the data for homogeneous webs of the same size as the heterogeneous ones. One obvious difference between the heterogeneous and homogeneous systems is the overall lower proportion of stable webs, but this need not always be the case. Previous work has shown that heterogeneity can increase as well as decrease linear stability in generalized models [45].

Figure 6 shows the dependence of stability on the number of prey of the global predator. These results are analogous to those of figure 4, but for a heterogeneous system with $S=7$ local species. The plot shows the same trend as for the homogeneous system. A larger number of prey species of the global predator leads to a higher proportion of stable webs. Furthermore, the proportion of stable webs is always above that obtained without global predator. In our example system, the stabilizing effect of the global predator is thus stronger in the heterogeneous system than in the homogeneous system.

In general, our results for spatially heterogeneous systems agree qualitatively with those for homogeneous systems, suggesting that the insights gained from the study of homogeneous systems have a broad range of validity.

Discussion

Our study has shown that global top predators can have a considerable stabilizing effect on meta-foodwebs. This result was obtained for a broad range of models and parameter values, as we used a generalized modelling approach that does not depend on the specific functional form of the interaction terms. The average stability of meta-foodwebs with a global top predator was larger than that of systems where the top predator was made local. The difference is often of the order of a factor of two. This means that when a global predator is made local in a stable ecosystem, the chance is around 50 percent that the system becomes unstable. Furthermore, we found that the stability of meta-foodwebs increases when the global predator has more prey species, and in many cases exceeds that of a system that has no such top predator. This latter effect appears to be considerably stronger when the system is heterogeneous, i.e., when different habitats have different equilibrium biomass densities.

Ripple et al.[17] mention a variety of functions that large predators fulfill in ecosystems around the world. Top predators are important at controlling the abundance in the lower tropic levels, this includes mesopredators[46] as well as herbivores. A decline in big top predators often results directly in a release of mesopredators. Indirect effects can cascade through the whole foodweb. Due to the stronger mesopredator populations the abundance of prey species on lower tropic levels may decrease. Despite the release of intermediate predators a lower number of top predators can lead to a limited control of herbivores and therefore to a decrease in plant abundance. Large top predators are also thought to play an important role for the resilience of ecosystems against invading species. Therefore a loss of the big top predators can have severe consequences especially in the context of climate change.

The extent to which habitats are fragmented has important effects on the functioning and stability of ecosystems [47], and especially on top predators[46]. Far-ranging predators are negatively affected by landscape fragmentation because they are exposed most strongly to edge effects, for instance due to conflicts with humans [48]. Our results show that the stabilizing effect of large top predators depends on their capability to range freely. But this is hindered by landscape fragmentation. Even if a top predator was able to survive the fragmentation of its habitat it could lose its stabilizing capabilities. Previous results also show how habitat fragmentation can completely change the trophic role of top predators[49].

Taken together, the results highlight the importance of far-ranging top predators for ecosystem stability. The decline of those predators in natural systems is a serious concern. The extinction of predatory mammals like the Eurasian brown bear, lynx and wolf on the British isle are only one of many examples. While some of these extinctions are caused by climatic and environmental changes, anthropogenic impacts are the biggest factor. The ecological effects of such island extinctions cascade across ecological networks and take centuries to be fully realized [50]. The protection or reintroduction of large top predators is therefore of great importance. However, our study suggests that their large home ranges must also be restored if they shall have a considerable stabilizing effect.

Methods

Generalized Modeling

We start with the general model equations (1) and (2), assuming that they have a steady state. We normalize all biomass densities and functions to this steady state. For example, if the population density of a local species at the steady state is denoted as $X_{i}^{k*}$ then the normalized population density is $x_{i}^{k}=\frac{X_{i}^{k}}{X_{i}^{k*}}$, and the normalized growth function is $g(x_{i}^{k})=\frac{G(X_{i}^{k})}{G(X_{i}^{k*})}$. This process yields the normalized meta-foodweb model equations

[TABLE]

where $i=1,\ldots,S$ and $k=1,\ldots,N$, and

[TABLE]

with $i=1,\ldots,S^{\prime}$. The factors in front of the normalized functions are expressed in terms of scale parameters, the interpretation of which are listed in Table 1. Apart from the $\alpha$, which are the biomass flow rates, the other scale parameters are dimesionless and lie in the interval $[0,1]$. In our model, the scale parameters are related to the structure of the spatial topology and of the local foodweb, as detailed in the next section.

The first step in a linear stability analysis consists in calculating the Jacobian. In addition to the scale parameters, this Jacobian depends also on exponent parameters, which are the derivatives of the normalized functions with respect to the normalized biomass densities [38]. The exponent parameters (or elasticities) for the dynamics within a patch are

[TABLE]

The exponent parameters related to dispersal between patches are

[TABLE]

where $i\neq j$. The interpretation of the exponent parameters can be found as well in Table 1. Note that different steady states are characterized by different values of the scale and exponent parameters and thus can have different stability properties.

Although the functions and the steady states are unknown, the exponent parameters have an clear meaning in the context of the system, and their typical range of values can therefore be obtained from general considerations (see Gross et al.[38, 5]). This range of values is given in Table 2 for three types of models. The first model is shaped after systems with a Holling type 2 functional response and conventional, diffusive dispersal. The Standard model uses the largest meaningful range of exponent parameters, as specified in previous publications[5, 43, 37, 45]. The last model implements an adaptive type of dispersal, where the prey avoids its predators and the predators follow their prey.

Generation of meta-foodwebs

The spatial topologies used in this paper were random geometric graphs (RGGs)[44]. These graphs were generated by randomly placing $N$ patches on a unit square and connecting all pairs of patches that had a distance smaller than

[TABLE]

with the parameter $r^{\prime}$ radomly chosen between 1.4 and 2 for each RGG. This was done to cover an ensemble with a wider range of mean degrees. For the $N=10$ patch RGGs used in most of the simulations mean degrees range from approximately 1.8 to 9. The RGGs were generated with periodic boundary conditions to avoid edge effects. Only connected graphs were retained.

In most of the simulations we used spatial webs with $N=10$ patches. Test runs with larger numbers of patches showed no qualitative change in the results. A larger number of patches leads to a denser spectrum of the Laplacian eigenvalues, while the range of the eigenvalues is comparable.

If not specified otherwise we used the niche model [42] to create the feeding relations between the different species. This model assigns to each species a random niche value $n_{i}$ between 0 and 1 using the uniform distribution. We used a connectance of 0.25 and retained only connected niche webs. For systems with a global species, the species with the largest niche value was identified as the global predator. Niche values are correlated with bodymass, and we therefore implemented allometric scaling of the biomass turnover rate

[TABLE]

The structures of the spatial habitat network and the local foodweb imposes relations between several scale parameters: The total biomass turnover rate consists of the turnover due to dispersal between the spatial patches and the turnover due to local in patch dynamics including the interaction with the global predator. To obtain homogeneous steady states we assume linkwise dispersal. This means all the links between patches have the same biomass turnover rate. Therefore the biomass turnover due to dispersal in a patch is proportional to the number of links. As the number of incoming and outgoing links is identical for undirected spatial graphs, the biomass gain and loss by dispersal cancel each other out ($\nu=\rho$). We fulfil all these conditions by using for each species a random factor $a_{i}\in[0,1]$ to fix the local biomass turnover

[TABLE]

and the biomass turnover through one spatial link

[TABLE]

The remaining scale parameters are obtained via the following considerations: In our model, the basal species gain biomass exclusively due to primary production ($\delta=0$). The only biomass input of the top and intermediate species is provided by feeding on their respective prey species ($\delta=1$). In order to fix the values of the scale parameters $\beta$ and $\chi$ in a consistent manner, we assign first a random link strength $s_{ij}$ to each feeding link between predator $i$ and prey $j$ using the standard normal distribution. This also includes the feeding links of the global predator. We interpret this link strength as the biomass flowing through link $ij$ per unit time. The scale parameters $\beta$ and $\chi$ give the relative contribution of the different feeding links to the biomass loss resp. gain to a node of the foodweb, and they are therefore obtained from the link strengths $s_{ij}$ via the relations

[TABLE]

where for local species the sums are calculated over only the local species. A sample foodweb with 5 species is shown in figure 8. In order to distinguish between local and global losses the scale parameter $\xi$ is introduced. The fraction of loss by predation due to the global predator (index $g$) is calculated by

[TABLE]

with the sum over all species including the global predator. The fraction of biomass loss due to predation $\sigma$ is randomly chosen between 0 and 1 for all prey species separately, for top species $\sigma=0$.

Finally, the non-constant exponent parameters are drawn randomly from uniform distributions in the intervals given in table 2, using the mentioned three models.

Structure of the Jacobian

The Jacobian of the model represented by the equations (1) and (2) or, equivalently, by the equations (5) and (6), has the structure (3), which is illustrated in terms of smaller blocks of sizes $S\times S$, $S\times S^{\prime}$, and $S^{\prime}\times S^{\prime}$ in figure 7. The first block $\mathbf{J}_{X}$ contains the information about the dynamics of the $S$ local species and the matrix $\mathbf{J}_{Y}$ captures the dynamics of the $S^{\prime}$ global species, while the matrices $\mathbf{J}_{XY}$ and $\mathbf{J}_{YX}$ capture the interaction between the spatial web and the global species due to global predation.

The first block $\mathbf{J}_{X}$ itself can be constructed out of $N\times N$ block matrices of size $S\times S$. We denote their constituents as ${\mathbf{P}_{X}}^{k}$, which capture the local dynamics within patch $k$, ${\mathbf{C}_{X}}^{kl}$ and ${\mathbf{\hat{C}}_{X}}{}^{kl}$, which capture the contributions to the Jacobian due to dispersal from patch $l$ to $k$, and ${\mathbf{O}_{X}}^{kk}$ and ${\mathbf{O}_{X}}^{kl}$, which are nonzero when there is an indirect coupling of species within patch $k$ or between patches $k$ and $l$ due to nonlinear predation by global predators. With this notation, the first block of the Jacobian takes the form

[TABLE]

The structure of the spatial network web enters via the positions of the coupling matrices $\mathbf{C}_{X}^{kl}$ and $\mathbf{\hat{C}}_{X}^{kl}$. In the case of a spatially homogeneous steady state all generalized parameters are are equal between the different patches. Therefore the different matrices occuring in the Jacobian do in fact not depend on the patch indices, and the matrices $\mathbf{C}_{X}$ and $\mathbf{\hat{C}}_{X}$ become equal. Using the identity matrix $\mathbf{I}$, the Laplacian matrix of the spatial patch network $\mathbf{L}$, and the matrix $\mathbf{N}$, that is 0 on the diagonal and 1 everywhere else, the first block of the Jacobian can be written more compactly as

[TABLE]

where $\mathbf{P}^{\prime}_{X}=\mathbf{P}_{X}+{\mathbf{O}_{X}}^{kk}$. The remaining entries of the matrix in figure 7 are $\mathbf{J}_{Y}$, which captures the dynamics of the global species, and the matrices ${\mathbf{C}_{XY}}^{k}$ and ${\mathbf{C}_{YX}}^{k}$, which capture the coupling of the spatial patches to the global species due to predation of local species by global species.

Calculation of the eigenvalues

The eigenvalues are calculated in a similar way as in [30]. The system has two types of eigenvectors, the homogeneous and the inhomogeneous ones. The homogeneous eigenvectors have the same components for every spatial patch and a non-zero vector component corresponding to the global predator. The components of the inhomogeneous eigenvectors differ between the spatial patches while the vector component corresponding to the global predator vanishes. In the following, we calculate these two types of eigenvectors separately.

We start with the inhomogeneous eigenvectors by calculating the eigenvalues of the first block of the Jacobian $\mathbf{J}_{X}$. In the next step we show that certain eigenvalues of $\mathbf{J}_{X}$ are also eigenvalues of the full Jacobian $\mathbf{J}$.

The eigenvalues of $\mathbf{J}_{X}$ have to be calculated from equation (13). The main difference to the system without global predators is given by the coupling matrices $\mathbf{O}_{X}^{kl}$ and $\mathbf{O}_{X}^{kk}$ and is caused by the indirect coupling between local species due to global predation. The Laplacian $\mathbf{L}$ and the matrix $\mathbf{N}$ can be written as

[TABLE]

with $g_{i}$ being the degree of the vertex with index $i$, $A_{kj}$ being the adjacency matrix, and $\delta_{ij}$ being the Kronecker delta. Both matrices commute, so there is a common basis of eigenvectors,

[TABLE]

This allows us to replace $\mathbf{L}$ and $\mathbf{N}$ in (13) by their eigenvalues, giving the following reduced eigenvalue equation

[TABLE]

The full eigenvector can than be written as

[TABLE]

and the complete eigenvalue equation of the first block of the Jacobian is then

[TABLE]

Next we show that certain eigenvalues of $\mathbf{J}_{X}$ are also eigenvalues of $\mathbf{J}$. With a given eigenvector $\mathbf{v}$ corresponding to the eigenvalue $\lambda_{\kappa\mu}$ of $\mathbf{J}_{X}$ we find

[TABLE]

The constructed vector $(\mathbf{v},0)^{T}$ is an eigenvector if

[TABLE]

The Jacobians block in question $\mathbf{J}_{YX}$ has the structure

[TABLE]

Therefore it follows

[TABLE]

Accordingly we only need to show that

[TABLE]

We multiply the eigenvector equation from the left side with a row vector filled with 1,

[TABLE]

Because $(1,\dots,1)\cdot\mathbf{L}=0$, either $\kappa=0$ or $(1,\dots,1)\cdot\mathbf{p}=0$. Consequently the condition holds for the eigenvectors with $\kappa\neq 0$.

The multiplicity of the eigenvalue $\kappa=0$ of the Laplacian is the number of connected components of the considered spatial graph. In this paper we focus only on connected graphs, therefore the multiplicity of $\kappa=0$ is 1. The matrix $\mathbf{N}$ has the eigenvalue $\mu=N-1$ exactly once, and the only other eigenvalue is $\mu=-1$. With this eigenvalue, we have

[TABLE]

We multiply this equation again from the left side with a row vector filled with 1 and rearrange the equation,

[TABLE]

Therefore the condition given by equation (23) holds also for eigenvectors with $\mu=-1$.

The eigenvalues of $\mathbf{J}_{X}$ with $\kappa\neq 0$ and the eigenvalues with $\mu=-1$, which we have calculated using equation (16) are also eigenvalues of $\mathbf{J}$. These eigenvalues and eigenvectors correspond to the inhomogeneous eigenmodes of the spatial system.

Next we look at the spatially homogeneous eigenvectors. The Laplacian $\mathbf{L}$ and the matrix $\mathbf{N}$ share the common eigenvector $\mathbf{p}=(1,...,1)^{T}$ with the eigenvalue $\kappa=0$ and $\mu=N-1$ respectively. This eigenvector corresponds to a homogeneous behavior of the spatial system which is influenced by the global species. Therefore we are unable to calculate the eigenvalue of the Jacobian corresponding to the homogeneous eigenmode by taking only $\mathbf{J}_{X}$ into account. We calculate the corresponding eigenvalues by simplifying the eigenvalue equation of the Jacobian to the homogeneous case by using eigenvectors constructed as

[TABLE]

where

[TABLE]

and $\text{dim}(\mathbf{1})=N$. In order to show that $\mathbf{v}$ is indeed an eigenvector, we multiply it with the complete Jacobian:

[TABLE]

In summary, we reduced the eigenvector equation of the spatially homogeneous eigenmodes to

[TABLE]

The dimension of this eigenvalue equation does not depend on the size of the spatial web. The spatial web enters only via the number of patches $N$. In addition the equation does not depend on the spatial coupling $\mathbf{C}_{X}$ between local species due to dispersal.

This eigenvalue equation gives us the spatially homogeneous eigenmodes, where all the patches behave the same way. The local populations do not exchange biomass between patches. The global predators only react to the total available biomass in the complete system and are therefore insensitive to the topology of the spatial web. Hence the coupling between global and local species is described by the local eigenmode of the system.

Ensemble sizes

For the systems in figures 1, 2 and 3 we used ensembles of $10^{6}$ webs per class.

In figure 4 we used an ensemble of more than $10^{7}$ webs. The stability baseline of the system without the global predator in figure 4 a was calculated with an ensemble of $10^{6}$ webs. In figure 4 b only bins with sufficient statistics are shown. Statistics gets worse for more extreme combinations of prey numbers and niche interval sizes.

For the systems in figure 5 we used ensembles of various sizes ranging from $3.1\cdot 10^{6}$ for the more stable systems to over $1.3\cdot 10^{7}$ for the least stable systems. Note that for the three least stable classes not even one single stable web was found.

In figure 6 an ensemble of more than $3.5\cdot 10^{7}$ webs was used. The stability baseline of the system without the global predator was calculated with an ensemble of $5\cdot 10^{6}$ webs.

The proportion of stable webs (PSW) is calculated by

[TABLE]

where $N$ is the total number of tested webs and $N_{\text{stable}}$ is the number of stable webs found in the given sample set. The standard deviation of the PSW is estimated by

[TABLE]

Error bars are only shown when the standard derivation is visible on the scale of the plot. In figure 3 the error bars would be smaller than the line width used in the plot.

Data Availability

This study did not involve any underlying data. All the information needed to generate the data analyzed during this study is included in this published article.

Acknowledgements

We acknowledge support by the German Research Foundation (DFG research unit 1748, contract number Dr300/13-2) and by the EPSRC (EP/N034384/1). For coverage of the publication costs, we acknowledge support by the German Research Foundation and the Open Access Publishing Fund of Technische Universität Darmstadt.

Author information

Affiliations

Institute of Condensed Matter Physics, Technische Universität Darmstadt,

Darmstadt, Germany

Andreas Brechtel & Barbara Drossel

UC Davis, Department of Computer Science, 1 Shields Av

Davis, Ca 95616, USA

Thilo Gross

Contributions

All authors contributed to and reviewed the manuscript. A.B. performed the simulations and prepared the figures.

Competing Interests

The authors declare that they have no competing interests.