Effect of localization on the stability of mutualistic ecological networks

The relationships between the core–periphery architecture of the species interaction network and the mechanisms ensuring the stability in mutualistic ecological communities are still unclear. In particular, most studies have focused their attention on asymptotic resilience or persistence, neglecting how perturbations propagate through the system. Here we develop a theoretical framework to evaluate the relationship between the architecture of the interaction networks and the impact of perturbations by studying localization, a measure describing the ability of the perturbation to propagate through the network. We show that mutualistic ecological communities are localized, and localization reduces perturbation propagation and attenuates its impact on species abundance. Localization depends on the topology of the interaction networks, and it positively correlates with the variance of the weighted degree distribution, a signature of the network topological heterogeneity. Our results provide a different perspective on the interplay between the architecture of interaction networks in mutualistic communities and their stability.

The parametrization of the interaction matrix w ij = γ 0 B ij /k δ i allow us to consider different scenarios: A) Mean field case (δ = 0) where there is no trade-off between species mutualistic strength and species degree; B) AInteraction strength-degree trade-off (δ = 0.5), where a weak w ij term is associated with a strong w ji term (e.g. a strong dependency of the plant on the pollinator and a weak dependency of the pollinator on the plant); C) Positive strength-degree relation (δ = −0.5), where we set a positive correlation between species interaction strength and degree.
0.06 0.08 0.10 0.12 0.14 0.16 0.18             Figure 4: Localization in mutualistic networks: mean field case. Comparison between IP R (v 1 in red, u 1 in blue) for each of the 59 empirical pollinator networks (horizontal solid lines), and Box-Whisker plots for the IP R of the 1000 networks generated through the null model 1                    Figure 10: Localization in mutualistic networks: mean field case. Comparison between IP R (v 1 in red, u 1 in blue) for each of the 59 empirical pollinator network (horizontal solid lines), and Box-Whisker plots for the IP R of the 1000 networks generated through the null model 2 (δ = 0). 9            Table 1: Correlations between network topological and spectral properties under NM1. Correlations ρ(x, y) measured using Spearman Rank Test (for parametrization δ = 0 using Holling Type I model with a=40 and b=0.05). rIP R refers to null model 1   (1), and compute the S × S community matrix Φ describing the ecological community dynamics around one of its stationary stable points. In this case Φ will depend on both the interaction weights (or strenghts) w ij , and the species biomass at equilibrium x * . As interaction weights are unknown from the data, we parametrize the matrix W describing the interaction strengths so to explore different kinds of ecological scenarios. Specifically we set (assuming no inter-species competition): where k i is the degree (number of mutualistic partners) of species i, γ 0 is the basal mutualistic strength, δ K is the Kronecker delta, while d i describes the species i self-limitation effect (for simplicity here we set d i =1 ∀i = 1, ..., S and γ 0 = 1.). δ ∈ [0, 1] sets the relationship between mutualistic interaction strength and species degree, and we have investigated different parametrizations corresponding to different ecological scenarios (see Supplementary Figure 1): A) Mean field case (δ = 0) where there is no trade-off between species mutualistic strength and species degree (2; 1)); B) Interaction strength-degree trade-off (δ = 0.5), where a weak w ij term is associated with a strong w ji term (e.g. a strong dependency of the plant on the pollinator and a weak dependency of the pollinator on the plant) (3; 4); C) ticular, we identify localization patterns by computing the rIP R, i.e. the ratio between the IP R of each real empirical network and the IP R of the corresponding random null model: rIP R i = IP R i / IP R ran i . The average · is taken among 1000 realizations of Φ ran . If rIP R is significantly larger than one, then the system is localized. Otherwise we say that the system is not localized. Finally the amplitude associated with the asymptotic propagation of the perturbation through the ecological network is defined as (1) in the main text and methods section). In all cases we set d = 1, and we chose the parameters of the Gamma distribution so to have an average species abundance of x = 1 and variance σ 2 x < σ 2 φ (in particular we set a = 40, b = 0.025 and γ 0 = 1). We thus obtain 59 community matrices Φ (one for each empirical pollinator network) for each of the three different ecological scenarios δ = 0.5, δ = 0. and δ = −0.5.

Null Model 1
As shown in Figures S2, localization patterns for δ = −0.5 are qualitatively the same of those presented in the main text (Figure 2a-c) for δ = 0 and δ = 0.5. Figures S3-S5 show the statistics of the localization patterns for each of the 59 networks and for different parametrizations through the Box-Whisker plots. The ends of the whiskers represent the minimum and maximum, whereas the ends of the box are the first and third quartiles and the white bar denotes the median. Single points represent outliers realizations. We finally study the effect of localization of the attenuation of the asymptotic amplitude A 1 . For each of the 59 networks, and for the different parametrization we calculate A 1 and compare it to the distribution of A ran 1 for 1000 random networks generated by null model 1. As Figures S6-8 show, localization indeed attenuates the perturbation amplitude (A 1 < A ran 1 ).

Null Model 2
As already shown in the main text, null model 2 generate most of the times random networks that have similar localization patterns to those observed in empirical pollinator networks. These results suggest that it is the heterogeneous degree distribution of empirical systems (16) the responsible of their higher localization: once we constraint the degree distributions to be fixed, then null model 2 generates very similar localization patterns of those observed in empirical mutualistic networks. Nodes weighted degrees (or strength s i -see below) also play a crucial role. In fact, a network with binary core-pheriperhy structure, but having anti-nested distributed weights (8), will not be localized because, contrarily to its degree distribution, the weighted degree distribution will be homogeneous (see section below). Also in this case we study the effect of localization of the attenuation of the asymptotic amplitude A 1 for each of the 59 networks, and for the different parametrization we calculate A 1 and compare it to the distribution of A ran 1 for 1000 random networks generated by null model 2. As Figures S12-14 show, as a consequence of the fact that the null mode 2 generate similar localization patterns to those observed in real data, then we find no significant attenuation of the asymptotic amplitude (A 1 ≈ A ran 1 ). Figures S9-S11 show the statistics of the localization patterns for each of the 59 networks and for different parametrizations through the Box-Whisker plots with respect null model 2. Only very few empirical networks are localized with respect this null model.
Finally we test if the correlation between the rIP R calculated through null model 1 and the network topological properties shown in the main text ( Figure 4) and Table 2, still hold when we calculate the rIP R using null model 2. As Supplementary Table 3 shows, we found that the correlations between rIPR (calculated through null model 2), size S and connectance C are now not significant, and thus we can conclude that it is indeed the heterogeneity in the weighted degree distribution the key structural aspect of ecological networks that is related to localization.

Crucial role of the strength distribution.
In this section we show that, when considering weighted networks, what determines localization of the leading eigenvector(s) is not necessarily the degree distribution, but rather the distribution of nodes strengths (s i = j w ij ). We consider two different networks: a random networkk (Erdös-Renyi (17)) and a scale-free (Barabasi-Albert (18)) network, with same size S and number of links L. The second network is much more nested (as a consequence of its heterogeneous degree distribution) than the first one. However, we arranged the weights in such a way that the random network has a larger variance in the strengths with respect to the heterogeneous network. We call anti-nested(8) this arrangement of weights (the binary matrix is nested, while the weighted matrix is less nested than its random counterpart). In this case the scale-free nested network displays a lower localization than the random network (i.e. a network with random connection and weights drawn from the same distribution). An example of this result is shown in Supplementary  Figure 15. 6 Results for different type of perturbations.
We here present here the results of the effect of localization on stability (in terms of the asymptotic amplitude of the perturbation A 1 ) for two type of perturbations: a) A noise ξ D that is independent of species characteristics.
In particular, we model ξ D as drawn from a normal distribution N (1, ζ). b) A noise ξ E that is species dependent, i.e. proportional to each of the species degree (ξ E (i) ∝ k i ξ D (i)). For simplicity, in this and in the following sections we used the model independent parametrization (see section 2). Results are robust with respect to different parametrizations of the community matrix.
7 Which species are more localized?
In this section, we study the topological properties of the nodes (species) which have the largest components of the eigenvectors. To accomplish this, we quantify the number of localized species by setting a threshold θ (node i is localized if v 1 (i)(u 1 (i)) > θ) and look at the centrality properties of these nodes. The centrality of a node, is a measure that assigns a given "importance" to that node in the network (17). There are several ways to denote this importance. The simplest one is the degree centrality k i , the number of connections that a node i has. But one may want to consider also "the quality" of the neighbors of a given node, not just its number. Then the eigenvector centrality has been introduced (17), where the centrality of node i is given by the i-th component of the principal eigenvector (i.e. relative to the maximum eigenvector of the graph matrix). We set θ = 1.5/ √ S (v 1 (i) ∼ u 1 (i) = 1/S would correspond to the extended, non localized case). We also compute the PageRank (PRC) measure, which is used by the Google web search corporation as a central part of their web ranking technology (17). We note that all the centrality measures are calculated on the symmetric (binary) adjacency matrix, and thus do not depend on the parametrization of the ecological networks. In all cases we find that the localized nodes are on average the most central nodes of the networks, i.e. nodes belonging to the core of the core-periphery (nested) structure of the mutualistic species interaction networks. These results suggest that the peripheral nodes, i.e. the specialized species, are more "protected" by perturbations affecting the system dynamics. The full statistic of the correlation among these quantities for the parametrization δ = 0. can be found in Table 1 (results for δ = −0.5 are qualitatively the same -see Table 2).
Both pair of vectors {w H , u 1 } and {u 1 , v 1 } are strongly correlated (see Supplementary Figure 17), indicating the ability of the network to contain the spread of the perturbation within the network, i.e., those species that are more affected by the perturbation at time t → 0 + (components i with large w H (i)), are those (or very close) to those that are most hit by the perturbation asymptotically (components i with large u 1 (i)), that in turn are those species that are most affected asymptotically by the perturbation after its spread through the network (components i with large v 1 (i)).