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.

Ecological networks may be viewed as a set of species (nodes) connected by interspecific interactions (competition, predation, parasitism and mutualism), represented by the links. Even though interaction strengths are largely unknown, the architecture of the ecological interaction networks has been thoroughly investigated, showing its important role in shaping and regulating community dynamics and in structuring diversity patterns12345678. Several studies recognized the strong impact of the non-random structures of empirical interaction networks on both the resilience (time to return to the steady state after a small perturbation) and the persistence (number of coexisting species at equilibrium) of ecological communities91011121314, and much theoretical effort has been made to understand the relationship between stability and complexity in ecological communities, one of the most debated issues in ecology15161718. In mutualistic networks, where species beneficially interact with each other, a core–periphery structure has been observed ubiquitously19. The network core refers to a central and densely connected set of nodes, while the periphery denotes a sparsely connected non-central set of nodes, which are linked to the core. It has been posited that the architecture of mutualistic networks minimizes competition and increases biodiversity7, community stability (resilience) and persistence20, but other studies have demonstrated that structured mutualistic ecological networks may be less stable than their random counterparts1421. It has also been shown that community stability decreases as community size increases, and that this result holds even for more realistic ecological interactions with a mixing of interaction types (‘hybrid communities')22. Most of the aforementioned studies focused either on the resilience of the system—measured by the maximum real part of the eigenvalues of the community matrix141521—or on the number of species that persist when starting from non-stationary conditions78. However, both approaches have important limitations. Indeed, the maximum real part of the community matrix eigenvalues only describes the rate of recovery from perturbations in the long time limit, providing no information on the transient response. Perturbations can grow significantly before decaying, possibly impacting species' fate (Fig. 1a). A system at its stable stationary state that experiences such initial amplifications of the perturbations is called reactive2324. On the other hand, persistence (measured as the fraction of initial species with positive stationary population density16) is strongly sensitive to initial conditions, the system's distance from stationarity and the choice of model and parameters82526. To garner a better understanding of the effect of perturbations on ecological communities, one should also study how the components of the leading eigenvectors (that is, the right and left eigenvectors associated with the eigenvalue having the largest real part) are distributed, that is, study the localization of the system. In condensed matter physics, localization, also known as Anderson localization27, is the absence of diffusion of waves in a disordered medium, and it describes the ability of waves to propagate through the system. Other approaches (for example, Markov chain models28, or the inverse community matrix29) can be used to study how disturbances propagate in species interaction networks and what their effects are on other species (that is, how many other species do they affect and what is the magnitude of this effect). However, it has been shown that small variations in the interaction strengths may lead to very different model predictions3031. Our theoretical framework may be considered as a complementary methodology to gain information on the general relation patterns between the interaction network architecture and the ability of perturbations to propagate within the system. Our goal in this work is to determine the degree of localization of eigenvectors in mutualistic ecological networks as a function of the network size, structure and interaction strengths, and to study the impact of localization on the perturbation amplitude and spreading within the system. Here we show that localization may be a useful mechanism that impacts on the stability of ecological networks. In fact, localization attenuates (asymptotically) and reduces perturbation propagation through the network. We find that mutualistic ecological networks are indeed localized and localization patterns are correlated with some network topological properties; in particular, heterogeneity in the weighted species degrees promotes localization in the network. Furthermore, the observed localization increases with the size of the ecological communities, highlighting a trade-off between the asymptotic resilience of the system and the attenuation of perturbations.

Figure 1

Propagation of the perturbation through the network.

(a) Trajectory of a perturbation through time. Reactivity (λ) measures whether perturbations grow before decaying; asymptotic resilience λ1 indicates whether perturbations eventually decay; and the asymptotic perturbation amplitude describes the intensity of the perturbation for large time. The principal right eigenvector determines which species will be affected most by the perturbation after its propagation, while the left principal eigenvector controls which species are the most sensitive to the initial perturbation. The weighted degree heterogeneity affects the localization pattern in the network: (b) is a regular graph where each node is connected to six other nodes, while (c) is a power-law scale-free graph2 of the same size and with similar connectance. In both cases, edge weights are randomly extracted from a Gamma distribution. The size and the colour of the nodes indicate the absolute values of the corresponding component of the leading right eigenvector. In b, all species are equally perturbed. In contrast, in c, only few species are affected.

Results

Theoretical framework

The mutualistic interactions of an ecological community can be encoded in a bipartite binary graph represented by its adjacency matrix B containing S nodes (species) that are partitioned into two disjoint sets, one containing the animals (insect pollinators), the other the plants. Each of the L (undirected) edges connects two nodes, one in the set of animals (of size A) and the other in the set of plants (of size P), that is, B=1 if insect k and plant l interact. S=A+P is the total number of species in the community. We analyse 59 bipartite binary networks available from the interaction web database9, and we construct the S × S community matrix Φ describing the linearized system dynamics, by assigning to each animal–plant interaction a positive ‘weight' (see Methods and Supplementary Fig. 1). Let x(t) be the S-component vector describing the abundance of the S populations at time t. The propagation of a given small perturbation =(ξ1,ξ2,...,ξ) acting on the system at stationarity will lead to small departures, δx(t), from the stationary state x* and can be studied by the linearized system of coupled differential equations , where δx(t)=x(t)−x* (with δx(0)=), which in turn can be studied in terms of the eigenvectors and eigenvalues of Φ, known as community matrix (see Methods). In particular, the asymptotic behaviour of the perturbed systems can be analysed in terms of the largest eigenvalue λ1 and corresponding left and right eigenvectors u1 and v1, that is, for large t, where is the amplitude associated with the asymptotic propagation of the perturbation through the ecological network (Fig. 1 and Methods).

Clearly, the relaxation of a system to its equilibrium state after a perturbation is not uniquely controlled by the leading eigenvalue of the community matrix. All the eigenvalues contain information on the timescales involved in the relaxation, while the corresponding eigenvectors determine how the perturbation spreads and relaxes in different species. The leading right eigenvector, in particular, sets the relative vulnerability of species and how they are affected by perturbations in the long run. If this eigenvector is localized, that is, if only few components/species have non-negligible values, then a perturbation after its propagation involves only few species. On the other hand, the left leading eigenvector indicates which species are most hit by the perturbation before its propagation. It also plays an important role in modulating the amplitude of the perturbation (that is proportional to ·u1—see Methods). Consider for example a 4 × 4 community matrix for which, and . If =(0,0,1,1), then u1·ξ=0 and the perturbation decay time will be very fast (controlled by λ2, rather than λ1—see equation (1) in the Methods section). On the other hand, if =(1,1,0,0), then the system asymptotic recovery time will be longer (proportional to 1/λ1), but only species 2 and 3 will be affected at these timescales. As a general trend, we will show that localization mainly depends on the heterogeneity of the network weighted degrees (or strengths s=(s1,s2,...,s)): in the case of high variability in these strengths, the system display localization (Fig. 1b,c). The behaviour immediately after the perturbation (that is, in the limit t→0+) can be analysed by studying its reactivity23, defined as the maximum amplification rate over all initial perturbations, and immediately after the perturbation. It can be shown23 that the reactivity λ can be computed as the maximum eigenvalue of H=(Φ+Φ)/2, the symmetric part of Φ. If λ1<0 and λ>0 then the equilibrium point is stable but reactive. Because λ≥λ1 (ref. 23), the reactivity can also be used to develop an early warning signal for systems approaching a non stable stationary state24. If the eigenvector w corresponding to λ, is also localized, then it means that the perturbation magnitude on these localized species will tend to grow (Fig. 1a), that is, in the short time, these species will be the most affected by the perturbation.

Localization patterns

We compare localization patterns of 59 empirical mutualistic networks and two corresponding random null models. In the first null model, we randomize the interactions while keeping the networks connected. In the second null model, we randomize the interactions, but we also constrain the network degrees sequence {k1,k2,...,k} to be as in the corresponding empirical networks (see Methods). To measure localization, we use the inverse participation ratio (IPR)27, the classical way to quantify how many relevant components are observed in the leading eigenvectors (see Methods). The degree of localization increases as IPR increases. If IPR is one, then only one component of the eigenvector is non-zero. We quantify the presence of localization by computing the rIPR defined as the ratio between the IPR of each real empirical network and the IPR of the corresponding random null model.

As Fig. 2a,b shows, most of the empirical networks are significantly more localized in both the right and left leading eigenvectors with respect to null model 1, while they have the same level of localization of null model 2 (Fig. 2d,e,i and Table 1). These results suggest that it is the core–periphery network structure of empirical systems (a manifestation of heterogeneous degree distributions) that is responsible for their higher localization: once we constrain the degree distributions to be fixed (that in the case of mutualistic networks are most likely approximate truncated power laws2), then null model 2 generates localization patterns very similar to those observed in empirical mutualistic networks (Table 1). Nodes strength s (or weighted degrees) also play a crucial role. In fact, an adjacency network with core–periphery structure, but having ‘anti-nested' distributed weights13, will not be localized because, contrary to its degree distributions, the weighted degree distribution will be homogeneous (see also Supplementary Methods, section 5 and Supplementary Fig. 15). The localization of w displays the same patterns (Fig. 2c,f), and we found that species that are the most affected by the perturbation at short timescales are also those that absorb most of the perturbation asymptotically—indicating a limiting capability of the the perturbation to propagate through the network. In fact, the position of the localized components for w is most likely to be the same of those for v1 and u1 (see Fig. 3 and Supplementary Methods, section 7).

Figure 2

Localization patterns and effect on the asymptotic amplitude.

(a–c) Localization of the leading eigenvectors for null model 1 versus empirical mutualistic networks (where , with q1=v1, u1 or w, with right eigenvector v1 in red, left eigenvector u1 in blue and reactive eigenvector w in green). Points represent average value of 1,000 randomizations, bars indicate the standard deviation. 1:1 line represents the value of the empirical mutualistic networks. (d–f) Same for null model 2. (g) Number of (localized) components as a function of the community size: a significant correlation is observed (Spearman Rank Test=0.715). (h) Effect of the localization on the asymptotic amplitude for the simulated perturbation ξ on empirical mutualistic networks with respect to null model 1 . In ≈85% of the cases, the perturbation amplitude is significantly attenuated (P value≤0.05). (i) P values of the observed values of localization and asymptotic amplitudes in empirical networks with respect to null model 1. Parameters here are δ=0.5 and γ0=1.

Table 1

Statistics of localization patterns.

	δ=0.5	δ=0	δ=−0.5	ΝΜ ♯
IPR/〈IPRran〉>1
u1	≈86%	≈76%	≈39%	1
v1	≈61%	≈64%	≈42%	1
wH	≈83%	≈72%	≈50%	1
u1	≈20%	≈20%	≈19%	2
v1	≈27%	≈27%	≈11%	2
wH	≈24%	≈20%	≈19%	2
A1/A1ran<1
	≈81%	≈81%	≈70%	1
	≈17%	≈15%	≈14%	2

Statistics of localization patterns for three different ecological scenarios (described by δ) given by the fraction of localized empirical networks with respect to generated null models (NM) for different parametrization (rIPR>1 and P value<0.05) and corresponding asymptotic attenuation of the perturbation ( and P value<0.05). v1, u1 and wH are the right, the left and the reactive eigenvectors, respectively. We found that most of the empirical networks are indeed localized with respect to NM1, but not with respect to NM2.

Figure 3

Localization in Norfolk mutualistic insect-grasslands community.

Example of localization in a real mutualistic community35 with S=78 species. The 61 insects and the 17 plants species are sorted according to degree. The flower species Leucanthemum vulgare is the most localized in each of the three eigenvectors, v1, u1 and w and corresponds to the species with the highest species degree and strength. The parameters are δ=0.5 and γ0=1, yielding a stable but reactive equilibrium (λ1<0, λ>0).

Relation between localization and network properties

Localization patterns in empirical mutualistic communities depend on both the size and the connectance of the species interaction network (Figs 2g and 4). While the leading eigenvalue λ1—the one with the largest real part which controls the relaxation time—increases for increasing community size15 (assuming that γ0 does not scale22 with S), we observe (Fig. 4a) an interesting strong positive correlation between community diversity (network size) and localization (rIPR). We also note that in empirical pollination communities, networks size is negatively correlated with network connectance14. Network connectance, in turn, is negatively correlated with localization (Fig. 4b): the higher the connectance, the higher is the ability of perturbations to propagate through the network (a general property observed also in financial networks32, and socio-environmental interdependent systems3334), and thus the lower the level of localization. The strong positive correlation between network size and localization leads to a trade-off between localization and asymptotic resilience (in terms of λ1—see Table 2). This result may shed light on the celebrated complexity-diversity paradox1715: the less an ecological community is resilient, the more it is localized and asymptotic perturbation is attenuated (Fig. 5). We finally note that localization is positively correlated with the variance of the weighted degree distribution, see Fig. 4c. This correlation reflects the fact the localization is a manifestation of the heterogeneity of the network topology. In the Supplementary Tables 1–3, we also calculate the correlations between localization (rIPR) and the topological properties of the networks under different parametrizations and with respect to null model 2. We found that in this latter case the correlations are not significant. We can thus conclude that it is indeed the heterogeneity in the weighted degree distribution, which is the key structural aspect of ecological networks that is related to localization.

Figure 4

Relation between localization and network topological properties.

Relation between localization (rIPR is the ratio between the IPR of the empirical network and the IPR of the corresponding random null model 1) and (a) networks size (S); (b) network connectance measured as the fraction of observed and possible links (C=L/(S(S−1)) with L=number of links); (c) network heterogeneity (variance of the weighted degree distribution). The parametrization used here is δ=0.5 and γ0=1. Results are reported for both right (in red) and left (in blue) leading eigenvectors (v1 and u1).

Table 2

Correlations between network topological and spectral properties.

ρ(x,y)	x	y	P value
≈0.760	S	rIPR[v1]	<10−4
≈0.754	S	rIPR[u1]	<10−4
≈0.800	S	rIPR[wH]	<10−4
≈−0.662	C	rIPR[v1]	<10−4
≈−0.754	C	rIPR[u1]	<10−4
≈−0.749	C	rIPR[wH]	<10−4
≈−0.769	S		<10−4
≈0.806	C		<10−4
≈−0.477	λ1		<10−4
≈0.578	λ1	rIPR[v1]	<10−4
≈0.390	λ1	rIPR[u1]	<10−4
≈0.468	λ1	rIPR[wH]	<10−4
≈0.460	σs2	rIPR[u1]	<10−4
≈0.535	σs2	rIPR[v1]	<10−4
≈0.45	Degree [k]	v1	<10−4
≈0.46	Strength [s]	v1	<10−4
≈0.67	Degree [k]	u1	<10−4
≈0.61	Strength [s]	u1	<10−4
≈076	Page rank	v1	<10−4
≈0.94	Eig. centrality	v1	<10−4

Correlations ρ(x,y) measured using Spearman Rank Test (parametrization δ=0.5) using Holling Type I model with a=40 and b=0.05). rIPR refers to null model 1. v1, u1 and wH are the right, the left and the reactive eigenvectors, respectively. k and s denote the network degree and strength vectors.

Figure 5

Relation between size and relative asymptotic amplitude.

(a) Relationship between the relative asymptotic amplitude (with respect to null model 1) and network size: attenuation increases for increasing community size (that is, decreases with S). (b) Asymptotic amplitude for perturbed empirical communities is compatible with that one generated by null model . Indeed null model 2 generates networks with the same level of localization of the empirical pollinator networks (rIPR≈1), and thus attenuation is not observed. Red points correspond to the parametrization δ=0.5, while blue squares represent the mean-field case (δ=0).

Attenuation of perturbation propagation and amplitude

To make analytical progress and to better understand the impact of the localization on the amplitude of the perturbation spreading throughout the network, we analyze the mean-field case (δ=0), and we assume that the perturbation vector is a constant, that is, =ξ01 (a S vector of unit components) and v1≈u1. Under this assumption, we are able to prove that and thus is maximum when which corresponds to a state of minimum localization of Φ, while is minimum when v1,=δ for some i, which corresponds to the fully localized case. Thus, the mean-field approximation with constant perturbation suggests that localization in the system reduces the amplitude of the principal mode of the perturbation wave. Indeed, our numerical simulations (Figs 2h and 5a and Supplementary Methods, section 4.1) confirm that a localized structure leads to a decrease in the principal amplitude of the perturbation also beyond the mean-field case (that is, δ≠0). Moreover, following the localization trend, the perturbation damping increases with the size of the system: the larger the ecological network, the stronger is the attenuation due to the system localization (Fig. 2d). Also the reverse is true: if a network is not significantly more localized than its corresponding null model, then no attenuation is observed (Fig. 5b, Table 1 and Supplementary Methods, section 4.2).

Discussion

We have developed a comprehensive theoretical framework to evaluate the relationship between the species interaction network architecture and the impact of a given perturbation on ecological mutualistic networks. Localization has thus two beneficial effects on ecological network robustness: (a) only a very low proportion of species in the community are significantly affected by a perturbation spreading throughout the network, and (b) localization leads to an attenuation of the perturbation effects on the system. These results are robust with respect to variation of the parameters (see Supplementary Methods, sections 1, 2 and 4, and Supplementary Figs 2–14) and thus hold for very general parametrization of the interaction strengths (that in general are unknown—see Supplementary Methods, section 1.1). We thus have shown that the eigenvectors of the community matrix play a crucial role in determining the impact and the propagation of the perturbation through the system. We found that the positions of the localized components of the principal eigenvectors strongly correlates with nodes degree centrality, species strength s, eigenvector centrality and page-rank centrality (Table 2, Supplementary Methods, section 7 and Supplementary Fig. 17): the proposed framework thus allows one to identify those species which are affected the most by a given perturbation. Interestingly, these are species with many mutualistic interactions, and on average with higher population abundances with respect to specialist species14. For example, in Fig. 3 we show the eigenvectors components of the leading eigenvectors for the community and reactivity matrix associated with the insect-grasslands ecological community in Norfolk35. For each of the 61 insects and the 17 plants, we can calculate the corresponding values of v1, u1 and w. Species within each class (pollinators and flowers) are then sorted according to their degree (number of interacting partners they have). We note that the flower species Leucanthemum vulgare, the species with the highest species degree and a high density, is the most localized species in the community for both v1, u1 and w, and it is the one that is likely to absorb most of a potential perturbation affecting the whole community.

A general emerging pattern observed for the mutualistic communities analysed in this work is that, while these systems are less resilient for increasing biodiversity (May's result15), localization and the corresponding perturbation attenuation increase with increasing species diversity. In other words, these mutualistic systems experience a trade-off between resilience and localization: small communities are faster in recovering their stable state after a perturbation, but they are less localized and the perturbation will have an impact on most of the species. On the other hand, large communities are less resilient (that is, need larger time to return to a stable state), but only few species will be affected by the perturbation, and amplitude of the perturbation will be attenuated while spreading through the network. The proposed theoretical analysis can be easily applied to networks with other interaction types, and it is illustrative of the potential of this new metric.

Methods

Data parametrization

To be able to describe different ecological scenarios, and following recent models in the literature2636, we parameterize the weighted interaction matrix as for i≠j and W=−d (see Supplementary Methods, section 1), where B is the adjacency matrix of the species interaction matrix, indicating presence (B=1) or absence (B=0) of interactions among species, k=∑B is the number of mutualistic partners of species i (species degree), s=∑W is the species strength (or weighted degree) and γ0 is a parameter describing the basal mutualistic strength, while δ a trade-off parameter controlling the relation between mutualistic interaction strength and species degree. Following a Holling Type I population dynamics model, we then build the community matrix Φ as , where denotes the stationary population abundance of species i, and we model it as random variable drawn from a Gamma distribution so to have an average species population abundance 〈x*〉=1 and of standard deviation (Supplementary Methods, section 1). By varying the parameter δ we investigate: (a) The architecture with constant interaction strength (‘mean-field' case714, δ=0); (b) The architecture with interaction strength-degree trade-off (δ>0), for example, specialist species interact stronger than the generalists one)3738. (c) Architecture where generalist species interact stronger than specialist species (δ<0). Using this parametrization, for a fixed δ the stability and reactivity of ecological communities can be controlled by the value of the basal mutualistic strength γ0, and the intra-specific competition d.

Null models

We generate two different random null models (NM) for Φ (Φ−ran). (NM1) We assign the L links in the adjacency matrix at random while keeping the network connected, and then parametrize it in the same way we do for empirical networks. (NM2) We assign the L links at random, but constraining the degree sequence (k1,k2,...,k) to be the same of the corresponding adjacency matrix B and then parametrize it in the same way we do for empirical networks. Our results are compared with 1,000 realizations of each of the null models. For all other details we refer to the Supplementary Methods, section 3.

Perturbation analysis

The effect of a given perturbation =(ξ1,ξ2,...,ξ) acting on the system at time t=0 will propagate in time obeying with initial condition x(0)=. The solution of the latter equation can be written in terms of the eigenvectors and eigenvalues of Φ:

where u, v and λ( are respectively the left, the right eigenvectors and the corresponding eigenvalues of the linearization matrix Φ. We ordered the eigenvalues so that 0>λ1>Re[λ(2)]>...>Re[λ(] (we note that in our case, as φ≥0, the Perron–Frobenius theorem holds and λ1=Re[λ(1)]). For simplicity, we will denote by the amplitude associated with the α-th mode of the perturbation.

Localization and effect on stability

We measure the localization using the inverse participation ratio IPR27, that is, , where q1=v1, u1 or w. In particular, we identify localization patterns by computing the rIPR, that is, the ratio between the IPR of each real empirical network and the IPR of the corresponding random null model: . The average 〈·〉 is taken among different realizations of Φ−ran. If rIPR is significantly larger than one, then the system is localized. Otherwise we say that the system is not localized. We can also quantify the number of localized species by setting a threshold θ and count the fraction of species with a leading eigenvector component larger than that threshold, that is, v1(i)>θ or u1(i)>θ. We set would correspond to the extended, non localized case). We quantify how the architecture of the ecological networks affects the impact of a simulated perturbation on the system by comparing the outputs , λ1, λ, v, u, w and ρ with respect to the corresponding random null models. We consider different type of perturbations. In Supplementary Methods (section 6 and Supplementary Fig. 16), we present results for (a) a noise which is independent of species characteristics, that is, drawn from a normal distribution of mean 1 and variance ζ2; (b) a noise that is species dependent, that is, proportional to the degree of each species (ξ(i)∝kξ(i)). In the main text we show results for a perturbation combining both types of noise, that is, all=+.

Additional information

A link to the Mathematica notebook with main functions needed to compute localizations and effect on stability of mutualistic ecological communities is here provided: https://github.com/suweis/Effect-of-Localization-on-the-Stability-of-Mutualistic-Ecological-Networks.

Additional information

How to cite this article: Suweis, S. et al. Effect of localization on the stability of mutualistic ecological networks. Nat. Commun. 6:10179 doi: 10.1038/ncomms10179 (2015).