Extinction-induced community reorganization in bipartite networks

We study how the community structure of bipartite mutualistic networks changes in a dynamic context. First, we consider a real mutualistic network and introduce extinction events according to several scenarios. We model extinctions as node or interaction removals. For node removal, we consider random, directed and sequential extinctions; for interaction removal, we consider random extinctions. The bipartite network reorganizes showing an increase of the effective modularity and a fast decrease of the persistence of the species in the original communities with increasing number of extinction events. Second, we compare extinctions in a real mutualistic network with the growth of a bipartite network model. The modularity reaches a stationary value and nodes remain in the same community after joining the network. Our results show that perturbations and disruptive events affect the connectivity pattern of mutualistic networks at the mesoscale level. The increase of the effective modularity observed in some scenarios could provide some protection to the remaining ecosystem.

We study how the community structure of bipartite mutualistic networks changes in a dynamic context. First, we consider a real mutualistic network and introduce extinction events according to several scenarios. We model extinctions as node or interaction removals. For node removal, we consider random, directed and sequential extinctions; for interaction removal, we consider random extinctions. The bipartite network reorganizes showing an increase of the effective modularity and a fast decrease of the persistence of the species in the original communities with increasing number of extinction events. Second, we compare extinctions in a real mutualistic network with the growth of a bipartite network model. The modularity reaches a stationary value and nodes remain in the same community after joining the network. Our results show that perturbations and disruptive events affect the connectivity pattern of mutualistic networks at the mesoscale level. The increase of the effective modularity observed in some scenarios could provide some protection to the remaining ecosystem.

I. INTRODUCTION
Mutualistic interactions between species are often represented as bipartite networks, where the interactions occur between two groups of species (generically resources and consumers), but not within the groups. Empirical mutualistic networks exhibit a number of macro-scale structural features such as nestedness [1,2], where specialists interact with proper subsets of the species that generalists interact with; modular organization, that captures the block structure [3,4]; and stability, which can be measured as the largest eigenvalue of the appropriate matrix [5]. Biological systems, and in general any complex system, are expected to withstand the loss of elements, either by random failure or driven by a directed perturbation (e.g., environmental change or a targeted attack) [6,7]. In the context of ecology, loss of biodiversity as a consequence of environmental perturbations disrupts ecosystems and their functioning significantly. The emergence of modularity is crucial for community ecology because such a compartmentalized structure can greatly influence dynamics, as the compartments buffer the spread of perturbation across the network [8,9].
Here, we analyze how species extinction affects the structure of mutualistic networks. Our study focuses on consumer removals, as, as pollinators have a higher immediate extinction risk than plants and also loss of a pollinator species may cause the co-extinction of plants that depend on them [10][11][12]. We present a detailed analysis of the community structure in response to the loss of pollinator species, using an empirical mutualistic network.
The study of the changes in community structure under different extinction scenario sheds light on the fragility of ecological communities to species extinctions. We further emphasize our results by showing how in a model * Author for correspondence: somaye@ifisc.uib-csic.es of mutualistic network growth [13] modularity and nestedness remain basically unchanged, in contrast to our results when extinction mechanisms are at play in a real bipartite network.

II. SPECIES EXTINCTION IN EMPIRICAL BIPARTITE NETWORKS
We analyzed a plant-pollinator interaction network, sampled in Mallorca (Balearic Islands). The dataset was collected from a dune marshland located at sea level in the northeast of the island (Son Bosc; SB hereafter). The authors of [14] sampled insect-flower visitation events during the consecutive flowering season, from April to July on randomly selected flowering plants. A total of 696 flower visits between 80 plants and 162 pollinators were recorded (Fig. 1). networks. These are composed of two different kinds of nodes: resources, here the n plants, and consumers, the m insect pollinators. The interactions between resources and consumers are represented by the incidence matrix A(n × m), whose entries A ij are equal to 1 if there is a mutualistic relation between nodes i and j, and A ij = 0 otherwise. In this work, we only consider the existence of an interaction but not its weight.
At each time step, an extinction event is modeled by removing a pollinator for node removal scenarios or an interaction for the interaction removal scenario. If a node loses all of its links, it becomes extinct. We simulate the loss of nodes and links with four different scenarios: (1) uniformly at random [15]; (2) directed extinction, in which the removed pollinator is chosen with a probability proportional to her number of links (degree) [10,16,17]; (3) generalist scenario, in which pollinators are sequentially removed from the most to the least connected pollinator (in case of a draw one of them is chosen at random); (4) specialist scenario, in which nodes were sequentially removed from the least-degree pollinator to the most-degree [18]; and (5) random interaction extinctions, in which links were removed randomly to model the disappearance of an interaction. To balance for the different number of nodes and interactions, we measure time as the fraction of nodes removed or the fraction of links removed. That is, for scenarios where nodes are removed, each event represents a time step of 1/m, while for interaction removal, a time step corresponds to 1/l, where l is the number of interactions. When all the nodes or all the links are removed, the time is equal to 1.

a.
b.
c. d. In order to compare visually the different extinction scenarios, we plot in Fig. 2 alluvial maps showing the community structure for a single realization of the differ-ent extinction dynamics. We have used the open-source library BiMat [19] to compute modularity and community structure in bipartite networks. We show the resulting community structure for the same fraction of time as described above for all the scenarios. For random extinction and especially for specialist extinctions, the communities remain similar during the initial steps. This is in contrast with directed extinction and random interaction extinction, where the communities change much more from the initial condition. The dynamics of the community structure in the scenario of directed extinctions behaves similarly as in the generalist extinction scenario (not shown). As expected, networks also lost more nodes in every level of directed extinction.
Note also that for a low fraction of extinction events directed extinction and interactions extinction perform very similarly. This is related to the fact that choosing an edge at random is similar to selecting nodes proportional to their degree.
We quantify these observations by measuring the modularity (Q) of detected communities, that is, densely connected non-overlapping subsets of nodes. The modularity of a bipartite network given a partition is defined as [20]: where A ij is the incidence matrix of the network, p ij is the null model matrix describing the expected probability of interactions between two types of nodes given their degrees, k i is the degree of resource node i and d j the degree of consumer node j; g i and h i are the community indices of nodes i and j and |E| is the number of links in the network. After a certain number of extinctions, eventually, bipartite networks break in a set of disconnected components. Thus, to consider the breakup of the network, we introduce the effective modularity Q e , which is calculated as the product of the relative size of the largest connected component S and the modularity Q: Q e = SQ. In Fig. 3 we show the effective modularity for the different extinction scenarios averaged over 500 independent realizations of the dynamics. We observe two behaviors: On the one hand a transition-like behavior, where there exists a critical fraction of extinction events for which the effective modularity sharply decreases, as happens for the generalist, directed and interaction scenarios, with critical fractions of extinction events approximately 0.35, 0.6 and 0.9 respectively. For the generalist scenario, the effective modularity even increases initially. On the other hand, in the random and specialist scenarios, the effective modularity decreases smoothly until the bipartite network is extinct.
In order to gain more insight into the structural reorganization of the network we measure several other garding the size of the largest component ( Fig. 4.a), which monitors the fragmentation of the network, we observe that the network is more sensitive to generalist species extinctions and collapses faster than the other scenarios, followed closely by the directed extinctions scenario, in line with previous knowledge on robustness under targeted attacks in networks [15,[21][22][23][24][25]. Random node extinctions and specialist extinctions behave very similarly, with a smooth decay of the largest connected component. Last, random interaction extinctions keep the largest connected component larger than in any other scenario until a fraction of extinctions equal to 0.8, where the system rapidly collapses. The number of communities ( Fig. 4.b) is increased initially in all the scenarios, but then remains constant around 10 communities for random node extinctions and specialist extinctions, decaying fast when the fraction of extinction events is almost 1 due to network decomposition. For random interaction extinctions, it behaves similarly in the beginning but toward a fraction of extinctions of 0.6, the number of communities increases rapidly to up to 20 before dropping fast again because of network decomposition. Finally directed and generalist extinctions behave similarly, with a steady increase in the number of communities until a certain fraction of extinction events, where the number of communities decays. This certain fraction is 0.55 and 0.7 for the generalist and directed scenarios respectively and they reach 33 and 21 communities respectively.
The next architectural pattern that we consider here is nestedness, which can be described as the tendency of specialists to interact with proper subsets of the nodes interacting with generalists [26,27]. There are several indices for quantifying nestedness depending on whether binary or weighted interaction data are provided. The most commonly used methods are: NTC (Nestedness temperature calculator) [28], SR (spectral radius of the adjacency matrix) [29] and NODF (Nestedness metric based on overlap and decreasing fill) [30]. Here we use the NODF metric to estimate nestedness.
In all the extinction scenarios, nestedness values decrease with extinction events from the very beginning (Fig. 4.c), due to the decrease of the largest degree of the bipartite network, which is positively correlated with NODF [31]. Therefore, the fastest decrease is found for generalist extinctions, followed by directed extinctions. Random node and interaction extinctions behave similarly, decaying more smoothly, almost in a linear fashion, to reach 0 nestedness when t = 1. Last, the specialist extinctions scenario is the one keeping the network more nested, related to the fact that this scenario is the one diminishing the largest degree the less.
The structural changes in the community structure of the bipartite network can be quantified with the community persistence, i.e., the probability that two nodes remain in the same community if they were initially in the same community, P i,j (M i = M j , t|M i = M j , t 0 ). We then compute the averages over all node pairs to get the mean persistence P M = P (M i = M j , t|M i = M j , t 0 ) . As illustrated in Fig. 4.d, community persistence decays initially fast and as more nodes (or links) are extinct for any scenario. In the random extinction scenario, the persistence decays at a slower rate. In random interaction extinctions, the persistence decays quickly after the extinction of a small fraction of interactions and then still decreases at a lower rate until the extinction of around 90% of the interactions. We observe an increase of the persistence in the directed and generalist scenarios. This increase is due to the breakup of the bipartite networks where the few interactions remaining corresponds to interactions that originally were identified in the same community.
The variability of community structure is captured with the versatility, V . Versatility is a metric of nodal affiliation which describes how closely each node is assigned with a community: V = 0 indicates that a node is always assigned to the same community; while V 0 determines that it is assigned to different communities depending on the realization [32]. The versatility of a node j is defined as: with a(i, j) = 1 if i and j are in the same community. 0 otherwise.
where a(i, j) is the expected value of a(i, j) averaged over different realizations evolved to the same fraction of extinction events. A high value of versatility reflects thus a loose community structure, as nodes might be assigned to one or other community, while a low versatility value stands for a robust community structure with well-defined communities. For versatility (Fig. 5) the results show that The versatility defined in Eq. 2 averaged over 500 realizations is represented in a color scale from V = 0 (red) to V = 3 (blue), the x-axis is the fraction of extinction events, and the y-axis represents each species. Each column corresponds the five extinction scenarios: (from left to right) random extinction, directed extinction, generalist extinction, specialist extinction and random interaction extinction. (Bottom) versatility averaged over all species and 500 realizations of the extinction sequences for the different scenarios.
random node and interaction extinctions behave similarly, with a decreasingly less defined community structure up to 75-80% of extinction events (growing versatility), followed by a decrease in versatility, associated with a more solid community structure. For directed extinctions, the structure evolves rapidly to a not well-defined community structure (high versatility) and around 25% of extinction events starts to build a more solid community structure, as versatility decays. For the generalist and specialist extinction scenarios, the picture is a bit more complex. Due to the semi-deterministic nature of the extinction sequences in these scenarios, at some points all of the realizations reach the same configuration, and thus the same community structure, giving rise to 0 versatility. These points are reversed in both scenarios as the sequences are basically reversed. Between these points, the versatility grows because the community structure is less defined as realizations reach different configurations.

III. STRUCTURES IN GROWING BIPARTITE NETWORKS
Additionally, we generate bipartite networks using an evolutionary model of mutualistic webs, through speciation and divergence of weights [13] and then perform numerical simulations to detect the community evolutions (Fig. 6).
In this model, nodes are considered to be either present or absent, with no role to be played by population size. Some properties such as the heterogeneity in degree distribution or the nestedness of ecological mutualistic networks are captured at the same time by the model.
The community structure is stationary during network growth, with the nodes in each community most probably remaining in the same community (Fig. 7). The modularity and the nestedness remain low and constant during the network evolution (Fig. 8), in contrast with their response when the nodes are removed in any of the scenarios presented above.

IV. CONCLUSIONS
We have studied the evolution of the community structure of an empirical ecological bipartite network in the context of extinction of consumer nodes (in this case pollinators). To do so we have introduced 5 different extinction scenarios to account for different extinction dynamics: 1) random node extinction, 2) directed extinctions, 3) generalist extinctions, 4) specialist extinctions and 5) random interaction extinctions. First, we qualitatively observe that during the initial steps the community structure is not affected too much under random node extinctions and specialist extinctions, in contrast to what happens under generalist, directed and random link extinction mechanisms. We next quantify the changes in the organization of the network and its community structure under the different extinction scenarios with a battery of measures. We show that the community structure is reorganized as a function of the fraction of extinction events, signaled for example by the high versatility values at certain moments of the dynamics. Besides that, and most importantly, our result for the effective modularity shows potential for evaluation of risks of ecosystems, if we know under which kind of extinction dynamics the network is suffering. Considering that modularity is a desirable characteristic in ecological ecosystems, as it can buffer the spread of perturbations [8,9], we can conclude that random node and specialist extinctions are always detrimental for the system, and the response of the system is approximately proportional (in terms of loss of modularity) to the percentage of loss of species. For random interaction extinctions modularity remains basically unchanged until 90% of interactions are removed, where the network collapses and modularity suddenly decays. This, therefore, is not such a detrimental extinction dynamics. For directed extinctions, we have also a mostly constant modularity until a drop towards 0 appears when the network fragments. This happens for a fraction of around 0.6 species gone extinct. The problem here is that the response is strongly non-linear and we can go from a relatively high value of the effective modularity to a very low one with just a few species going extinct. This is the same in the case of generalist extinctions, but even more dangerous, as the drop comes at around 35% of species going extinct. Nevertheless, below that critical value for generalist extinctions the effective modularity actually increases, which may be an effective, but dangerous, way of endowing the network with a with higher modularityjust removing a few of the most generalist nodes, without taking too many, as we may totally dismantle the network and end up with low modularity. We postulate, based on our results, that the disappearance of a few generalist species in mutualistic ecosystems might be a way of protecting the system against the spread of perturbations, but that this is a dangerous game, given that if too many are gone extinct, the modularity suddenly drops. Ecosystems are typically robust to the removal of a small fraction of the species, extinction of only a single species positioned at the core of the community cause significant to total network collapse [33] For a growing bipartite network model, in contrast, the community structure is more stationary than in the extinction scenarios of the real mutualistic networks in the sense that the modularity slightly changes as the network grows in comparison to the variation of the modularity with extinctions.
It rests a challenge to widen the results to the weighted case. The results would depend on whether the weights are distributed homogeneously or heterogeneously. Typically, the weight in mutualistic networks captures the fraction of visits of a pollinator to a plant normalized with the total number of visits. As such, the total weight per plant always adds up to one. If the distribution of weights is homogeneous, our expectation is that the communities would behave similarly. A different case would be for weights distributed mostly along a subset of pollinators. In this case, the results can be greatly affected.