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

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 generalized 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. These results highlight the importance of top predator conservation.

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[5,9-11]. The topic of foodweb stability has gained renewed urgency by recent papers raising concerns that past perturbations have not only led to the direct loss of species but have also left the respective systems more fragile. This increases the impact of future perturbations, leading to further and possibly bigger losses[12,13]. Habitat loss, overexploitation, and numerous other stress factors have caused global n class="Disease">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. Although they constitute typically only a tiny fraction of the biomass[15], top predators have been found to have a strong impact on stability[5,16,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 importance of spatial structure for ecosystem dynamics is generally acknowledged[19]. Mathematical modelling has been focussing for many years on classical metapopulation models[20,21], island biogeography[22], and movement ecology[23]. However, spatially explicit studies of the dynamics of large foodwebs (meta-foodwebs) have only begun to appear recently[24-32]. Based on earlier studies[33-36] and the recent work of Pillai et al.[25,37] 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[32,38,39]. In a recent publication[32] we introduced an efficient method for identifying such instabilities in large spatially explicit models of foodwebs dynamics, by combining a generalized modelling approach[40,41] with a Master-Stability function method[42,43].

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, which is essentially a linear stability analysis that is based on generic expressions for the biomass turnover rates and elasticities at steady states[40]. In order to accommodate global top predators, we extend the master-stability function approach to meta-foodwebs[32]. This approach allows us to disentangle the spatial habitat network from the local foodweb when calculating the eigenvalues of the Jacobian[42].

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′ global species that can feed on all patches simultaneously. Below we particularly consider in examples, while in the mathematical work we leave S′ undefined for generality. We denote the biomass density of local species i in patch k by and the biomass density of the global species i by Y. The dynamical equations of the model have the following form: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 the rate of dispersal from patch l to k. Furthermore, 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 describes the n class="Disease">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[44] to generate foodwebs with realistic topologies. In this model each species is assigned a random niche value and a feeding interval the center of which lies below the niche value. A species feeds on all species with niche values inside its feeding interval. An example of a five-species niche web is shown in Fig. 1. 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, Table 1 and Fig. 1). The chosen intervals for the generalized parameters follow closely previous publications[39-41,45].

Figure 1

Illustration of the assignment of the scale parameters for a local foodweb with 5 species. In the circles, the species index and the niche value are given. On the links, the relative link strength s is given as well as the scale parameters and γ. According to the relations (15) the parameters β add up to 1 for a given prey (e.g., species number 2 being eaten by 3 and 4), and the parameters add up 1 for a given predator (e.g., species 4 feeding on 1 and 2).

Table 1

Constant values and distribution intervals of the exponent parameters.

	ϕX/ϕY	ψX/ψY	μX/μY	γX/γY	λX/λYX	ω X	ωˆ_X	κ X	κˆ_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]

In the case of the ‘Adaptive’ model the parameters and depend on the feeding links of the local foodweb. If species i is the prey of species j then is randomly assigned a value inside the interval and is otherwise set to 0. This means the dispersal of prey is stimulated if predators are present. Additionally the parameter is randomly assigned a value inside the interval [0, 1] if species j does feed on species i and is set to 0 otherwise. This corresponds to predators that prefer to disperse into patches with higher prey abundance.

We use random geometric graphs (RGG)[46] 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 formwhere the indices X and Y refer again to local and global species. Furthermore, the block J has the internal structurefamiliar from systems without global predator[32], provided that the feeding rates of the global predator depend linearly on prey biomass. Here, is a Kronecker product, I is an identity matrix of suitable size, P is the Jacobian of an isolated patch, L is the Laplacian matrix describing diffusion on the patch network, and C is the linearization of the dispersal terms[32]. 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 synchronization[42,43] but is applied to steady states as in one of our previous publications[32] 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′ (which is 1 or 0), of spatial patches N, and a set of intervals for the generalized parameters, see Table 1. Each system of the ensemble is generated by randomly creating a spatial network (RGG), a niche foodweb with S + S′ 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 2 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 3 shows that the proportion of stable webs remains almost unchanged when the global top predator is removed from the system.

Figure 2

Proportion of stable webs for 10-patch systems with a global top predator (red) and systems with a local top predator (blue), for the three models specified by the parameter intervals listed in Table 1, and for different species numbers . The horizontal line separates the small proportion of unstable webs with the leading eigenvalue corresponding to a homogeneous mode from the larger proportion corresponding to an inhomogeneous (pattern-forming) mode. The figure shows that making a global predator local decreases stability and that the majority of instabilities is due to pattern-forming instabilities.

Figure 3

Same as Fig. 2, but in the system without the global predator no local predator is added. The proportions of stable webs and both types of instabilities are very similar in the systems with and without the global predator. Comparison with the previous figure shows that adding a global predator does not decrease stability while adding a local predator does decrease stability.

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 n class="Chemical">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, γ, 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[40]. The reason is that the predator is close to saturation for small γ 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 Fig. 4b,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 γ of the top species for the two types of three-species systems is shown in Fig. 4.

Figure 4

Proportion of stable webs for 10-patch Holling type 2 systems as function of the sensitivity to the total amount of prey γ available to the top species (which is species number 3) of a 3-species omnivore module (top) and a 3-species food chain (bottom), with the top species being either global (red) or local (blue). Additionally, the mean proportion of stable webs of the 2 species system without the species 3 is shown as the horizontal dashed line. The data show that systems with global top species are more stable than those with local top species, and that the generalist module is more stable than the three-species food chain. With a generalist top species it is even more stable than the two-species system.

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 γ is larger than 0.25. With the specialist top predator stability is smaller than for the two-species system but approaches it for larger γ.

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 5a shows the proportion of stable webs for system with 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.

Figure 5

(a) Proportion of stable webs for a 10-patch Holling type 2 system with local species and a global predator for the different numbers of prey species of the global predator. The dashed line represents the proportion of stable webs of the system when the global species is removed. The system’s stability increases when the global species has more prey species. For more than 4 prey species the system with the global predator exceeds the stability of the system without it. (b) Proportion of stable webs for the same system as in (a), showing additionally he dependence on the niche interval size. The niche interval size is the difference between the largest and the smallest niche values of global predator’s prey species. White color corresponds to the proportion of stable webs of local system when the global species is removed, red color corresponds to higher and blue color to lower stability.

In Fig. 5b 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 1). 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 6 shows the proportion of stable webs. It is analogous to Fig. 2, 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[47].

Figure 6

Proportion of stable webs for heterogeneous five-patch systems with a global top predator (red) and systems with a local top predator (blue), for the three models specified by the parameter intervals listed in Table 1, and for different species numbers . The figure shows that making a global predator local decreases stability also for the more general, inhomogeneous meta-foodwebs. The lighter shaded bars indicate the proportion of stable webs of the corresponding homogeneous systems. In the model shown here, the heterogenoeus meta-foodwebs have a lower average stability than their homogeneous counterparts.

Figure 7 shows the dependence of stability on the number of prey of the global predator. These results are analogous to those of Fig. 5, but for a heterogeneous system with 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.

Figure 7

Proportion of stable webs for heterogeneous five-patch systems with a global top predator and local species for the Holling 2 specified by the parameter intervals listed in Table 1. The figure shows that a higher number of prey species leads to more stable systems. This is in agreement with the homogeneous case (see Fig. 5a). The stability of the system without the global top predator is indicated by the dashed line, showing that it is less stable than the system with global top predator.

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 three different classes of models and a broad range of parameter values. This broad scope of the investigation was made possible by using 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, as shown in Figs 2 and 6. 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. From Fig. 2, we concluded furthermore that the stabilizing effect consists to a large part in the suppression of pattern-forming instabilities.

In a more detailed investigation, 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, see Figs 4 and 7. This latter effect appears to be considerably stronger when the system is heterogeneous, i.e., when different habitats have different equilibrium biomass densities, as can be seen from a comparison of Fig. 7 with Fig. 4.

Ripple et al.[17] mention a variety of functions that large predators fulfill in ecosystems around the world. According to these authors, top predators are important at controlling the abundance in the lower trophic levels, this includes mesopredators[48] 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 trophic 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.

The extent to which habitats are fragmented has important effects on the functioning and stability of ecosystems[49], and especially on top predators[48]. 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[50]. Our results show that the stabilizing effect of large top predators depends on their capability to range freely, i.e., on their being global predators and not local predators. But landscape fragmentation makes far-ranging predators more local. In the light of our results this means that 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[51].

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[52]. 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 then the normalized population density is , and the normalized growth function is . This process yields the normalized meta-foodweb model equationswhere and , andwith . The total amount of biomass that is available for predation by a given predator is defined as:where the functions r are the normalized response of a predator to the abundance of a certain prey population. The factors in front of the normalized functions are expressed in terms of scale parameters, the interpretation of which are listed in Table 2. Apart from the α, which are the biomass flow rates, the other scale parameters are dimesionless and lie in the interval . 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.

Table 2

Generalized parameters used to describe the meta-foodweb.

Parameter	Interpretation
Exponent
ϕ	Sensitivity of primary production to own population density
γ	Sensitivity of predation to total available prey density
λ	Exponent of prey switching
ψ	Sensitivity of predation to predator density
μ	Exponent of closure
ω	Sensitivity of dispersal to own the species population density in the output patch
ω˜	Sensitivity of dispersal to own species population density in the target patch
κ	Sensitivity of dispersal to population density of some other species in the output patch
κˆ	Sensitivity of dispersal the population density of some other species in the target patch
Scale
α	Biomass flow
σ	Fraction of biomass loss due to predation
σ˜	Fraction of biomass loss due to respiration
ξ	Fraction of loss by predation due to global species
ξ˜	Fraction of loss by predation due to local species
β	Relative contribution to biomass loss due to predation by a certain predator
δ	Fraction of local growth by predation
δ˜	Fraction of local growth by primary production
χ	Relative contribution of population as prey to a certain predator
ν	Fraction of total biomass gain due to dispersal
ν˜	Fraction of total biomass gain due to predation and primary production
ρ	Fraction of total biomass loss due to dispersal
ρ˜	Fraction of total biomass loss due to predation and respiration

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[40]. The exponent parameters (or elasticities) for the dynamics within a patch are

The exponent parameters related to dispersal between patches arewhere . The interpretation of the exponent parameters can be found as well in Table 2. 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.[5,40]). This range of values is given in Table 1 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,39,45,47]. 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)[46]. These graphs were generated by randomly placing N patches on a unit square and connecting all pairs of patches that had a distance smaller thanwith the parameter r′ 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 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 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[44] to create the feeding relations between the different species. This model assigns to each species a random niche value n 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

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 (). We fulfil all these conditions by using for each species a random factor to fix the local biomass turnoverand the biomass turnover through one spatial link

The remaining scale parameters are obtained via the following considerations: In our model, the basal species gain biomass exclusively due to primary production (). The only biomass input of the top and intermediate species is provided by feeding on their respective prey species (). In order to fix the values of the scale parameters β and in a consistent manner, we assign first a random link strength s 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 β and 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 via the relationswhere for local species the sums are calculated over only the local species. A sample foodweb with 5 species is shown in Fig. 1. In order to distinguish between local and global losses the scale parameter is introduced. The fraction of loss by predation due to the global predator (index g) is calculated bywith the sum over all species including the global predator. The fraction of n class="Disease">biomass loss due to predation σ is randomly chosen between 0 and 1 for all prey species separately, for top species .

Finally, the non-constant exponent parameters are drawn randomly from uniform distributions in the intervals given in Table 1, 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 , , and in Fig. 8. The first block J contains the information about the dynamics of the S local species and the matrix J captures the dynamics of the S′ global species, while the matrices J and J capture the interaction between the spatial web and the global species due to global predation.

Figure 8

Structure of the Jacobian given in equation (3). The block J consists of the blue and green parts, the block J consists of the yellow parts, the block J consists of the orange blocks and the block J is shown in red.

The first block J itself can be constructed out of block matrices of size . We denote their constituents as , which capture the local dynamics within patch k, and , which capture the contributions to the Jacobian due to dispersal from patch l to k, and and , 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

The structure of the spatial network web enters via the positions of the coupling matrices and . 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 C and become equal. Using the identity matrix I, the Laplacian matrix of the spatial patch network L, and the matrix N, that is 0 on the diagonal and 1 everywhere else, the first block of the Jacobian can be written more compactly aswhere . The remaining entries of the matrix in Fig. 8 are J, which captures the dynamics of the global species, and the matrices and , 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 Brechtel et al.[32]. 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 J. In the next step we show that certain eigenvalues of J are also eigenvalues of the full Jacobian J.

The eigenvalues of J have to be calculated from equation (18). The main difference to the system without global predators is given by the coupling matrices and and is caused by the indirect coupling between local species due to global predation. The Laplacian L and the matrix N can be written aswith g being the degree of the vertex with index i, A being the adjacency matrix, and δ being the Kronecker delta. Both matrices commute, so there is a common basis of eigenvectors,

This allows us to replace L and N in (18) by their eigenvalues, giving the following reduced eigenvalue equation

The full eigenvector can than be written asand the complete eigenvalue equation of the first block of the Jacobian is then

Next we show that certain eigenvalues of J are also eigenvalues of J. With a given eigenvector v corresponding to the eigenvalue of J we find

The constructed vector is an eigenvector if

The Jacobians block in question J has the structure

Therefore it follows

Accordingly we only need to show that

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

Because , either or . Consequently the condition holds for the eigenvectors with .

The multiplicity of the eigenvalue 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 is 1. The matrix N has the eigenvalue exactly once, and the only other eigenvalue is . With this eigenvalue, we have

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

Therefore the condition given by equation (28) holds also for eigenvectors with .

The eigenvalues of J with and the eigenvalues with , which we have calculated using equation (21) are also eigenvalues of J. These eigenvalues and eigenvectors correspond to the inhomogeneous eigenmodes of the spatial system.

Next we look at the spatially homogeneous eigenvectors. The Laplacian L and the matrix N share the common eigenvector with the eigenvalue and 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 J into account. We calculate the corresponding eigenvalues by simplifying the eigenvalue equation of the Jacobian to the homogeneous case by using eigenvectors constructed aswhereand . In order to show that v is indeed an eigenvector, we multiply it with the complete Jacobian:

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

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 C 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 Figs 2, 3 and 4 we used ensembles of 106 webs per class.

In Fig. 5 we used an ensemble of more than 107 webs. The stability baseline of the system without the global predator in Fig. 5a was calculated with an ensemble of 106 webs. In Fig. 5b 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 Fig. 6 we used ensembles of various sizes ranging from for the more stable systems to over for the least stable systems. Note that for the three least stable classes not even one single stable web was found.

In Fig. 7 an ensemble of more than webs was used. The stability baseline of the system without the global predator was calculated with an ensemble of webs.

The proportion of stable webs (PSW) is calculated bywhere N is the total number of tested webs and Nstable is the number of stable webs found in the given sample set. The standard deviation of the PSW is estimated by

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