Extinction models of robustness for weighted ecological networks

Analysis of ecological networks is a valuable approach to understanding the vulnerability of systems to environmental change. The tolerance of ecological networks to co-extinctions, resulting from sequences of primary extinctions, is a widely-used tool for modelling network ‘robustness’. Previously, these ‘extinction models’ have been developed for and applied to binary networks (sometimes weighted by species abundance) and recently used to predict cascades of co-extinctions in plant-pollinator networks. There is a need for robustness models that can be applied to weighted networks (i.e. where interaction frequencies are recorded) and there is a need to understand how the structure of a network affects its robustness. Here, we developed a framework of extinction models that examine the robustness of networks to random and structurally correlated extinctions (i.e. where avalanches of co-extinctions can occur due to feedback between the two trophic levels). Our models can be applied to networks with binary and weighted interactions. We tested how the average and the range of robustness values is impacted by network structure and the impact of structurally-correlated extinctions sampling non-uniformly from the distribution of random extinction sequences. Our findings are illustrated with plant-pollinator networks. We found that the range of robustness values is driven by the structural heterogeneity of species (plants) in the sequence of primary extinctions. With structurally-correlated extinctions, the networks tested were more robust to ‘avalanche’ extinctions than expected due to the robustness of highly-connected species, but are less robust when preferentially selecting plant species with shared pollinators using a ‘random walk’ model. We found that weighted interactions tend to amplify observed effects and increase variation in robustness. Our framework of models is a new approach to assessing network robustness, permitting it to be calculated with weighted, as well as binary, interaction networks. Models which permit cascades of extinctions vary in their impact on robustness, thus highlighting the vital importance of understanding the model when interpreting robustness.


INTRODUCTION
Network analysis has become an important tool for ecologists seeking to understand the vulnerability of ecosystems to environmental change. Recent research has centred on network approaches for improving our understanding of plant-pollinator communities and extinctions, especially in the light of the widely documented recent declines in key insect pollinators such as honeybees, bumblebees and butterflies (Biesmeijer et al., 2006;Senapathi et al., 2015;Goulson, Lye & Darvill, 2008;Benton, 2006). These trends are concerning for biodiversity, ecosystem function and food security (Potts et al., 2010) as insect pollinators are known to play a vital role in providing ecosystem services (Bailes some pre-determined rule for species survival. Most models, thus far, have used simple rules for secondary extinctions e.g. species become extinct when all their existing links are lost. Network robustness can then be determined from the area under the curve of the proportion of the secondary type that survive against the proportion of the primary type removed (Burgos et al., 2007;see Fig. 1a).
These models have been adapted in various different ways. One approach has been to use specific sequences of primary extinctions, such as ordered by traits of the nodes (Dunne, Williams & Martinez, 2002;Memmott, Waser and Price, 2004;Pocock, Evans & Memmott, 2012;and Santamaria et al., 2016). Another development was to allow rewiring of edges (pollinators switching from one plant to another) based on interactions known from other data , with the aim of increasing the biological realism of extinction simulations -an issue that many of these papers acknowledge is lacking. Most recently, Vieira and Almeida-Neto (2015) allowed feedback in co-extinction between guilds, which they implemented stochastically, so permitting the possibility of cascades of extinctions.
Early extinction models showed that the robustness of communities to random primary extinctions increased with network connectance, the fraction of the possible interactions that were actually observed (Dunne, Williams & Martinez, 2002) and the resulting robustness was often interpreted in terms of network nestedness (Memmott, Waser and Price, 2004). Vieira and Almeida-Neto (2015) also found that cascades were more likely in highly connected networks.
However, more detailed investigation of the impact of network structure on robustness has been lacking. These studies form the foundation of extinction models for plant-pollinator communities on which we base our work.
6 Initially, many of the observed plant-pollinator networks had binary edge weights; interactions between pairs of species were either observed or not.
Increasingly researchers are measuring the frequency or importance of interactions to create weighted networks, yielding a better description of the interactions observed (Ings et al., 2009;Tylianakis et al., 2010) and accounting better for under-sampling biases (Bersier, Banašek-Richter & Cattin, 2002). But the lack of weighted robustness models means that the information on the interaction weights, when known, either has to be discarded (e.g. Pocock et al. 2012) or is used simply to weight the binary outcomes by node abundance . Therefore, there is a need for models that are suitable for measuring robustness on weighted networks. Here, we aim to address this need.
One of the features of these extinction models is that when using random sequences of primary extinctions on a single network, there is a broad distribution in the resulting robustness values (see Fig. 1b). Robustness is therefore a product both of structural heterogeneity of the network (eg Pastor et al., 2012) and of the method of producing extinction sequences; we will explore each of these contributions to robustness in this paper.
Here, our aim was firstly to develop a suite of extinction models which can be applied to weighted, as well as binary, bipartite interaction networks. Secondly, we sought to determine how robustness (assessed using our different models) is affected by network structure and how our models influence measured robustness. Throughout this paper we apply our models to mutualistic bipartite networks (specifically plant-pollinator networks) but discuss how they can be applied to any bipartite network (e.g. with trophic, uni-directional interactions).

MATERIALS AND METHODS
In this study, we examine the robustness of observed plant-pollinator networks that describe observed interactions between species in a community. A network has P plant nodes and A animal nodes, and contains E interactions between species, encoded in the × matrix M. Interactions may be binary (b) or weighted (w).
We illustrate our models and findings using a plant-pollinator network, with data collected by Memmott (1999), from Ashton Court, a site in Bristol, UK. We will refer to this as the Ashton Court (AC) network. This is a well-sampled network (Bluthgen, Menzel & Bluthgen 2006) with interactions recorded over a short period of time (1 month). The AC network is highly resolved: all plants were identified to species (P=25) and many pollinators were identified to species level (and the remaining pollinators identified to morphotype: A=79).
M AC contains 104 species, E=299, with connectance (proportion of realised interactions) of 0.151. Interactions in the AC network are weighted by the number of observations. The degree distribution of plant species is highly skewed as is often the case in plant-pollinator networks.
We also present results for two other networks to illustrate the generality of our model outcomes. These networks were collected in Shelfhanger (Sh), Norfolk, UK (Dicks, Corbet & Pywell , 2002); M Sh has P=16, A=36, E=89 and connectance 0.148; and in Ottowa, Canada (Small, 1976); M Ot has P=13, A=34, E=141 and connectance 0.319. Both networks have interactions weighted by the number of observations. We selected these networks because, like the AC network, they 8 describe northern hemisphere, temperate ecosystems, and have a similar size to the AC network, but differ in having lower and higher connectance respectively.

Model Development
We took as our starting point the extinction model of Memmott, Waser and Price (2004), who analysed the robustness of binary networks by making species of one type (in their case, pollinators) extinct in a random order, i.e. they used a random primary extinction sequence. From this, we developed two new extinction models that each include sub-sequences of plant extinctions that are determined by network structure. All three models (summarized in Fig. 2) can be applied both to weighted and to binary interaction data by introducing a threshold rule.
In this section, we first describe the features that are common to all our extinction models and then outline the distinctive features of each, highlighting the relationships between ours and previous extinction models.

Universal Model Features
Starting from the observed matrix M, a node of one guild (plants) is removed as a primary extinction. Extinctions result in the loss of interactions from M, monitored in the 'reduced' matrix C. The loss of interactions may, according to the rules of the particular model, result in the secondary loss of nodes of the other guild (pollinators). In our new models (see below) the rules admit the possibility of each secondary pollinator extinction giving rise to further knock-on plant extinction(s). These plant extinctions cannot be considered 'primary', but will take their place in what we shall continue to refer to as a 'primary extinction sequence' of the P plant species.
All models proceed until all plant nodes are removed and all species -plants and pollinators -are extinct. The robustness (R) of the network is calculated as where is the number of plant species that have gone extinct (from 0 to P) and ( ) is the number of pollinator species remaining in the network (from A to 0).
R is the normalized (0 < < 1) area under the curve of a graph of the proportion of plant nodes that have gone extinct against the proportion of surviving pollinator nodes (see Fig. 1a). Values of R closer to 1 indicate higher 'robustness' of the network to primary extinctions (Burgos et al., 2007 andAlbert, Jeong, &Barabási, 2000). We use (  The value of R is dependent on the specific sequence of primary extinctions, so running many random extinction sequences will, for all our models, produce a frequency distribution of values of R ( Fig. 1b) which we denote ( ).
A key model feature we adopted is a threshold rule for secondary extinctions, so that a node becomes extinction once it has lost a fraction T or more of its observed interactions (binary M), or of its observed total interaction weight (weighted M). Clearly the value of T that we choose is arbitrary. It must lie in the range 0 < ≤ 1 . [T=0 is a pathological case; all pollinator species become extinct after the first primary plant extinction; T = 1 generates the extinction rule for most previous models (Dunne, Williams & Martinez, 2002;Memmott, Waser & Price, 2004;Kaiser-Bunbury et al., 2010), although Vieira introduced a nodespecific threshold ≤1.] We generated distributions of robustness ( ) for a range of threshold values (0.1 to 1 in 0.1 intervals) for three observed plant-pollinator networks (Ashton Court, Shelfanger and Ottawa) to determine the effect of T on the results of our extinction models. We then chose a threshold of T = 0.5 for the remainder of the paper: i.e. a secondary extinction occurs when a node has lost at least half of its interactions (binary M) or weights (weighted M). It should be noted that the 'effective T' (Teff ) could be greater than T; for example, with a binary network and T = 0.5, a node linked to 5 others would go extinct after losing 3 edges, giving an effective T of 3/5=0.6. Since most pollinators are observed visiting a relatively small number of plants, the difference between the specified and the 'effective' threshold can be noticeable and so we calculate the node-averaged Teff in all cases.

New extinction model features
We present three distinct models, which we denote: 1. Secondary Only (SO), 2.
Deterministic Avalanche (DA) and 3. Random Walk (RW). Each model can be used with binary or weighted interaction data and is prefixed with 'b' or 'w' to indicate which.

Model 1. Secondary Only model (bSO and wSO)
In the Secondary Only model the order of primary plant extinctions is random.
All pollinator extinctions are secondary and determined by the threshold rule.

11
There is no spread of extinctions beyond the secondary extinction of pollinators.
The method is as follows: 1 Then calculate R according to equation 1.
Were = 1 used here, step 2c would never result in tertiary plant extinctions and no avalanches would occur, so the DA and SO models would be identical. The 'stochastic co-extinction model' developed by Vieira and Almeida-Neto (2015) is a special case of our bDA model where the threshold is applied stochastically and is node specific; specifically, extinctions of nodes at our step 2c occur with probability = 1-(remaining interactions) / (interactions at start). We adopt the term 'avalanche' for our spreading deterministic extinctions to differentiate them from the stochastic 'cascades' of Vieira and Almeida-Neto (2015), which occur once only, triggered by the first primary extinction.

Model 3. Random Walk model (bRW and wRW)
The RW model is similar to DA, in that a trigger can cause an avalanche of non- 6. Identify plant f as the new e and make it extinct 7. Loop through steps 2 to 6. If no neighbours exist in step 3, revert to step 1.
Repeat steps 1 to 7 until there are no species remaining in the network. Then calculate R according to equation 1. 14

Natural extensions of our models
We have developed these three models for application to binary and weighted mutualistic bipartite networks and with a random order of primary plant extinctions (i.e. the selection of the next extinction in step 1 of Models 1,2 and 3 is random). However, we note that these models can easily be modified to use ordered primary extinctions, where the choice of plant in step 1 of Models 1,2 and 3 is according to a pre-determined rule (based on node degree, biological plant trait etc). The models can also be applied to bipartite networks with unidirectional dependencies (no feedback between the trophic levels, e.g. trophic or host-parasitoid interactions), though in that case avalanches cannot occur.

Comparison of robustness distributions from the three extinction models
The distribution ( ) generated from a single network M will depend on the model used and whether the data are weighted or binary. If there are P plant species in the network, there are P! distinct plant sequences. The SO models sample uniformly from these possibilities (i.e. all sequences are equally likely).
The DA and RW models do not sample uniformly, because avalanches produce non-random sub-sequences determined by the structure of the network. Using binary and weighted versions of the Ashton Court (AC) network we generated 25,000 extinction sequences using each of the 3 models, in order to assess the effect of R on model choice. To create values of R that lie close to the theoretical maximum and minimum we ran bSO with plant extinctions in order of increasing and decreasing degree.

Assessing node and network-level determinants of variation in robustness
Having described the variation in robustness, we finally sought to assess the attributes which determine this variation under each model. We did this in two ways: by creating null versions of the AC network with defined characteristics, and by examining extinction rank (when a plant goes extinct in an extinction sequence). These calculations are designed to test the role of plant degree in determining the central tendency and spread of ( ) in our three extinction models.
Generating null networks to test effect of degree distribution To explore the effect of degree distribution ( ) on robustness we created three The bSO, bDA and bRW models were run on the observed AC network and on these three exemplar random networks to produce 10,000 extinction sequences, and a corresponding distribution of robustness ( ) for each.

Examining plant extinction rank in relation to R
To explore whether (for example) high-degree plants tend to go extinct toward the beginning of a primary extinction sequence, we defined the position in a sequence when a plant became extinct as its extinction rank (r), 1 ≤ ≤ . We ran each extinction model 25000 times, using binary and weighted versions of the AC dataset, and compared ℎ( ), the distribution of extinction rank for each species generated by the simulations. From h(r) for each species we calculated the median extinction rank (rm) and degree (k) and tested for correlation using the Spearman coefficient. We expected that, by definition, r would be equal across all species for the SO model but not for the DA or RW models, since avalanches and random walks will tend to select (or avoid) high-degree nodes preferentially.

Testing on other networks
We tested our models on the Shelfanger and Ottawa networks. For each network we generated 25,000 extinction sequences, using each of the 3 models, in binary and weighted form. We used a fixed threshold of = 0.5 for all cases as we are not directly comparing the networks, only seeking to confirm the generalities of the resulting ( ) distributions.

Varying the value of the threshold for secondary extinctions
Median robustness Rm increases non-linearly with T, and the least robust of our three networks at low T becomes the most robust at high T (Fig. 3a). However, this appears to be an artefact of the relationship between T and Teff, the nodeaveraged effective threshold (Fig. 3b), because Rm increases linearly with Teff and the three networks are increasingly robust in order of increased connectance, as found by Dunne, Williams and Martinez (2002), at all values of Teff (Fig 3c). The remainder of our results are presented for the AC network only (where Teff = 0.694 for our chosen T =0.5).

Robustness Distributions
The distributions ( ) produced by each of the 3 models for binary and weighted data (Fig. 4)  weighted, not binary, interaction data increases the IQR of ( ) in all cases.

Network randomisation test
Compared to the results of the binary extinction models for the AC network (Fig.   5a), we found that randomizing the degree distributions caused the distribution of robustness ( ) to be narrower (Fig. 5b-d), and this was especially so when the plant degree distribution is randomised (5c and 5d). This confirms that the observed, highly skewed, plant degree distribution of the AC network drives the broad robustness distributions we generate for this network. Note though that the Rm remain in the same order (RW<SO<DA) in every case, showing the consistency of effect from these models.

Extinction rank of plant species, and the effect on R
Plant degree is a predictor of the plant's extinction rank in the DA and RW models. In the SO models, the expected rank is constant for all plant species, irrespective of degree, because the extinction sequence is entirely random. In contrast, the observed extinction ranks of two example plant species from the DA and RW models are clearly skewed (shown in the insets of Fig. 6). In the DA models median extinction rank is positively correlated with plant degree ( robustness is high compared to the SO models (Fig. 4b cf. Fig. 4a). In contrast, in the RW models plants with high degree are more vulnerable to extinction (the model preferentially 'homes in' on well-connected plants) which results in overall lower robustness of the network (Fig. 4c cf Fig. 4a). Clearly the apparent robustness of a network is dependent on the model used to incorporate feedback and structurally correlated co-extinctions.

DISCUSSION
The robustness of an ecological network is a valuable metric, because from it we can learn several things and it provides a quantitative metric for describing the vulnerability of an ecological community to hypothetical extinction scenarios, i.e.
sequences of primary extinctions. Based on this we can compare the vulnerability of different networks to these scenarios and determine properties of the environment or of the network that correlate to robustness. However, robustness can also be used to say something about the network itself, in which case the breadth of a robustness distribution is important.
Extinction models that calculate robustness have been around for a long time but until recently ecologists have used simple versions of these models (see Fig. 1).
Here, we have created a framework of extinction models, building on those of Importantly we have used an extinction threshold of less than one, i.e.
pollinators can go extinct before all their plants go extinct. This addition has an ecological motivation -pollinators may decline to extinction due to declining food resources -and adds greatly to the flexibility of the model. Having T<1 20 allows us to create weighted versions of our models and permits the potential for feedback and, hence, avalanches of extinctions cascading between the trophic levels (as first shown by Vieira and Almedia-Neto, 2015).
All of our extinction models, in binary and weighted form produce a broad distribution of robustness values ( ) for every network we analysed, indicating that there is some aspect of the structure of the network which causes this variation. We found this to be the degree distribution of the nodes, specifically the plant degree distribution. Plant pollinator networks tend to have fewer plant species than pollinator species ( < ), so the potential for a skewed plant degree distribution in increased, making it more influential on robustness, as indeed was the case in our test network (Memmott 1999).
Though 'robustness' has in the past been used to suggest priorities for conservation or management (Pocock, Memmott andEvans. 2012, Devoto et al. 2012), extinction models are not an attempt to predict precisely how an ecosystem would collapse. They do, nonetheless, offer a means to quantify and compare the structure of ecological networks, although to do this we need to ensure we are comparing like-for-like. We found that the median robustness Rm is a linear function of (effective) threshold Teff and the effect is consistent across networks (Fig. 3). However, Teff has a non-linear relationship with the selected threshold value T, so if our models are to be used in the future to compare the robustness of ecological networks, it is important to note that Teff, not T, is the appropriate parameter to ensure comparability.
Plant-pollinator communities are increasingly described with weighted interactions and so our weighted extinction models for robustness are a valuable advance. We found (Fig. 4)  There are different ways in which feedback between trophic levels can be applied and we developed two illustrative models: the Deterministic Avalanche (DA) and the Random Walk (RW) models. These models (and others like the cascade model developed by Vieira and Almedia-Neto, 2005) may appear to be generating new outcomes but in reality, they simply produce a non-random sample of robustness values from those generated by a simple SO model. The Ashton Court dataset generated a huge range of R values, all of which can be realised in the Secondary Only models. The DA and RW models preferentially sample extinction sequences to produce skewed subsets of the SO outcomes (the P! extinction sequences are not all equally likely, and some will be impossible).
The Deterministic Avalanche Model preferentially samples nodes that are 1 step away from each other in the network and extinctions can 'fan out' from each 22 trigger (a randomly selected plant extinction). This mechanism can be likened to a 'breadth first' search on the network (see for example, Kolaczyk, 2009). On the other hand, the Random Walk model follows a path though the plant projection network F away from trigger in the manner of a "depth first" search (Kolaczyk, 2009 Overall, the SO, DA and RW models can be used to examine the effect of a range of ecological extinction scenarios. There are many possible applications of the SO model. The DA and RW models are more specific in their application and notably produce opposite effects on the distribution of robustness for a given network.
Therefore, this highlights the importance of selecting an appropriate model for the ecological circumstances. We expect that there will be future developments to refine these models to usefully predict the outcomes of possible ecological scenarios.
All of these extinction models are designed to be applied to real network data.
Therefore, it is vital to consider the quality and reliability of the data being used.
Pollination networks vary hugely in sampling method, period of collection and taxonomic resolution. We caution against comparing the outcomes of extinction models across multiple networks that have come from different sources without consideration of the data, methods and specific question. For example, CaraDonna et al. (2017) highlight the potential pitfalls of assuming that a static network, comprising all accumulated interaction data, is an appropriate representation of a community. Further work in understanding temporal variation and the description of fully-resolved plant-pollinator networks are key to improving the utility of extinction models.
We hope that by improving our understanding of extinction models at a fundamental network level, and by setting out different areas of model extension, our work will guide future developments in the analysis of the vulnerability of ecosystems to environmental change.