Ranking of critical species to preserve the functionality of mutualistic networks using the k-core decomposition

Data

We have analyzed the Web of life collection (Fortuna, Ortega & Bascompte, 2014), comprised by 89 mutualistic networks, with 59 communities of plants and pollinators and 30 of seed dispersers (http://www.web-of-life.es/). There are 57 communities with binary adjacency matrix (i.e., the interaction between the two species is recorded but not its strength), and 32 with weighted matrix, where the strength is accounted for. Network sizes range from 6 to 997 species, the minimum number of links is 6 and the maximum is 2,933.

Decomposition and k-magnitudes

The idea of core decomposition was first described by Seidman (1983) to measure local density and cohesion in social graphs. It has been successfully applied to visualize large systems and networks (Alvarez-Hamelin et al., 2005; Kitsak et al., 2010; Zhang et al., 2010; Barberá et al., 2015).

The k-core of a network is a maximal connected sub-network of degree greater or equal than k. That means that each node in the sub-network is tied to at least k other nodes in the same sub-network.

A simple algorithm to perform the k-core decomposition prunes links of nodes of degree equal or less than k (Batagelj & Zaversnik, 2003). The process starts removing links with one of their edges in a node of degree 1. This procedure is recursive and ends when all the remaining nodes have at least two links. The isolated nodes are the 1-shell. Then it continues with k = 2, and so on. After performing the k-decomposition, each species belongs to one of the k-shells (Fig. 1). The m-core includes all nodes of m-shell, m + 1-shell... with m ranging from 1 to the max k index of that particular network. For instance, the 2-core of Fig. 1 is the union of the 1- shell and the 2-shell.

Mutualistic networks are bipartite, with two guilds of species (plant–pollinator or plant-seed disperser in the studied collection). Links among nodes of the same class are forbidden. We will call these guilds A and B.

Based on the k-core decomposition, we define three k-magnitudes. In order to quantify the distance from a node to the innermost shell of the partner guild, we define k_radius. The k_radius of node m of guild A is the average distance to all species of the innermost shell of guild B. We call this set N^B. (1)${\displaystyle k_{\mathrm{radius}}^A\left(m\right)=\frac{1}{\mid N^B\mid }\sum _{j\in N^B}dis{t}_{mj}m\in A,\forall j\in B}$ where dist_mj is the shortest path from species m to each of the j species that belong to N^B. The minimum possible k_radius value is 1 for one node of the innermost shell directly linked to each one of the innermost shell set of the opposite guild.

To obtain a measure of centrality in this k-shell based decomposition, we define k_degree as (2)${\displaystyle k_{\mathrm{degree}}^A\left(m\right)=\sum _j\frac{a_{mj}}{k_{\mathrm{radius}}^B\left(j\right)}m\in A,\forall j\in B}$ where a_mj is the element of the interaction matrix that represents the link, considered as binary. If the network is weighted, a_mj will count as 1 for this purpose if there is interaction, 0 otherwise. It could be understood that k_degree(m) is a fine-grained measure, whereas degree would be a coarse-grained measure. k_degree(m) is like a “continuous” degree where each node i linked to node m adds the inverse of its k_radius(j). Generalists score high k_degree, whereas specialists, which have only one or two links, with similar k_radius, score lower k_degree. This magnitude reminds the definition of the Harary index (Plavšić et al., 1993) but only considering paths from the nodes tied from m to the nodes of the innermost shell.

Figure 2 shows how k_degree works for one particular network. There are many nodes with the same degree value (Fig. 2A), such as specialists with just one or two links, that from a ranking point of view are equivalent. On the contrary k_degree, maps the degree distribution onto a more continuous one (Fig. 2B), because of the contribution of the inverse of k_radius. In Fig. 2C the cumulative distributions of both indexes are overimposed over the degree scale (We perform the linear regression k_degree = α × degree and fit the k_degree distribution to the degree values).

Finally, we introduce k_risk as a way to measure how vulnerable is a network to the loss of a particular species: (3)${\displaystyle k_{\mathrm{risk}}^A\left(m\right)=\sum _j{a}_{mj}\left(k_{\mathrm{shell}}^A\left(m\right)-k_{\mathrm{shell}}^B\left(j\right)\right)+ϵk_{\mathrm{shell}}^A\left(m\right)m\in A,\forall j\in B,k_{\mathrm{shell}}^B\left(j\right)<k_{\mathrm{shell}}^A\left(m\right).}$

The k_risk of a given species is the sum of the k-shell differences to all the nodes of the other guild on a lower k-shells to which it is connected. Each one is weighted by the difference of the k indexes. The second element of Eq. (3) is meant to solve ties among species when they belong to different k-shells, and is a very small quantity (in our implementation we use 0.01, two orders of magnitude lesser than the sum).

In an intuitive way, if we remove one node strongly connected to others of lower k-shells, these species are in high risk of being dragged by the primary extinction. On the other hand, the extinction is much less dangerous for the species of higher k-shells linked to the same node, because they enjoy more redundant paths towards the network nucleus.

Applying the k-magnitudes to a network

Figure 3 is an small seed disperser network with five species of plants, four species of thrushes and eleven links. We call, by convention, guild A the set of plants, and guild B the set of birds. The k-core decomposition was performed with the R igraph package (Csardi & Nepusz, 2006). The maximum k index is 2. The four bird species belong to 2-shell; there are three plant species in 1-shell and two in 2-shell. In this example each species of 2-shell is directly tied to all species of the opposite guild 2-shell, but this is not a general rule.

The shortest path from plant species 2 to each of the four bird species of 2-shell is 1, because of the direct links. So, $k_{\mathrm{radius}}^A\left(2\right)$ is 1. The same reasoning is valid for plant species 1. The reader may check that the k_radius of bird species of 2-shell is 1 as well, measuring their shortest paths to plants species 1 and 2.

Computation of this magnitude is simple although a bit more laborious for 1-shell plant species. We work plant species 4 as an example. First, we find the shortest paths to each bird species of 2-shell. Shortest paths are depicted with different colors. Plant species 4 is tied to seed disperser species 1, so distance is 1. On the other hand, there is no direct link with bird species 2. Shortest path is pl4-disp1-pl2-disp2, and distance is 3. It is easy to check that distances from plant species 4 to bird species 3 and 4 are also 3. Once we have found the four distances, we compute $k_{\mathrm{radius}}^B\left(4\right)$ as the average of 1, 3, 3 and 3, that is 2.5.

The values of k_degree are straightforward to compute. For instance, the k_degree of disperser species 1 is: (4)${\displaystyle k_{\mathrm{degree}}^B\left(1\right)=\frac{1}{k_{\mathrm{radius}}^A\left(1\right)}+\frac{1}{k_{\mathrm{radius}}^A\left(2\right)}+\frac{1}{k_{\mathrm{radius}}^A\left(4\right)}+\frac{1}{k_{\mathrm{radius}}^A\left(5\right)}=2.8.}$

The last k-magnitude we defined was k_risk. We use again the disperser species 1 as example. Links to species of the same or upper k-shells are irrelevant to compute k_risk, so only plant species 4 and 5 are taken into account. (5)${\displaystyle k_{\mathrm{risk}}^B\left(1\right)=k_{\mathrm{shell}}^B\left(1\right)-k_{\mathrm{shell}}^A\left(4\right)+k_{\mathrm{shell}}^B\left(1\right)-k_{\mathrm{shell}}^A\left(5\right)+ϵk_{\mathrm{shell}}^B\left(1\right)=\left(2-1\right)+\left(2-1\right)+0.01x2=2.02.}$

This magnitude may seem counter-intuitive, because the k_risk of a highly connected species like plant 1 is 0.02, almost the same of that of peripheral plant 3 (0.01). This is because plant 1 has no ties with lower k-shell animal species. The k_risk ranks species to assess resilience, it has not an absolute meaning. It just tells us that it is more dangerous for the network to remove the disperser 1 than plant 1, and plant 1 than plant 3 (Table 1).

The k-magnitudes of the example network are shown in Table 1.

Extinction procedures

In order to rank the critical species to preserve the functionality of mutualistic networks and visualize eventual decomposition of the giant component (i.e., the highest connected component of a given network), we carried out two static extinction procedures. Static assumption implies that there is not rewiring (e.g., plants that have lost their pollinators are not pollinated by other insects), despite this kind of network reorganization is observed in nature (Ramos-Jiliberto et al., 2012; Goldstein & Zych, 2016; Timóteo et al., 2016). Nodes are ranked once, before the procedure starts, as in most of robustness assessments studies (Memmott, Waser & Price, 2004; Kaiser-Bunbury et al., 2010; Domínguez-García & Muñoz, 2015).

In the first method, one species is removed each step, in decreasing order according to the chosen index, no matter to which guild it belongs. Four ranking indexes are compared: k_risk, k_degree, degree and eigenvector centrality. The k indexes were computed with the R package kcorebip; degree and eigenvector centrality with the degree and evcent functions of the igraph package.

To estimate the damage caused to the network, the fraction of remaining giant component was used. The procedure stops when this ratio is equal or less than 0.5. To break ties, we ran 100 experiments for each network and index, shuffling species with the same ranking value. The percentage of removed species needed to get to 0.5 of the remaining giant component is used to measure the performance of the ranking. The lower the percentage of removed species, the more efficient the ranking is in destroying the network. The top performer scores the least average removal percentage (Fig. 4).

The second extinction procedure that we followed is more common in the literature. Only animal species are actively removed (primary extinctions); secondary extinctions happen when nodes become isolated (Memmott, Waser & Price, 2004).

The fraction of surviving plant species is measured as a function of the removed fraction of animal species (Figs. 5A and 5C) and the area under the curve is the value to compare performance. We averaged the results of 100 repetitions.

In this case, in addition to the four indexes of the first experiment, we include MusRank a non-linear ranking algorithm for bipartite networks (Tacchella et al., 2012), inspired by PageRank (Allesina & Pascual, 2009). This algorithm is not valid for the first extinction method. Domínguez-García & Muñoz (2015) showed that MusRank achieves excellent performance for this extinction procedure.

In the second extinction procedure, we also measured the fraction of remaining giant component (Figs. 5B and 5D). Extinction sequences are identical, the only difference is that both magnitudes are measured for each step.

First extinction method

k_risk was the ranking method with the lowest average species removal percentage to destroy half of the Giant Component in most of the networks (67 out of 89 networks) (see supplemental material, Table S1). Figure 6 shows the performance comparison of the four ranking criteria. There are some ties, more frequent when networks are small. Network size is the key factor to explain why the performance range is so wide.

As size increases, the removal percentage to break the giant component decreases (Fig. 7). When the network is big, the primary extinction of key nodes triggers an important amount of secondary ones. If the community has 100 or more species, k_risk is even a better predictor of the most damaging extinction sequence and outperforms the other indexes for 28 out of 32 networks.

Second extinction method

MusRank ranking method had the lowest area under the extinction curve for 85 of the 89 studied networks (Fig. 8), and in the other 4 the difference is so small that may be just an effect of the averaging procedure. So, MusRank is the optimal ranking index to destroy the network following this algorithm.

On the contrary, when the efficiency of the network destruction was measured through the area under the curve (AUC) of the surviving Giant Component fraction the MusRank index had the highest values, placing it as the least efficient ranking method according to this criterion (Fig. 9). In this case, k_degree is the most efficient index for 42 out of 89 networks. We must underline that the extinction sequences are the same, the only difference is the measured output.

We have worked out one example (Fig. 5) to explain this shocking difference in performance depending on the measured outcome. On the upper row (subplots A and B), the difference for both ranking indexes when measuring the giant component. While this magnitude decreases at a constant pace for MusRank, there is a sharp reduction of the component size when one third of animal species are removed following the k_degree ranking. On the lower row (subplots C and D), opposite results are obtained when accounting for the fraction of surviving plant species.

The destruction of this pollinator network sheds light on the root cause of the difference. The network has 36 pollinator and 16 plant species (Fig. 10A), two of them are outside the giant component. When the 13 top animal species ranked by MusRank are removed (pollinators 3, 1, 7, 15, 32, 6, 14, 33, 13, 31, 8, 16, 10), the community reaches the degraded structure of Fig. 10B. The size of the giant component is 27 (54% of the original), and there are 23 pollinator and 6 plant species.

If we remove the 13 top animal species ranked by k_degree (pollinators 1, 3, 7, 13, 15, 2, 11, 20, 12, 8, 6, 5, 10) instead, the community structure is that of Fig. 10C. Now, the size of the giant component is 19 (38% of the original), and there are 23 pollinator and 9 plant species. MusRank has killed more plant species, but the giant component is clearly smaller ranking by k_degree.