Study systems

Four natural food webs from different habitats are included in the analysis; Broadstone Stream56 (running water), Kelp bed57 (marine), Montane Forest57 (terrestrial) and Phytotelmata58 (freshwater) (Fig. 3, Supplementary Table 1). Nodes in the natural networks represent true species or trophic species, that is, functional or taxonomical groups that share resource and consumer nodes, and are here simply referred to as ‘species’. Species richness and connectance (proportion of possible trophic links that is realised—excluding any cannibalistic links) range from 17 to 28 and from 0.06 to 0.17, respectively (Supplementary Table 1). For the food webs Broadstone Stream and Montane Forest the lowest trophic level are the primary consumers (herbivores) due to low resolution at the primary producer level. Each species on the lowest trophic level competes directly with all other species on that level and intraspecific competition is present in all species.

Furthermore, we generate pyramidal model food webs with 12 and 24 species distributed over three trophic levels; primary producers, primary consumers and secondary consumers. The ratio between the numbers of species at the different trophic levels is 3:2:1. Trophic links are distributed randomly, but with some constraints: consumers must feed on at least one resource from the trophic level directly below themselves and cannibalism is not allowed. Secondary consumers can be omnivorous, that is, feeding not only on the primary consumers but also on the primary producers. Connectance is randomly drawn from the uniform interval [0.05 0.30]. Each primary producer species competes directly with all other primary producer species and intraspecific competition is present in all species.

An additional set of model food webs are generated using the PNM59. The PNM uses a Gaussian formulation for the probability that species i eats species j. That is, there is always a non-zero probability that a species feeds on another species, and the closer the niche position of the resource is to the optimal diet position of the consumer, the higher the probability. This allows for gaps in the diet width of consumers. There are no constraints on trophic links except that cannibalism is not allowed. The size of the PNM food webs ranges between 12 and 21 species and connectance ranges between 0.02 and 0.35. Intra- and inter-specific competition is modelled as in pyramidal food webs.

Finally, we generate a set of antagonistic bipartite model networks containing 12 basal and 12 herbivore species. We vary the degree of nestedness in these networks using the methods developed by Thébault and Fontaine49. Connectance, intra- and inter-specific competition is modelled as in pyramidal and PNM food webs.

The analysed systems represent a broad spectrum of antagonistic network types and topologies; bipartite and pyramidal networks, networks generated by the probabilistic niche model algorithm and natural food webs from different habitats. Niche models like the PNM have been widely used to generate antagonistic networks. We also include pyramidal networks in our analysis since these networks have relatively discrete trophic levels, in agreement with the pattern found in real food webs60.

Network dynamics

Network dynamics are described by a generalized Rosenzweig–MacArthur model22,61:

Here dN i /dt is the rate of change of the density of species i with respect to time in a network with S species, b i (t) is the per capita growth (mortality) rate for primary producers (consumers) of species i at time t and is the per capita effect of species j on the per capita growth (mortality) rate of species i. The interaction coefficients, , describes competitive and predator–prey (trophic) interactions. The functional response of predators to the densities of their prey is non-linear of type II (hyperbolic) or III (sigmoidal). Equation (1) can then be rewritten in a more detailed form as:

Here competitive interactions are described by the first sum; a ij and a ii (j=i) is the per capita strength of inter- and intra-specific competition, respectively. The second and third sum describe the loss and gain of species i’s density due to consumption. L(i) is the set of species being competitors to species i, C(i) is the set of species being consumers of resource species i, R(j) is the set of species being resources of consumer species j and R(i) is the set of species being resources of consumer species i. Further, h ij is the ‘preference’ (feeding effort) and α ij is the intrinsic attack rate of consumer species j on resource species i. T ij is the time needed for consumer i to catch and consume resource j (handling time), e ij is the rate at which resource j is converted into new consumer i (conversion efficiency) and q ij (capture exponent) is an interaction specific positive constant. The capture exponent (q) and handling time (T) determine the type of functional response: q=1 and T=0 gives a type I (linear) functional response; q=1 and T>0 gives a hyperbolic type II functional response and increasing q gradually from 1 to 2 gives an increasingly sigmoidal type III functional response.

Environmental variability is introduced in species’ intrinsic growth or mortality rates as

where is the mean per capita growth or mortality rate of species i and ɛ i (t) is a stochastic process defined by

Here ɛ ik are independent and uniformly distributed stochastic values on the interval [−1 1], yielding a variance in ɛ of 0.333. Extreme values are as likely as the mean value of 0. The correlation among species in their response to environmental fluctuations (ρ), that is, the correlation between time series of ɛ ik is generated as in ref. 22. We explore two scenarios: high (ρ=0.75) and low (ρ=0.25) correlation among species in their response to environmental variation. To get continuous time series of stochastic values, ɛ ik are then interpolated using a cubic spline interpolation method implemented in the Matlab function interp1. The resulting set of coupled ordinary differential equations are integrated over time using a variable order solver based on numerical differentiation formulas implemented in Matlab’s ode15s (model computer code available on request). We also indirectly account for demographic stochastic variation and other processes that affects small populations disproportionately by defining quasi-extinction thresholds (see Simulation).

Adaptive rewiring and adjustment of feeding efforts

It is plausible that extinction thresholds, set by processes like demographic stochasticity and Allee effects, could be higher than the density thresholds causing consumers to rewire. We therefore assume that rewiring and compensatory change in feeding efforts takes place in connection with the extinction of preferred prey species, and not before. Hence, an adaptive consumer adjusts its resource preferences (feeding efforts), h ij , when a resource in its original set of resources goes extinct. If the extinct resource was the preferred resource, another of its original resources is randomly chosen as preferred and h ij is updated for all its remaining resources so that h ij still sums to 1. If a non-preferred resource goes extinct, the remaining non-preferred resources will have their preferences increased while the preferred resource is unaffected. Finally, each time an adaptive consumer loses its last resource it rewires to a completely new resource from a trophic level below itself, thus forming a novel trophic link.

Consumers with similar behavioural, physiological and morphological traits are likely to use similar types of resources. Thus, the new resource is randomly chosen from the set of resources belonging to the consumer which has the most similar set of prey as that of the rewiring consumer itself. Similarity of consumers with respect to their sets of prey species is quantified using the Jaccard index of similarity (Fig. 1). We assume that rewiring entails a cost in the sense that the consumer is less effective in exploiting its new resource. Two cost scenarios are explored: a 75% (high cost) or a 25% (low cost) reduction of a consumer’s intrinsic attack rate, α ij , on the new resource species.

Parameterization

Mean intrinsic growth rates of primary producers, , are set to 1. Mean intrinsic mortality rates of consumers, , are randomly drawn from the uniform distribution [−0.01 −0.001] and ordered such that secondary consumers have lower intrinsic mortality rates than primary consumers (herbivores), that is, |b primary consumers |>|b secondary consumers |. This is because body mass of species often increases with trophic level62,63 and larger body mass often confers lower mortality rates64. Direct interspecific competition among primary producers, a ij , is set to 0.5. There is no direct competition among consumer species—they only compete indirectly through shared resources. Intraspecific competition, a ii , is set to 1 for primary producers and to 0.01 for consumers.

The intrinsic attack rate, α, and handling time, T, are assumed to be equal for all consumers and set to 1.5 and 0.2, respectively. The type of functional response is controlled by the capture coefficient q, which is set to 1 for a functional response of type II and 1.25 for type III. We vary the proportions of consumers with hyperbolic (type II) and sigmoidal (type III) functional responses as both types are likely to be present in a natural food web65 and because the type of functional response of predators can have important consequences for the risk of extinction of their prey66. Conversion efficiency, e, is set to 0.2 for feeding interactions between adjacent trophic levels and to 0.02 for omnivorous links which are considered less efficient.

Resource preferences (feeding efforts) of consumers, h ij , sum up to 1. Two scenarios with respect to prey preferences of consumers are analysed: (1) consumer species either show strong preference for one of their prey species (specialists) or (2) equal preference for each of their prey species (generalists). Specifically, in the specialist scenario consumers with more than one resource species show high preference (h ij,high =0.9) for one of them (randomly assigned) while the preferences for the other resources are low and equal; h ij,low =(1−h ij,high )/(no. of resources−1). In the generalist scenario consumers with more than one resource species show equal preference for each of them; h ij =1/(no. of resources).

To verify the generality of our results we perform an expanded analysis covering a larger parameter space using LHS design67 with randomised model parameters. This is done for a subset of the webs: 12-species pyramidal webs, PNM webs with 10 consumer species and the Phytotelmata food web. The selection of networks and scenarios used for this expanded analysis is based on representativeness and computational feasibility. Using LHS we first increase the parameter range analysed by setting the lower and upper interval limits of each of the 10 parameters to 0.5 and 2 times the constant parameter value used in the main analysis, respectively. Next, these larger intervals are divided into 10 equally sized sections (intervals). From each section a Latin Hypercube sample is randomly drawn to create a 10-dimensional Latin Hypercube that will hold 10 sets of parameter mean values, representing the larger parameter space (Supplementary Table 3). Model parameter values are then randomly drawn from uniform intervals with the upper and lower interval limits set to the LHS parameter mean ±10%. The parameter sets that do not result in persistent networks (having species densities above predefined quasi-extinction thresholds) for at least 20,000 time units during a deterministic simulation, within 10,000 trials, are discarded. The remaining parameter sets are used in the analysis. We use the lhs-package (version 0.1) in R68 to generate the Latin Hypercube.

Simulation

For the pyramidal model food webs, 12 and 24 species networks are simulated. Here, the number of consumer species is 6 and 12, respectively, (consumers constitutes 50% of the species in the webs). For the PNM food webs the number of consumers may not be equal in webs of the same size (total number of species). Here webs with 6 and 10 consumer species are simulated. The median number of species in these webs is 13 and 18 species, respectively. As in the pyramidal food webs, the consumers constitute 50% of all species in the bipartite model networks. Here networks with 24 species (that is, 12 consumer species) are simulated. The number of species in the 4 natural food webs are 17, 21, 27 and 28, respectively (Supplementary Table 1). It was only possible (computationally feasible) to generate persistent (see below) models of the natural food webs for the strong prey preference (specialist) scenario. For all simulated network types and scenarios, the fraction of adaptive (F Adap ) and non-adaptive (1-F Adap ) consumers as well as the fraction of consumer species having a functional response of type II (F FR2 ) or type III (1-F FR2 ) is varied between 0 and 1 (in discrete steps corresponding to one consumer species).

For each consumer trait combination within each scenario, 50 persistent replicates are generated. Networks are considered persistent if they remain feasible (having species densities above predefined quasi-extinction thresholds) for at least 20,000 time units during a deterministic simulation. As initial densities for this simulation we use the equilibrium densities from the corresponding system with a linear functional response (type I). These persistent networks are then integrated for another 10,000 time units (corresponding to ∼8000 generations for the primary producers and ∼30 generations for the longest lived consumer species in the theoretical webs), now with environmental stochasticity affecting intrinsic growth and mortality rates of species. As initial densities for these simulations we use the final densities from the deterministic simulation. Demographic stochasticity is implicitly modelled by defining quasi-extinction thresholds. The quasi-extinction thresholds are set to 0.002, 0.0002 and 0.00002 for primary producers, primary consumers and secondary consumers, respectively. However, the absolute levels are not of primary importance here since our aim is to compare relative extinction risks among the different scenarios explored, not to predict absolute extinction risks. If the density of a species falls below the quasi-extinction threshold it is considered extinct and is deleted from the system. A consumer species that has lost all prey species (and cannot rewire) will die from starvation. We assume that this starvation is not an instantaneous event, but a process that takes some time. During this time intermediate consumer can still interact with other species (that is, their predators), while top consumers cannot (at this point they will be isolated from all other species). Therefore intermediate consumers are kept in the system for 50 time steps following final prey loss, while top consumers are removed immediately.

As a complement to the simulations we also study a simple 4-species system using analytical methods to get further insight of possible mechanisms involved. In this analysis, we also explore how the strength of interspecific competition among primary producers and the strength of intraspecific competition in consumers modulate the effects of rewiring.

Statistical analysis

To investigate potential influences of food web topology on the effect of rewiring on food web persistence a regression tree analysis69 is performed. To this end, a set of network properties is calculated for each food web (Supplementary Table 4). For each food web type, only properties that show variation among replicates are included in the analysis. Depending on food web type some of these properties show no variation among replicate webs, for example, the proportion of basal species is constant for all pyramidal food webs. In addition, nestedness is only included in the analysis of the bipartite webs, as calculating nestedness for networks with more than two trophic levels is far from trivial (but see ref. 70). The fraction of adaptive consumers (F Adap ) and the fraction of consumers with a functional response of type II (F FR2 ) are included in the analysis to better understand the relative importance of these predictors versus the food web topology. Consequently, the number of potential explanatory variables varies between the food web types (22 for PNM, 17 for pyramidal and 12 for bipartite food webs, Supplementary Table 4). All analyses described below are performed for each food web type independently.

Inevitably, many network properties are related. To avoid an overly redundant set of explanatory variables we perform a correlation analysis. Network properties are removed so that no absolute value of pairwise correlation |ρ|>0.65 remains in the set. From clusters of network properties with high correlations between them, the most commonly used variable is used, for example, connectance among a larger group of link related properties. A screening analysis shows that our selection criteria of ‘convention’ did not result in the exclusion of any exceptionally powerful predictor (exceptional as compared with the power of the other properties within a certain correlation cluster). For the screening analysis we use generalized additive models (GAM) for estimating the explanatory power of each network property independently. A few exceptions to the cutoff rule of 0.65 are allowed in cases when two highly correlated properties are both deemed of great interest: for the bipartite food webs both nestedness (significant) and connectance (ρ=0.88) are included and for the PNM webs both average trophic height and connectance (ρ=0.75) are included in the regression analysis. However, by dropping either one of the properties from the set of explanatory variables and re-running the regression tree analysis with the diminished data set, we make certain that these exceptions do not negatively affect the performance of the analysis.

After the selection process, the set of explanatory variables involved in the regression tree analysis consist of eight network properties for the PNM webs, seven properties for the pyramidal webs and five properties for the bipartite webs (Supplementary Table 5). For each food web type four scenarios are analysed: specialist and generalist consumers; low and high rewiring cost (Supplementary Figs 6–8). All analyses are performed in R (ref. 68; version 3.1.1 for correlation and screening analysis and 3.0.1 for the regression tree analysis). Package mgcv (version 1.7-23) is used for the GAMs of the screening analysis, using function gam with basis cr and basis dimension, k, set to 4 (see documentation for mgcv). Package mvpart (version 1.6-1) is used for the regression tree analysis, using function mvpart applying the 1SE rule, allowing a minimum of 10 data points in a leaf and running 1000 10-fold cross validations (see ref. 69 and documentation for mvpart).