Introduction

Invasive species are one of the leading threats to biodiversity in the world1 and are known to impact species functioning across many trophic levels and ecological guilds2. Accurately assessing the impact of invasive species, however, is a complicated and resource-intensive process3. By their nature, ecological communities are complex and vary greatly in time. Even a herbivorous insect can cause indirect effects on other non-plant native species: for example, a seed feeding fly, released to control a pasture weed in the USA, caused a large increase in the abundance of predatory mice4. Species which feed on or compete with a wide range of species and taxa, therefore, may have larger and more widespread effects on large parts of the food web. Ecological networks which characterize and quantify the species (network nodes) and their interactions (links) can provide a powerful tool for studying the impacts of introduced species. Several authors have studied how the impacts of invasive species on ecological networks can be modelled and propagated across trophic levels5,6,7. However, the ultimate test of these methods is whether they can produce results reliable enough to inform conservation practice and policy.

Following7, we use a generalized model8, where the biomass of the species of the network change in time according to general functions, to explore the effect of alien species on an island food web. The generalized model allows analyses of ecologically realistic network interactions and measure the addition of virtual invasive species. The model finely resolves the species, which leads to a system that is an order of magnitude larger than previous systems studied with this approach, allowing for a detailed analysis6,7.

Our objective is to determine how the effects of each of four virtual introduced species, the “cyber-invasives”, which represent stylised classes of invaders, propagate through an ecologically realistic network. As an example, we use the data based on the real ecosystem of Flat Holm, a small offshore island.

In the methods, we explain how to construct an ecological network and we demonstrate how the impact of the introduction of an invasive species on each native species can be computed. In the results, we apply this methodology to the Flat Holm network. We find that species in the same taxonomic group respond similarly to invasive species. However, differences within a given taxonomic group can occur. Then, we comment on the predictions given by the model about the perturbations caused by different cyber-invasives. Finally, we discuss the limits of the study, the potential improvements of the mathematical model and the overall results, along with their implications on ecological network analysis.

Methods

Flat Holm Island is a 35 ha offshore island in the Bristol Channel, UK, where very few invasive vertebrate species are present (in particular, there are no invasive non-native rats, Rattus spp.). We use a list of species observed on the island and data from the scientific literature to characterise their diets (supplementary information). We omitted some species because they do not trophically interact with the main food web, such as pollinators and seed-dispersers. This yielded a food web of 227 species: 23 birds, 2 reptiles, 133 invertebrates, 68 plants and 1 fungus (Fig. 1). The system is thus significantly bigger than the the 37-species food web that was previously the largest web studied with this methodology9.

Figure 1
figure 1

The Flat Holm food web. Species are grouped according to their taxonomy: plant, invertebrate, bird, fungus and reptile. Left: The full system. The dots represent the species. Right: The simplified system. Numbers on nodes/links show the number of species/interactions that are aggregated in the formulation of the simplified system.

Full size image

Generalized modelling (GM) is a universal approach to the investigation of dynamical systems8. A generalized food web model which describes the dynamics of ecological networks, in which we consider trophic interactions as a flow of biomass was proposed in8. We define Xi as the biomass density of species i. The dynamics of each species i are described by an ordinary differential equation of the form

$$\frac{{\rm{d}}}{{\rm{dt}}}{X}_{i}={G}_{i}({X}_{i},{T}_{i})+{S}_{i}({X}_{i})-{L}_{i}({X}_{i},{U}_{i})-{M}_{i}({X}_{i}),$$
(1)

where the terms represent the biomass gain by predation G, the gain by primary production S, the biomass loss by predation L, the loss due to other causes of mortality M. Ti and Ui are the total amount of prey and predators that are effectively available to species i. The key insight of generalized modelling is that the local stability of feasible steady states (that is where the species have positive biomass in the states under consideration) in such systems can be investigated highly efficiently without restricting the functions to specific functional forms. The stability of steady states is governed by the system’s Jacobian matrix that contains certain derivatives of Eq. (1).

The Jacobian matrix only contains the partial derivatives of these functions taken at the steady state, e.g. ∂Gi/∂Xi|*, where |* means taken at the steady state. Though these derivatives depend on the unknown steady state, we can assess the stability by normalizing each interaction function by the biomass taken at the steady state, e.g. \({G}_{i}^{\ast }/{X}_{i}^{\ast }\).

$${\frac{{\rm{\partial }}{G}_{i}}{{\rm{\partial }}{X}_{i}}|}_{\ast }={\frac{{G}_{i}^{\ast }}{{X}_{i}^{\ast }}\frac{{\rm{\partial }}{\rm{l}}{\rm{o}}{\rm{g}}{G}_{i}}{{\rm{\partial }}{\rm{l}}{\rm{o}}{\rm{g}}{X}_{i}}|}_{\ast }$$
(2)

The right-hand side of Eq. (2) has two parts, the first matches the scale parameters which define the biomass flows at the steady state (i.e. \({G}_{i}^{\ast }/{X}_{i}^{\ast }\) is the per capita loss rate of species i due to predation at the steady state). The second value, the logarithmic derivative, represents the elasticities or the exponent parameters, which measure the non-linearity of the interaction function at the steady state (for instance a value of 1 represents a linear relationship and 2 a quadratic relationship)8. Using the same notations as in8 for these specific and now interpretable parameters (Table 1), the off-diagonal elements of the Jacobian matrix J can be written as

$${J}_{n,i}={\alpha }_{n}[{\rho }_{n}{\gamma }_{n}{\chi }_{n,i}{{\rm{\Lambda }}}_{n,i}-{\sigma }_{n}({\beta }_{i,n}{\psi }_{i}+\sum _{m=1}^{N}{\beta }_{m,n}{{\rm{\Lambda }}}_{m,i}({\gamma }_{m}-1){\chi }_{m,i})]$$
(3)

and its diagonal elements as

$${J}_{i,i}={\alpha }_{i}[(1-{\rho }_{i}){\varphi }_{i}+{\rho }_{i}({\gamma }_{i}{\chi }_{i,i}{{\rm{\Lambda }}}_{i,i}+{\psi }_{i})-(1-{\sigma }_{i}){\mu }_{i}-{\sigma }_{i}({\beta }_{i,i}{\psi }_{i}+\sum _{m=1}^{N}{\beta }_{m,n}{{\rm{\Lambda }}}_{m,i}[({\gamma }_{m}-1){\chi }_{m,i}+1])]$$
Table 1 Generalized model parameters as defined in8 and21 and the corresponding rules for the parametrization, see supplementary information.
Full size table

Intuitively, the terms in Eq. (3) with ρ correspond to the biomass gain part of the model, i.e. Si and Gi functions (ρi = 0 for primary producers and 1 for predators). Terms with σ are related to the loss in biomass of species i, i.e. Mi and Li functions respectively, either due to natural mortality (μ) or due to predation (σi = 0 if i is a top-predator). The sums over m take into account the feeding competition within and between species, the predation of species m changes when its prey availability is modified.

The size of the model makes it hard to determine all parameters individually. Instead, we use an ensemble approach where different realizations of the model are studied, which are generated by drawing parameters randomly from plausible ranges, identified by biological reasoning. Following7, we choose the values of the parameters such that the system exhibits realistic scaling of biomass turnover and non-linear functional responses and prey switching behaviour. To stabilize the relatively large system to the point where stable steady states occur, we assume quadratic mortality of top predators (supplementary information).

The direct effect of a cyber-invasive j on a native species i is recorded in a perturbation matrix K, defined as

$${K}_{i,j}\,:=\,{\frac{{\rm{\partial }}{A}_{i}}{{\rm{\partial }}{Y}_{j}}|}_{\ast }$$
(4)

where Ai(X) = dXi/dt is the right-hand side of the ODE system and Yj is the biomass of an additional population j7. The perturbation matrix specify the direct effect on the prey of the cyber-invasive but not the relative indirect effects on the biomass of all species. To estimate these indirect effects we compute the impact matrix I, given by

$${I}_{i,j}\,:=\,{\frac{{\rm{\partial }}{X}_{i}^{\ast }}{{\rm{\partial }}{Y}_{j}^{\ast }}|}_{0}$$
(5)

where the derivative is evaluated in the limit of vanishing biomass Yj of the invasive species that enters the system (|0). As long as the perturbation is small, we can make a linearity approximation and write

$${\frac{{\rm{\partial }}{X}_{i}^{\ast }}{{\rm{\partial }}{Y}_{j}^{\ast }}|}_{0}={\frac{{\rm{\partial }}{X}_{i}}{{\rm{\partial }}{A}_{i}}|}_{\ast }{\frac{{\rm{\partial }}{A}_{i}}{{\rm{\partial }}{Y}_{j}}|}_{\ast }$$
(6)

that is

$${\bf{I}}=-{{\bf{J}}}^{-1}{\bf{K}}$$
(7)

where J−1 is the inverse of the Jacobian matrix (this is a corollary of the implicit function theorem10). Although the approximation neglects higher order effects of the non-linearity, basic effects are captured in the Jacobian matrix. The result obtained below and in previous papers9 suggest that even (or perhaps especially) in large systems reliable qualitative predictions can be made in this way.

The perturbation matrix provides information about the direct effect of the perturbation; whereas the impact matrix provides information about the indirect effect of the perturbation, once it has spread through all of the network. The perturbation matrix of the full system (227 species) can be found on line11. For the simpler system, assuming that each cyber-invasive consumes equally all of its prey species, the perturbation matrix is

$$\begin{array}{c}\,\,\,\,\,\,\,\,P\\ {\bf{K}}=\begin{array}{c}Carnivore\\ Herbivore\\ Insectivore\\ Omnivore\end{array}(\begin{array}{r}0\\ -0.1\\ 0\\ -0.1\end{array}\end{array}\,\,\,\begin{array}{r}I\\ 0\\ 0\\ -0.1\\ -0.1\end{array}\,\,\,\begin{array}{r}B\\ -0.1\\ 0\\ 0\\ -0.1\end{array}\,\,\,\begin{array}{r}F\\ 0\\ -0.1\\ 0\\ 0\end{array}\,\,\begin{array}{c}\,\,R\\ \begin{array}{r}-0.1\\ 0\\ 0\\ -0.1\end{array}\,)\end{array}$$
(8)

where −0.1 means a decrease by 10% in the targeted prey biomass. As detailed data on the relative strength of direct impacts is not available we use a generic value to indicate a negative impact.

The impacts are quantified as the indirect relative loss or gain of biomass in the steady state, normalized with respect to the unperturbed steady state (before the arrival of the invasive) per unit of direct effect. These numbers should be read as proxy values with negative number indicating losses and greater absolute values indicating stronger impacts. Thus, invasive species in this model do not cause extinctions, but instead, cause knock-on effects on the biomass of some species in the network: the model tells us the direction and relative magnitude of the perturbation.

We use the generalized model to numerically generate an ensemble of 106 steady states for the Flat Holm food web. For this purpose we draw the generalized model parameters randomly from the identified ranges. The generalized modelling procedure guarantees that all of these states are feasible in the sense that for each state we can write a plausible food web model such all species are stationary and have positive biomass densities. We examine the Jacobian matrices of all steady states in the ensemble to determine their stability and discard all unstable states, which reduces the number of states in the ensemble to ≈9000. For these stable steady states, we then use Eq. (7) to estimate the impact of the four cyber-invasives, where K is chosen to reflect the characteristics of each invasive: we study the effects of a cyber-insectivore, which feeds on all invertebrates (behaving like a hedgehog or a large shrew), a cyber-herbivore, which feeds on all plants (simulating a highly generalist vertebrate herbivore, such as a goat or rabbit), a cyber-carnivore, which feeds on all vertebrates (behaving as a generalist predator such as a cat), and a cyber-omnivore, which feeds on a subset of plants, invertebrates, birds and reptiles. The diet of the cyber-omnivore is based on the rat Rattus rattus and Rattus norvegicus. Rats are very common invaders on islands and known to cause significant impacts across a wide range of taxa. We inferred the likely direct impact of rats from the literature (e.g. Twigg 197512, Atkinson 198513 and Towns et al.14) and from field knowledge of prey taxa likely to appeal to rats (Varnham et al., unpublished data).

For comparison we also used the same methodology to analyse a simpler system where each group of species (birds, plants, etc.) is represented by a single network node (Fig. 1 Right, supplementary information).

Ethical approval

We were not required to complete an ethical assessment prior to conducting our research.

Results

Considering the simple 5-species model first, the method predicts that the directly affected species and their predator are negatively impacted, while the prey and competitors of directly affected species are positively impacted (Fig. 2). These results are consistent with ecological expectations.

Figure 2
figure 2

Impact of the four cyber-invasives on both networks. The y-axis is the positive or negative impact of the invasive species. Each bar is one native species/taxonomic group. Species are grouped in five taxonomic groups: birds (blue), plants (green), reptiles (red), invertebrates (orange) and fungus (purple).

Full size image

Comparing the predictions of the 5-species model with those of the larger network, the general pattern of impact is similar. However, the full model reveals greater detail with some species experiencing an impact that is the opposite than others in the same group. The method thus has the potential to explain patterns observed in the real world, or highlight particular species as sensitive indicators of invasive impact.

In case of the cyber-herbivore the analysis reveals some issues in the underlying data. The model does not resolve the resources for all of the invertebrate species. In the generalized model, these species, therefore, are assigned an abstract carrying capacity which is not affected by the herbivore. The introduction of the herbivore then benefits these species by reducing predation pressure (Fig. 2). We verified that those invertebrate species for which the model predicts a positive impact from the cyber-herbivore are indeed those who were not assigned links to a plant resource. We therefore do not expect these positive impacts to occur in the real world unless an invertebrate indeed utilizes a resource that is not negatively affected by the herbivore.

The value of the detailed model is revealed in the analysis of the cyber-carnivore. While the 5-species model predicts a mild positive impact on plants, the detailed model shows that a few plants are very strongly positively affected while the majority experiences a moderate negative impact (Fig. 2). For the reptiles both models predict only a slight impact by the cyber-carnivore as an increase in prey biomass counteracts direct predation. However, in the 5-species model the net impact is predicted to be positive, while it is negative in the full model.

The cyber-omnivore (or invasive rat) has the most pervasive impact of all the invasives considered, causing large responses throughout the network (Fig. 2). The cyber-omnivore has a strong negative effect on most birds while those birds that are not in the cyber-omnivore’s diet are positively impacted. The reduction in birds and also in reptiles leads to a positive impact on almost all herbivores, which then leads to a negative impact on plants.

To understand the widespread impact of the cyber-omnivore and the cyber-carnivore we also computed the overall sensitivity and influence of species which provides a proxy values for their role in propagating generic perturbations (formal definition in7, supplementary information). A sensitive species is more likely to be strongly perturbed by any disturbances propagating through the network, while an influential species may propagate the perturbation to the other species it is linked to in a more efficient way. Birds and reptiles proved most influential, with reptiles also showing a relatively high sensitivity to perturbations. The generalist bird species, having more links with the rest of the network than the others, are the most influential ones (Fig. 4 in supplementary information).

Discussion

In summary, our results suggest that invasive species can have a wider impact on native ecological networks than previously thought, and that these effects are both direct and indirect. In spite of the many assumptions made for the parametrization and the high-dimensionality of the data, the simplistic model predicts the overall pattern of impacts almost as well as the full model, but loses much of the detail. Moreover, the results highlight the role of birds (and more generally top predators) as important transmitters of indirect effects. The catastrophic impact of invasive species on bird faunas observed in our model is well documented, including the loss of at least five species from Lord Howe Island (Australia) following the introduction of the black rat Rattus rattus13, the loss of at least ten species in Hawaii following the introduction of avian malaria15 and the loss of ten of the 12 species of forest dwelling bird on Guam following the arrival of the brown tree snake Boiga irregularis16. The findings from the field following invasive species introductions in temperate and tropical countries (e.g. New Zealand, UK offshore islands and Antigua) match the model’s outputs for birds, reptiles and mammals13,14,17 (and references therein). Despite being a relatively small component of the overall system, birds play a major role as a transmitter of cascading effects to other species in the network.

There are two main limitations with our approach. First, the methodology constructing the food webs leads to a larger proportion of bird links than other species. This over-representation highlights a common problem with empirically-derived food webs, where some taxa are easier to identify than others and for which more dietary data exist18. Second, a major caveat of the theoretical approach is that it provides a linearised approximation and captures some aspects of non-linearity in the Jacobian matrix, but cannot accommodate higher order effects that occur far from the stationary states and may play a role in the response to strong and sudden perturbations. However this is a greater concern in small simple systems than in the large food web considered here. The analysis used here has the advantage of enabling a much safer analysis than simulation based approaches due to its advantageous numerical and statistical properties. Previous papers, including for instance9, have shown that this approach allows to predict sequences of extinctions based on very limited data. In the absence of better data and tools, the full network method used here presents a reasonable approach to arrive at predictions which take many properties and features of the real system into account.

An obvious improvement to the model considered here would be to incorporate non-trophic interactions. To a certain extent these interactions are taken into account. For instance indirect effects such as exploitative and apparent competition are captured. Direct (non-trophic) competition or mutualism between species of the same trophic group would constitute a competition link in the full model (which is not captured), but appears as a self-limitation term in the simplified model (which we take into account). The good agreement between the results of the full and simplified models therefore suggests that omitting such competition links in the full model should not have a significant effect on the result. Mutualistic interactions, particularly between members of different trophic groups, e.g. pollination and seed dispersal could have a significant effect and we will focus on incorporating them in future versions of the model. However, we feel that these interactions that modelling these interactions naively as biomass flows could lead to unrealistic consequences, and thus more careful modelling work is necessary.

In reality, the impact of alien species on native communities is only likely to increase19. Ecologists need more and better tools if they are to predict their impact and plan their control, and simulation models such as the one presented here could be a part of this toolbox. Conservationists could allocate additional resources to surveying those species which are predicted by this model to be particularly (and unexpectedly) at risk. While pre- and post-eradication surveys are increasingly integrated into eradication projects (e.g.20), these have historically focused on a few taxa–usually highly visible, charismatic and easy to survey species such as birds. While birds are predicted to be particularly susceptible to invasive species, as these models show other species are also likely to be affected. Moreover, it is not just alien predators of vertebrates that can be detrimental, herbivore and insectivore can have significant indirect negative effects. More and better collaboration between theoretical and empirical ecologists, with the questions addressed coming from conservation practitioners, will likely provide the most rapid progress in this field.