## Abstract

A classic measure of ecological stability describes the tendency of a community to return to equilibrium after small perturbations. While many advances show how the network architecture of these communities severely constrains such tendencies, one of the most fundamental properties of network structure, i.e. degree heterogeneity—the variability of the number of links associated with each species, deserves further study. Here we show that the effects of degree heterogeneity on stability vary with different types of interspecific interactions. Degree heterogeneity consistently destabilizes ecological networks with both competitive and mutualistic interactions, while its effects on networks of predator–prey interactions such as food webs depend on prey contiguity, i.e. the extent to which the species consume an unbroken sequence of prey in community niche space. Increasing degree heterogeneity tends to stabilize food webs except those with the highest prey contiguity. These findings help explain why food webs are highly but not completely interval and, more broadly, deepen our understanding of the stability of complex ecological networks.

## 1. Introduction

Understanding the intricate relationship between the structure and dynamics of complex ecological systems has been one of the key issues in ecology [1–4]. Equilibrium stability of ecological systems, a measure that considers an ecological system stable if it returns to its equilibrium after a small perturbation, has been a central research topic for over four decades [1,5–15]. Empirical observations suggested early on that communities with more species are more stable, i.e. a positive diversity–stability relationship [16]. Yet, these intuitive ideas were challenged by the pioneering work of May [1,2], who rigorously proved a negative diversity–stability relationship based on linear stability analysis of randomly constructed ecological communities. The contradiction between these findings forms the eminent ‘diversity–stability debate’ [6,7].

May's seminal work considered community matrices *M* of size *S* × *S*, where *S* is the number of species, the off-diagonal element captures the impact that species *j* has on species *i* around a feasible equilibrium point ** x*** of an unspecified dynamical system describing the time-dependent abundance vector

**(**

*x**t*) of the

*S*species [1,2]. Since empirical parametrization of the exact functional form of

**(**

*f***(**

*x**t*)) is extremely difficult for complex ecological systems, May assumed that

*M*'s are randomly drawn from a distribution with mean zero and variance

_{ij}*σ*

^{2}with probability

*C*and are 0 otherwise. Hence

*σ*represents the characteristic interaction strength, and

*C*is the ratio between actual and potential interactions in the ecological system. For simplicity, all diagonal elements were chosen to be −1, which represents the intrinsic dampening time scale of each species so that if disturbed from equilibrium it would return with such a dampening time by itself. May found that for random interactions drawn from a normal distribution , a randomly assembled system is stable (i.e. all the eigenvalues of the community matrix

*M*have negative real parts) if the so-called ‘complexity’ measure . This implies that higher complexity (e.g. due to larger

*C*or

*S*) tends to destabilize community dynamics [1,2].

The reason May's result continues to be influential more than 40 years later is not because complex ecological systems have to be unstable, but because real ecological systems presumably have some specific structures that allow them to be stable despite their complexity [4]. In other words, nature must adopt what May called ‘devious strategies’ to cope with this diversity–stability paradox. One of such strategies is the existence of particular interspecific relationships observed in nature, e.g. predator–prey, competition and mutualism (see electronic supplementary material, §I). Recently, Allesina & Tang refined May's result and provided analytic stability criteria for these three types of interspecific interactions [14]. They found remarkable differences between predator–prey interactions, which are stabilizing, and mutualistic and competitive interactions, which are destabilizing. These newer findings allow many different strategies employed by nature to be more powerfully tested with the refined stability criteria as a reference point.

### 1.1. Degree heterogeneity

The above results rely on the key assumption that the underlying network structure is almost completely random. Indeed, the construction of the community matrix *M* follows almost the same procedure as the classical Erdős–Rényi (ER) random graph model [17]. However, just like many other real-world complex systems, the underlying networks of ecological systems are far from random. Instead, they often display importantly non-trivial topological features, e.g. *degree heterogeneity* (the distribution of species among different levels of generality) [18–21], *nestedness* (the tendency for the specialists' niches to be contained within generalists' niches) [22,23] and *modularity* (the tendency of species to be linked to more species inside groups of species than to species outside those groups) [24], and the correlation between interaction strengths [25]. It has been recently shown that the network architecture favouring stability fundamentally differs between trophic and mutualistic networks [12]. For example, a highly connected and nested architecture promotes community stability in mutualistic networks [11,23], whereas the stability of trophic networks is enhanced in modular and weakly connected architectures.

Despite these and other remarkable results, much about the relationship between network structure and stability of complex ecological systems remains unknown. One fundamental issue concerns the impact of degree heterogeneity, which determines the distribution of species among those with different numbers of predator and prey of species. This distribution tightly constrains the amount and type of niche overlap and competition between species [26]. Typically, the degree heterogeneity of a network can be represented by , where is the mean degree of the network, *k _{i}* is the degree of species

*i*(e.g. species

*i*'s total number of predator or prey species), represents the raw second moment and

*P*(

*k*) is the degree distribution of the network. Note that by definition

*ξ*≥ 1. Also, the higher the degree heterogeneity, the larger the value of

*ξ*. For example, scale-free (SF) networks [27] that have power-law degree distributions have higher degree heterogeneity than ER random networks [17] with Poisson degree distributions. A key advance in understanding complex networks over the last 15 years has been how significantly degree heterogeneity affects many network properties and dynamical processes, from error and attack tolerance [28,29], to graph spectra [30], epidemic spreading [31], interdependent fragility [32] and controllability [33]. It is reasonable to expect that degree heterogeneity would also affect stability of complex ecological systems. Indeed, studying the stability of ecological networks with degree heterogeneity has been highlighted as one of the main challenges in the recent list of open problems [34] including determining level of the controversial property of nestedness [23,35]. To approach this issue in a systematic fashion, we follow May's model-independent framework [1], helping us avoid the difficulty of parametrizing the exact dynamics of complex ecological systems.

## 2. Results

### 2.1. Impact of degree heterogeneity on the stability of ecological networks with random interactions

For an ecological system with bidirectional random interactions, we generate its underlying network topology first, using three different network models (see Material and methods): multi-modal [36], ER [17] and SF [27], with given mean degree *k* and effective connection probability . Then we construct the community matrix *M* as follows: (i) set all the diagonal elements *M _{ii}* = −

*d*with

*d*an arbitrary constant, representing the intrinsic damping time scale of each species; (ii) draw the off-diagonal element

*M*from a normal distribution whenever there is a link from species

_{ij}*j*to species

*i*. The real part of

*M*'s most positive eigenvalue is given by (see electronic supplementary material, §III.A, for the proof). Re(

*λ*

_{m}) has to be negative to ensure the equilibrium stability, yielding the stability criterion 1.1

Thus, any factor that increases (or decreases) Re(*λ*_{m}) will destabilize (or stabilize) the ecological system, respectively. Increasing *ξ* tends to destabilize the ecological system with random interactions. Note that for random networks generated by the ER model, *ξ* → 1 for large 〈*k*〉. Hence for large *S*, equation (1.1) naturally recovers May's classic result [1]. We also consider general directed networks with random interactions (not necessarily bidirectional) and obtain a similar stability criterion (see electronic supplementary material, §III.B).

Figure 1*a* shows the impact of degree heterogeneity on the stability of ecological systems with random interactions. The underlying architectures are constructed from different network models with tuneable degree heterogeneity *ξ* (electronic supplementary material, Sec. II). We find that when the complexity measure is fixed, Re(*λ*_{m}) grows monotonically as *ξ* increases, implying that larger *ξ* tends to destabilize an ecological system with random interactions. Moreover, our numerical results suggest that , which agrees well with our analytical prediction (electronic supplementary material, §III.A). This finding can be further illustrated by the distribution of *M*'s eigenvalues. We find that for increasing *ξ*, the radius of the cycle encompassing all the eigenvalues becomes larger, hence Re(*λ*_{m}) increases, tending to destabilize the ecological system. Interestingly, we also find that as *ξ* increases the eigenvalue distribution becomes more non-uniform with very high density of eigenvalues concentrated around the centre of the circle (figure 1*b–d*).

### 2.2. Impact of degree heterogeneity on the stability of ecological networks with predator–prey interactions

For ecological systems with predator–prey interactions, i.e. whenever *M _{ij}* > 0 then

*M*< 0, we generate the underlying network using five different models (see Material and methods): multi-modal, ER, SF, cascade [37] and niche model [38], with given mean degree 〈

_{ij}*k*〉. When there is a directed edge from species

*j*to

*i*, we draw the off-diagonal element

*M*from the half-normal distribution and

_{ij}*M*from . We still set all the diagonal elements

_{ji}*M*= −

_{ii}*d*. (We also considered the presence of a small fraction of positive diagonal elements in the community matrix

*M*and we found that they do not change our main conclusions, see electronic supplementary material, §I.D for details.) As shown in figure 2

*a,b*, for simple model networks (3-modal, ER and SF), large degree heterogeneity hampers the stability of predator–prey ecological systems, yet moderate heterogeneity is stabilizing. This non-monotonic behaviour can be further illustrated by the eigenvalue distribution (figure 2

*c*). We find that as degree heterogeneity

*ξ*increases, not only does the eigenvalue distribution become non-uniform but also the boundary that encompasses all eigenvalues changes from ellipse to bow-tie shape (figure 2

*c*). This alteration of boundary shape induces the non-monotonic behaviour of Re(

*λ*

_{m}) for varying

*ξ*, which cannot be explained by any existing theories [1,2,14,39].

In empirical food webs, there exist trophic hierarchies and niche structures [40,41]. To reproduce these structural properties, we employ the widely used cascade [37] and niche [38] models. When the number of species *S* and the connection probability *C* are fixed, we tune the degree heterogeneity *ξ* to assess its impact on stability. As shown in figure 3*a*, for networks generated by the niche model degree heterogeneity *ξ* tends to destabilize predator–prey ecological systems, in contrast with networks generated by the cascade model for which degree heterogeneity is stabilizing (figure 3*c*). We also calculate the eigenvalue distributions of the community matrix associated with the niche model and the cascade model, finding that there is no drastic shape alteration of the boundary encompassing all the eigenvalues (figure 3*b,d*).

A primary distinction between the cascade and niche models is *prey contiguity,* which describes the tendency for species to consume a contiguous sequence of prey in tropic niche space of the whole community [38,41] (figure 4*a*). The remarkable difference between figure 3*a,c* prompts us to systematically explore the effect of prey contiguity. We adopt the relaxed niche model [42] to generate underlying networks with tuneable prey contiguity *g* such that *g* = 0 (or *g* = 1) corresponds to the original cascade (or niche) model, respectively. We find that, indeed, the impact of degree heterogeneity on stability depends on the prey contiguity *g* (figure 4*b*). In particular, degree heterogeneity is stabilizing for *g* < *g** while destabilizing for *g* > *g**, where *g** ≈ 0.85 (figure 4*c*). Interestingly, the contiguity of most empirical food webs [41,42] lies in the regime where degree heterogeneity favours community stability, suggesting that the existing degree heterogeneity in real-world food webs might have been shaped by stability (see electronic supplementary material, table S1, for the source and structural properties of each empirical food web). To test this hypothesis, for each empirical food web we generate a model network whose degree heterogeneity and prey contiguity are equal to that of the empirical network. As shown in figure 4*d* the stability of empirical food webs is well approximated by the corresponding model networks, indicating that degree heterogeneity and prey contiguity together largely determine the stability of predator–prey networks (see also electronic supplementary material, §IV.C).

In real food webs, the interaction strengths are thought to be skewed towards weak links and negative interactions are usually stronger than positive ones [43–56]. We explored the impact of skewness of interaction strengths by generating predator–prey networks with the relaxed niche model and assigning each link a weight that is randomly drawn from a heavy-tailed distribution instead of from a uniform range. As shown in electronic supplementary material, figure S12*a*, for predator–prey networks with skewed interaction strengths, the degree heterogeneity is stabilizing when the prey contiguity *g* is not too high (i.e. not higher than 0.95). We also consider the impact of the asymmetry between negative and positive interactions. To do so, we draw the weights of each pair of interactions (*M _{ij}*,

*M*) independently from a bivariate distribution [56], ensuring that the absolute values of negative interaction strengths are much larger than that of positive ones. We find that taking the asymmetry into account does not change our result shown in figure 4

_{ji}*c*(see electronic supplementary material, figure S12

*b*). It is worth noting that our finding also holds for other distribution with different values of correlation 〈

*M*〉 (see electronic supplementary material, figure S11).

_{ij}M_{ji}To our knowledge, there had previously been no explanation for the prominent empirical finding that feeding niches of predators in real food webs are close to contiguous but not completely so (i.e. the prey contiguity *g* is close to but not exactly 1). This curious aspect of food webs may be explained by our result that degree heterogeneity tends to destabilize completely interval (*g* = 1), but not near interval ( as in figure 4*c*) food webs. Note that this result cannot be obtained by a recent theoretical approximation [39] based on less systematic analyses of degree heterogeneity.

The discovery that a slight reduction of prey contiguity leads to such a qualitatively different stabilizing effect deserves further study. Our analyses (electronic supplementary material, §IV.C and figure S14) suggest that a network becomes more stable if, for each species, we rewire a fraction of its links to connect the species below its original feeding niche. By contrast, if we rewire a fraction of its links to connect the species above its original feeding niche, the network becomes more unstable and the degree heterogeneity is destabilizing. Since the niche hierarchy typically mimics a trophic hierarchy from plants at the base up through top-level carnivores, these results indicate that links that connect species more closely to the energetic base of the food web appear to be driving the stabilizing effect.

### 2.3. Impact of degree heterogeneity on the stability of ecological networks with mixed interactions of competition and mutualism

For ecological systems with a mixture of competitive and mutualistic interactions where *M _{ij}* and

*M*always have the same sign, we construct the community matrix

_{ji}*M*following a similar approach as in the case of predator–prey interactions. We find that for all multi-modal (figure 5

*a*), ER and SF (figure 5

*b*) model networks the Re(

*λ*

_{m}) grows monotonically as we increase

*ξ*, implying that degree heterogeneity tends to destabilize systems with competitive and mutualistic interactions. This is further demonstrated by the eigenvalue distribution of the community matrix

*M*(figure 5

*c*). We find that, similar to the predator–prey ecological systems, when the degree heterogeneity

*ξ*increases, not only does the eigenvalue distribution become non-uniform, but also the boundary that encompasses all the eigenvalues changes from ellipse to bow-tie shape (figure 5

*c*). Despite the alteration of the eigenvalue distribution, the real parts of the eigenvalues always expand as

*ξ*increases. This explains the monotonic impact of degree heterogeneity on the stability of ecological systems with a mixture of competitive and mutualistic interactions. Note that for symmetric networks with only mutualistic interactions the degree heterogeneity is also destabilizing [26].

## 3. Discussion

Here, we systematically explore, with extensive numerical and analytical calculations, the impact of degree heterogeneity on the stability of three types (predator–prey, competitive or mutualistic) of ecological networks. We do not, however, determine whether each type itself is stabilizing or not. Our results suggest that for ecological networks with random interspecific interactions or a mixture of competitive and mutualistic interactions, degree heterogeneity tends to destabilize ecological systems. For the important type of ecological networks, those involving predator–prey interactions (e.g. food webs), the impact of degree heterogeneity on stability depends on the prey contiguity. Degree heterogeneity tends to destabilize predator–prey systems with completely contiguous feeding niches while stabilizing more realistically structured systems with less than complete contiguity. This helps explain why real food webs are incompletely contiguous and have high degree heterogeneity while the stabilizing effect of breaks in contiguity that connect species more closely to biotic energy producers helps explain the mechanism behind the stabilizing effect of degree heterogeneity.

The structure of ecological networks has been recognized as one key ingredient contributing to the coexistence between high complexity and stability in real ecological systems [4,6,7,35]. Our results demonstrate that, depending on the type of interspecific interactions, degree heterogeneity of ecological networks has fundamentally different impacts on the community stability. This broadens previous results [12,14] suggesting that different types of ecological networks are restricted to assuming different architectures. Given the very few studies [26,39] on effects of degree heterogeneity on ecological network stability, our findings offer a novel perspective on how nature produces diverse, yet stable ecological systems [1]. Moreover, since we use a model-independent framework and linear stability analysis, our findings are not limited to ecological networks, but instead hold for any system of differential equations resting at an equilibrium point. For example, our results could be relevant in economic [57,58], organizational [59] and biological [60] systems where competition, cooperation and consumer–resource interactions occur.

Our results also raise several open questions and opportunities for future work. For example, there could be uncertainties in the interaction strengths of ecological network. Hence it is worth exploring the effect of the interaction strength uncertainties on the stability of degree-heterogeneous ecological networks by applying the pseudospectra of matrices [61–64]. Moreover, we used the relaxed niche model to demonstrating the effect of prey contiguity on local stability. It would be interesting to explore and compare the effect of prey contiguity in other models with varying contiguity [65,66].

## 4. Material and methods

To construct the community matrix *M*, we first generate the network using different models with a tuneable degree heterogeneity (electronic supplementary material, §II) and then assign the strength of each edge according to its interaction type [14] (see also electronic supplementary material, §I). In the following, we briefly describe the network models used in our study.

### 4.1. Multi-modal networks

A *q*-modal network has only *q* different node degrees [36]. The number of nodes with degree *k*_{(h)} is *S*/*q*, where *h* = 1,2, … , *q* and *S* is the network size. To generate a *q*-modal network, we select two different nodes *i* and *j* with probabilities *k _{i}*/(

*S*〈

*k*〉) and

*k*/(

_{j}*S〈k〉*), respectively, where 〈

*k*〉 is the mean degree, and add an edge between them unless such edge exists already. This edge-adding process repeats until there are

*S〈k〉*/2 edges in the network. To tune the degree heterogeneity, we fix the mean degree 〈

*k*〉 and change the variance of the

*q*degrees.

### 4.2. Erdős–Rényi networks

An ER random network is generated using the following steps [17]: (1) start with *S* isolated nodes; (2) select a node pair, and generate a random number between 0 and 1. If the random number exceeds *p*, where , connect the selected node pair with an edge, otherwise leave them disconnected; (3) repeat step 2 for each of the *S*(*S*−1)/2 node pairs. The node degrees of an ER network follow a binomial distribution with mean degree 〈*k*〉 = *pS* and second moment . Hence its degree heterogeneity is . Consequently, for a network with *S* nodes and mean degree 〈*k*〉, its degree heterogeneity is uniquely determined.

### 4.3. Scale-free networks

The node degrees of a SF network follow a power-law distribution and its degree heterogeneity is tuneable by changing the value of the exponent *γ* [27]. To generate SF networks with various degree heterogeneity, we employ the static method [67]: (1) start with *S* isolated nodes and assign each node *i* a weight , where 0 < *α* < 1; (2) select two different nodes *i* and *j* with probabilities *η _{i}*/(

*S〈η〉*) and

*η*/(

_{j}*S〈η〉*), respectively, where 〈

*η*〉 is the average value of the nodes' weights, and added an edge between them unless such edge exists already; (3) repeat step 2 until

_{i}*〈k〉S*/2 edges have been added into the network. The process yields a network with the power-law degree distribution where the exponent is

*γ*= 1 + 1/

*α*. Hence, we can tune its degree heterogeneity by changing the value of

*α*.

### 4.4. Cascade model

In real food webs, directed cycles appear much less common than food chains [68], meaning that the species are ordered in a natural cascade or hierarchy such that species tend to prey on species below it and tend to be preyed on by those above it in the hierarchy. The original cascade model [37] generated networks with strict hierarchical structure devoid of cycles. Since it cannot generate networks with tuneable degree heterogeneity, we resolved this issue by employing a recent variant [69] of cascade model: first, order the species with an index *i* = 1,2, … , *S*, and assign each species *i* a weight ; then, add directed edges using the static method described above. The difference is that, in step 2 of the static method, we add a directed edge from species *i* to *j* only when *i* < *j*, ensuring the predation hierarchy. Note that the networks generated with this process have skewed, though not strictly power-law, degree distributions (electronic supplementary material, figure S5), and their degree heterogeneity is tuneable with varying *α* (see electronic supplementary material, §II.D)*.*

### 4.5. Niche models

The original niche model [38] generates food webs with complete prey contiguity. To explore the impact of degree heterogeneity on the stability of food webs with varying values of prey contiguity, we employ the relaxed niche model [42] which modifies the original niche model by adding a parameter to control the niche widths relative to their maximum possible widths. When *g* = 1, the niche widths are at their narrowest and the model is identical to the original niche model. As *g* decreases towards zero, feeding ranges are widened and species that fall within the niches have lower probability of being consumed, so that non-interval networks can occur. When *g* = 0, niches are as wide as possible and the relaxed niche model mimics the cascade model. The feeding range of each species is drawn from a skewed distribution so that the degree heterogeneity can be tuned by changing the skewness of the feeding range distribution (electronic supplementary material, §II.E).

## Data accessibility

The data and code that support the findings of this study can be downloaded at: http://www.gangyanlab.com/Code&Data/EcoStability.zip.

## Authors' contributions

Y.-Y.L. conceived the project. G.Y., N.D.M. and Y.-Y.L. designed and performed the research, analysed the results and wrote the manuscript.

## Competing interests

We declare we have no competing financial interests.

## Funding

Y.-Y.L. gratefully acknowledges the support from the John Templeton Foundation (award no. 51977). N.D.M. acknowledges the support from the US National Science Foundation and the Department of Energy. G.Y. acknowledges the support from the Thousand Young Talents Program in China.

## Acknowledgement

We thank Hai-jun Zhou, R. J. Williams, Ingo Scholtes and Thilo Gross for valuable discussions.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3793447.

- Received March 12, 2017.
- Accepted May 24, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.