## Structured Abstract

### Introduction

Several major developments in theoretical ecology have relied on either dynamical stability or numerical simulations, but oftentimes, they have found contradictory results. This is partly a result of not rigorously checking either the assumption that a steady state is feasible—meaning, all species have constant and positive abundances—or the dependence of results to model parameterization. Here, we extend the concept of structural stability to community ecology in order to account for these two problems. Specifically, we studied the set of conditions leading to the stable coexistence of all species within a community. This shifts the question from asking whether we can find a feasible equilibrium point for a fixed set of parameter values, to asking how large is the range of parameter values that are compatible with the stable coexistence of all species.

### Rationale

We begin by disentangling the conditions of global stability from the conditions of feasibility of a steady state in ecological systems. To quantify the domain of stable coexistence, we first find its center (the structural vector of intrinsic growth rates). Next, we determine the boundaries of such a domain by quantifying the amount of variation from the structural vector tolerated before one species goes extinct. Through this two-step approach, we disentangle the effects of the size of the feasibility domain from how close a solution is to its boundary, which is at the heart of previous contradictory results. We illustrate our method by exploring how the observed architecture of mutualistic networks between plants and their pollinators or seed dispersers affects their domain of stable coexistence.

### Results

First, we determined the network architecture that maximizes the structural stability of mutualistic systems. This corresponds to networks with a maximal level of nestedness, a small trade-off between the number and intensity of interactions a species has, and a high level of mutualistic strength within the constraints of global stability. Second, we found that the large majority of observed mutualistic networks are close to this optimum network architecture, maximizing the range of parameters that are compatible with species coexistence.

### Conclusion

Structural stability has played a major role in several fields such as evolutionary developmental biology, in which it has brought the view that some morphological structures are more common than others because they are compatible with a wider range of developmental conditions. In community ecology, structural stability is the sort of framework needed to study the consequences of global environmental change—by definition, large and directional—on species coexistence. Structural stability will serve to assess both the range of variability a given community can withstand and why some community patterns are more widespread than others.

## A structural approach to species interactions

What determines the stability of ecological networks? Rohr *et al.* devised a conceptual approach to study interactions between species that emphasizes the role of network structure (see the Perspective by Pawar). Using the example of mutualistic networks of communities of plants and their pollinator species, they show how the structure of networks can determine the persistence of the interactions. Network structures and architectures observed in nature intrinsically match the most stable solution. This approach has promise for application to questions of ecological community stability under global change.

*Science*, this issue 10.1126/science.1253497; see also p. 383

## Abstract

In theoretical ecology, traditional studies based on dynamical stability and numerical simulations have not found a unified answer to the effect of network architecture on community persistence. Here, we introduce a mathematical framework based on the concept of structural stability to explain such a disparity of results. We investigated the range of conditions necessary for the stable coexistence of all species in mutualistic systems. We show that the apparently contradictory conclusions reached by previous studies arise as a consequence of overseeing either the necessary conditions for persistence or its dependence on model parameterization. We show that observed network architectures maximize the range of conditions for species coexistence. We discuss the applicability of structural stability to study other types of interspecific interactions.

A prevailing question in ecology [particularly since May’s seminal work in the early 1970s (*1*)] is whether, given an observed number of species and their interactions, there are ways to organize those interactions that lead to more persistent communities. Conventionally, studies addressing this question have either looked into local stability or used numerical simulations (*2*–*4*). However, these studies have not yet found a unified answer (*1*, *5*–*12*). Therefore, the current challenge is to develop a general framework in order to provide a unified assessment of the implications of the architectural patterns of the networks we observe in nature.

## Main approaches in theoretical ecology

### Dynamical stability and feasibility

Studies based on the mathematical notions of local stability, *D*-stability, and global stability have advanced our knowledge on what makes ecological communities stable. In particular, these studies explore how interaction strengths need to be distributed across species so that an assumed feasible equilibrium point can be stable (*1*–*4*, *13*–*17*). By definition, a feasible equilibrium point is that in which all species have a constant positive abundance across time. A negative abundance makes no sense biologically, and an abundance of zero would correspond to an extinct species.

The dynamical stability of a feasible equilibrium point corresponds to the conditions under which the system returns to the equilibrium point after a perturbation in species abundance. Local stability, for instance, looks at whether a system will return to an assumed feasible equilibrium after an infinitesimally small perturbation (*1*–*3*, *13*). *D*-stability, in turn, looks at the local stability of any potential feasible equilibrium that the system may have (*15*–*17*). More generally, global stability looks at the stability of any potential feasible equilibrium point after a perturbation of any given amplitude (*14*–*17*). A technical definition of these different types of dynamical stability and their relationship is provided in (*18*).

In most of these stability studies, however, a feasible equilibrium point is always assumed without rigorously studying the set of conditions allowing its existence (*5*, *14*, *15*, *19*). Yet in any given system, we can find examples in which we satisfy only one, both, or none of the feasibility and stability conditions (*3*, *16*, *17*, *19*). This means that without a proper consideration of the feasibility conditions, any conclusion for studying the stable coexistence of species is based on a system that may or may not exist (*3*, *5*, *19*).

To illustrate this point, consider the following textbook example of a two-species competition system
(1)where *N*_{1} and *N*_{2} are the abundances of species 1 and 2; β_{11} and β_{22} are their intraspecific competition strengths; β_{12} and β_{21} are their interspecific competition strengths; and α_{1} and α_{2} are their intrinsic growth rates. An equilibrium point of the system is a pair of abundances and that makes the right side of the ordinary differential equation system equal to zero.

Although the only condition necessary to guarantee the global stability of any feasible equilibrium point in this system is that the interspecific competition strengths are lower than the intraspecific ones (β_{12}β_{21} < β_{11}β_{22}), the feasibility conditions are given by and (*3*, *4*, *19*). This implies that if we set, for example, β_{11} = β_{22} = 1, β_{12} = β_{21} = 0.5, α_{1} = 1 , and α_{2} = 2, we fulfill the stability condition but not the feasibility condition, whereas if we set β_{11} = β_{22} = 0.5, β_{12} = β_{21} = 1, and α_{1} = α_{2} = 1, we can satisfy the feasibility condition but not the stability one. To have a stable and feasible equilibrium point, we need to set, for instance, β_{11} = β_{22} = 1, β_{12} = β_{21} = 0.5, and α_{1} = α_{2} = 1 (a graphical illustration is provided in Fig. 1).

The example above confirms the importance of verifying both the stability and the feasibility conditions of the equilibrium point when analyzing the stable coexistence of species (*3*–*5*, *19*). Of course, we can always fine-tune the parameter values of intrinsic growth rates so that the system is feasible (*16*, *17*). This strategy, for example, has been used when studying the success probability of an invasive species (*20*). However, when fixing the parameter values of intrinsic growth rates, we are not anymore studying the overall effect of interspecific interactions on the stable coexistence of species. Rather, we are answering the question of how interspecific interactions increase the persistence of species for a given parameterization of intrinsic growth rates. As we will show below, this is also the core of the problem in studies that are based on arbitrary numerical simulations.

### Numerical simulations

Numerical simulations have provided an alternative and useful tool with which to explore species coexistence in large ecological systems in which analytical solutions are precluded (*3*). With this approach, one has as a prerequisite to parameterize the dynamical model, or a least to have a good estimate of the statistical distribution from which these parameters should be sampled. However, if one chooses an arbitrary parameterization without an empirical justification, any study has a high chance of being inconclusive for real ecosystems because species persistence is strongly dependent on the chosen parameterization.

To illustrate this point, we simulated the dynamics of an ecological model (*6*) with three different parameterizations of intrinsic growth rates (*21*). Additionally, these simulations were performed over an observed mutualistic network of interactions between flowering plants and their pollinators located in Hickling Norfolk, UK (table S1), a randomized version of this observed network, and the observed network without mutualistic interactions (we assume that there is only competition among plants and among animals). As shown in Fig. 2, it is possible to find a set of intrinsic growth rates so that any network that we analyze is completely persistent and, at the same time, the alternative networks are less persistent.

This observation has two important implications. First, this means that by using different parameterizations for the same dynamical model and network of interactions, one can observe from all to a few of the species surviving. Second, this means that each network has a limited range of parameter values under which all species coexist. Thus, by studying a specific parameterization, for instance, one could wrongly conclude that a random network has a greater effect on community persistence than that of an observed network, or vice versa (*10*–*12*). This sensitivity to parameter values clearly illustrates that the conclusions that arise from studies that use arbitrary values in intrinsic growth rates are not about the effects of network architecture on species coexistence, but about which network architecture maximizes species persistence for that specific parameterization.

Traditional studies focusing on either local stability or numerical simulations can lead to apparently contradictory results. Therefore, we need a different conceptual framework to unify results and seek for appropriate generalizations.

### Structural stability

Structural stability has been a general mathematical approach with which to study the behavior of dynamical systems. A system is considered to be structurally stable if any smooth change in the model itself or in the value of its parameters does not change its dynamical behavior (such as the existence of equilibrium points, limit cycles, or deterministic chaos) (*22*–*25*). In the context of ecology, an interesting behavior is the stable coexistence of species—the existence of an equilibrium point that is feasible and dynamically stable. For instance, in our previous two-species competition system there is a restricted area in the parameter space of intrinsic growth rates that leads to a globally stable and feasible solution as long as ρ < 1 (Fig. 3, white area). The higher the competition strength ρ, the larger the size of this restricted area (Fig. 3) (*19*, *26*). Therefore, a relevant question here is not only whether or not the system is structurally stable, but how large is the domain in the parameter space leading to the stable coexistence of species.

To address the above question, we recast the mathematical definition of structural stability to that in which a system is more structurally stable, the greater the area of parameter values leading to both a dynamically stable and feasible equilibrium (*27*–*29*). This means that a highly structurally stable ecological system is more likely to be stable and feasible by handling a wider range of conditions before the first species becomes extinct. Previous studies have used this approach in low-dimensional ecological systems (*3*, *19*). Yet because of its complexity, almost no study has fully developed this rigorous analysis for a system with an arbitrary number of species. An exception has been the use of structural stability to calculate an upper bound to the number of species that can coexist in a given community (*6*, *30*).

Here, we introduce this extended concept of structural stability into community ecology in order to study the extent to which network architecture—strength and organization of interspecific interactions—modulates the range of conditions compatible with the stable coexistence of species. As an empirical application of our framework, we studied the structural stability of mutualistic systems and applied it on a data set of 23 quantitative mutualistic networks (table S1). We surmise that observed network architectures increase the structural stability and in turn the likelihood of species coexistence as a function of the possible set of conditions in an ecological system. We discuss the applicability of our framework to other types of interspecific interactions in complex ecological systems.

## Structural stability of mutualistic systems

Mutualistic networks are formed by the mutually beneficial interactions between flowering plants and their pollinators or seed dispersers (*31*). These mutualistic networks have been shown to share a nested architectural pattern (*32*). This nested architecture means that typically, the mutualistic interactions of specialist species are proper subsets of the interactions of more generalist species (*32*). Although it has been repeatedly shown that this nested architecture may arise from a combination of life history and complementarity constraints among species (*32*–*35*), the effect of this nested architecture on community persistence continues to be a matter of strong debate. On the one hand, it has been shown that a nested architecture can facilitate the maintenance of species coexistence (*6*), exhibit a flexible response to environmental disturbances (*7*, *8*, *36*), and maximize total abundance (*12*). On the other hand, it has also been suggested that this nested architecture can minimize local stability (*9*), have a negative effect on community persistence (*10*), and have a low resilience to perturbations (*12*). Not surprisingly, the majority of these studies have been based on either local stability or numerical simulations with arbitrary parameterizations [but, see (*6*)].

### Model of mutualism

To study the structural stability and explain the apparently contradictory results found in studies of mutualistic networks, we first need to introduce an appropriate model describing the dynamics between and within plants and animals. We use the same set of differential equations as in (*6*). We chose these dynamics because they are simple enough to provide analytical insights and yet complex enough to incorporate key elements—such as saturating, functional responses (*37*, *38*) and interspecific competition within a guild (*6*)—recently adduced as necessary ingredients for a reasonable theoretical exploration of mutualistic interactions. Specifically, the dynamical model has the following form
(2)where the variables *P _{i}* and

*A*denote the abundance of plant and animal species

_{i}*i*, respectively. The parameters of this mutualistic system correspond to the values describing intrinsic growth rates (α

*), intra-guild competition (β*

_{i}*), the benefit received via mutualistic interactions (γ*

_{ij}*), and the saturating constant of the beneficial effect of mutualism (*

_{ij}*h*), commonly known as the handling time. Because our main focus is on mutualistic interactions, we keep as simple as possible the competitive interactions for the sake of analytical tractability. In the absence of empirical information about interspecific competition, we use a mean field approximation for the competition parameters (

*6*), where we set and (

*i*≠

*j*).

Following (*39*), the mutualistic benefit can be further disentangled by , where *y _{ij}* = 1 if species

*i*and

*j*interact and zero otherwise,

*k*is the number of interactions of species

_{i}*i*, γ

_{0}represents the level of mutualistic strength, and δ corresponds to the mutualistic trade-off. The mutualistic strength is the per capita effect of a certain species on the per capita growth rate of their mutualistic partners. The mutualistic trade-off modulates the extent to which a species that interacts with few other species does it strongly, whereas a species that interacts with many partners does it weakly. This trade-off has been justified on empirical grounds (

*40*,

*41*). The degree to which interspecific interactions

*y*are organized into a nested way can be quantified by the value of nestedness

_{ij}*N*introduced in (

*42*).

We are interested in quantifying the extent to which network architecture (the combination of mutualistic strength, mutualistic trade-off, and nestedness) modulates the set of conditions compatible with the stable coexistence of all species—the structural stability. In the next sections, we explain how this problem can be split into two parts. First, we explain how the stability conditions can be disentangled from the feasibility conditions as it has already been shown for the two-species competition system. Specifically, we show that below a critical level of mutualistic strength (), any feasible equilibrium point is granted to be globally stable. Second, we explain how network architecture modulates the domain in the parameter space of intrinsic growth rates, leading to a feasible equilibrium under the constraints of being globally stable (given by the level of mutualistic strength).

### Stability condition

We investigated the conditions in our dynamical system that any feasible equilibrium point needs to satisfy to be globally stable. To derive these conditions, we started by studying the linear Lotka-Volterra approximation (*h* = 0) of the dynamical model (Eq. 2). In this linear approximation, the model reads(3)where the matrix *B* is a two-by-two block matrix embedding all the interaction strengths.

Conveniently, the global stability of a feasible equilibrium point in this linear Lotka-Volterra model has already been studied (*14*–*17*, *43*). Particularly relevant to this work is that an interaction matrix that is Lyapunov–diagonally stable grants the global stability of any potential feasible equilibrium (*14*–*18*).

Although it is mathematically difficult to verify the condition for Lyapunov diagonal stability, it is known that for some classes of matrices, Lyapunov stability and Lyapunov diagonal stability are equivalent conditions (*44*). Symmetric matrices and *Z*-matrices (matrices whose off-diagonal elements are nonpositive) belong to those classes of equivalent matrices. Our interaction strength matrix *B* is either symmetric when the mutualistic trade-off is zero (δ = 0) or is a *Z*-matrix when the interspecific competition is zero (ρ = 0). This means that as long as the real parts of all eigenvalues of *B* are positive (*18*), any feasible equilibrium point is globally stable. For instance, in the case of ρ < 1 and γ_{0} = 0, the interaction matrix *B* is symmetric and Lyapunov–diagonally stable because its eigenvalues are 1 – ρ, (*S _{A}* – 1)ρ + 1, and (

*S*– 1)ρ + 1.

_{P}For ρ > 0 and δ > 0, there are no analytical results yet demonstrating that Lyapunov diagonal stability is equivalent to Lyapunov stability. However, after intensive numerical simulations we conjecture that the two main consequences of Lyapunov diagonal stability hold (*45*). Specifically, we state the following conjectures: Conjecture 1: If *B* is Lyapunov-stable, then *B* is *D*-stable. Conjecture 2: If *B* is Lyapunov-stable, then any feasible equilibrium is globally stable.

We found that for any given mutualistic trade-off and interspecific competition, the higher the level of mutualistic strength, the smaller the maximum real part of the eigenvalues of *B* (*45*). This means that there is a critical value of mutualistic strength () so that above this level, the matrix *B* is not any more Lyapunov-stable. To compute , we need only to find the critical value of γ_{0} at which the real part of one of the eigenvalues of the interaction-strength matrix reaches zero (*45*). This implies that at least below this critical value γ_{0} < , any feasible equilibrium is granted to be locally and globally stable according to conjectures 1 and 2, respectively. We can also grant the global stability of matrix *B* by the condition of being positive definite, which is even stronger than Lyapunov diagonal stability (*14*). However, this condition imposes stronger constraints on the critical value of mutualistic strength than does Lyapunov stability (*39*).

Last, we studied the stability conditions for the nonlinear Lotka-Volterra system (Eq. 2). Although the theory has been developed for the linear Lotka-Volterra system, it can be extended to the nonlinear dynamical system. To grant the stability of any feasible equilibrium (*P _{i}* > 0 and

*A*> 0 for all

_{i}*i*) in the nonlinear system, we need to show that the above stability conditions hold on the following two-by-two block matrix (

*14*,

*43*)(4)

*B*differs from

_{nl}*B*only in the off-diagonal block with a decreased mutualistic strength. This implies that the critical value of mutualistic strength for the nonlinear Lotka-Volterra system is larger than or equal to the critical value for the linear system (

*45*). Therefore, the critical value derived from the linear Lotka-Volterra system (from the matrix

*B*) is already a sufficient condition to grant the global stability of any feasible equilibrium in the nonlinear case. However, this does not imply that above this critical value of mutualistic strength, a feasible equilibrium is unstable. In fact, when the mutualistic-interaction terms are saturated (

*h*> 0), it is possible to have feasible and locally stable equilibria for any level of mutualistic strength (

*39*,

*45*).

### Feasibility condition

We highlight that for any interaction strength matrix *B*, whether it is stable or not, it is always possible to find a set of intrinsic growth rates so that the system is feasible (Fig. 2). To find this set of values, we need only to choose a feasible equilibrium point so that the abundance of all species is greater than zero ( and ) and find the vector of intrinsic growth rates so that the right side of Eq. 2 is equal to zero: and . This reconfirms that the stability and feasibility conditions are different and that they need to be rigorously verified when studying the stable coexistence of species (*3*, *16*, *17*, *19*). This also highlights that the relevant question is not whether we can find a feasible equilibrium point, but how large is the domain of intrinsic growth rates leading to a feasible and stable equilibrium point. We call this domain the feasibility domain.

Because the parameter space of intrinsic growth rates is substantially large (**R*** ^{S}*, where

*S*is the total number of species), an exhaustive numerical search of the feasibility domain is impossible. However, we can analytically estimate the center of this domain with what we call the structural vector of intrinsic growth rates. For example, in the two-species competition system of Fig. 4A the structural vector is the vector (in red), which is in the center of the domain leading to feasibility of the equilibrium point (white region). Any vector of intrinsic growth rates collinear to the structural vector guarantees the feasibility of the equilibrium point—that is, guarantees species coexistence. Because the structural vector is the center of the feasibility domain, then it is also the vector that can tolerate the strongest deviation before leaving the feasibility domain—that is, before having at least one species going extinct.

In mutualistic systems, we need to find one structural vector for animals and another for plants. These structural vectors are the set of intrinsic growth rates that allow the strongest perturbations before leaving the feasibility domain. To find these structural vectors, we had to transform the interaction-strength matrix *B* to an effective competition framework (*45*). This results in an effective competition matrix for plants and a different one for animals (*6*), in which these matrices represent respectively the apparent competition among plants and among animals, once taking into account the indirect effect via their mutualistic partners. With a nonzero mutualistic trade-off (δ > 0), the effective competition matrices are nonsymmetric, and in order to find the structural vectors, we have to use the singular decomposition approach—a generalization of the eigenvalue decomposition. This results in a left and a right structural vector for plants and for animals in the effective competition framework. Last, we need to move back from the effective competition framework in order to obtain a left and right vector for plants (and ) and animals (and ) in the observed mutualistic framework. The full derivation is provided in the supplementary materials (*45*).

Once we locate the center of the feasibility domain with the structural vectors, we can approximate the boundaries of this domain by quantifying the amount of variation from the structural vectors allowed by the system before having any of the species going extinct—that is, before losing the feasibility of the system. To quantify this amount, we introduced proportional random perturbations to the structural vectors, numerically generated the new equilibrium points (*21*), and measured the angle or the deviation between the structural vectors and the perturbed vectors (a graphical example is provided in Fig. 4A). The deviation from the structural vectors is quantified, for the plants, by , where and are, respectively, the angles between α^{(}^{P}^{)} and and between α^{(}^{P}^{)} and . α^{(}^{P}^{)} is any perturbed vector of intrinsic growth rates of plants. The deviation from the structural vector of animals is computed similarly.

As shown in Fig. 4B, the larger is the deviation of the perturbed intrinsic growth rates from the structural vectors, and the lower is the persistence of the community as defined by the fraction of surviving species. This confirms that there is a restricted domain of intrinsic growth rates centered on the structural vectors compatible with the stable coexistence of species. The greater the tolerated deviation from the structural vectors within which all species coexist, the greater the feasibility domain, and in turn the greater the structural stability of the system.

As shown in Fig. 4B, the greater the deviation of the perturbed intrinsic growth rates from the structural vectors, the lower is the fraction of surviving species. This confirms that there is a restricted domain of intrinsic growth rates centered on the structural vectors compatible with the stable coexistence of species. The greater the tolerated deviation from the structural vectors within which all species coexist, the greater the feasibility domain, and in turn the greater the structural stability of the system.

### Network architecture and structural stability

To investigate the extent to which network architecture modulates the structural stability of mutualistic systems, we explored the combination of alternative network architectures (combinations of nestedness, mutualistic strength, and mutualistic trade-off) and their corresponding feasibility domains.

To explore these combinations, for each observed mutualistic network (table S1) we obtained 250 different model-generated nested architectures by using an exhaustive resampling model (*46*) that preserves the number of species and the expected number of interactions (*45*). Theoretically, nestedness ranges from 0 to 1 (*42*). However, if one imposes architectural constraints such as preserving the number of species and interactions, the effective range of nestedness that the network can exhibit may be smaller (*45*). Additionally, each individual model-generated nested architecture is combined with different levels of mutualistic trade-off δ and mutualistic strength γ_{0}. For the mutualistic trade-off, we explored values δ ∈ [0, ..., 1.5] with steps of 0.05 that allow us to explore sublinear, linear, and superlinear trade-offs. The case δ = 0 is equivalent to the soft mean field approximation studied in (*6*). For each combination of network of interactions and mutualistic trade-off, there is a specific critical value in the level of mutualism strength γ_{0} up to which any feasible equilibrium is globally stable. This critical value is dependent on the mutualistic trade-off and nestedness. However, the mean mutualistic strength shows no pattern as a function of mutualistic trade-off and nestedness (*45*). Therefore, we explored values of with steps of 0.05 and calculated the new generated mean mutualistic strengths. This produced a total of 250 × 589 different network architectures (nestedness, mutualistic trade-off, and mean mutualistic strength) for each observed mutualistic network.

We quantified how the structural stability (feasibility domain) is modulated by these alternative network architectures in the following way. First, we computed the structural vectors of intrinsic growth rates that grant the existence of a feasible equilibrium of each alternative network architecture. Second, we introduced proportional random perturbations to the structural vectors of intrinsic growth rates and measured the angle or deviation (η_{(}_{A}_{)}, η_{(}_{P}_{)}) between the structural vectors and the perturbed vectors. Third, we simulated species abundance using the mutualistic model of (*6*) and the perturbed growth rates as intrinsic growth rate parameter values (*21*). These deviations lead to parameter domains from all to a few species surviving (Fig. 4).

Last, we quantified the extent to which network architecture modulates structural stability by looking at the association of community persistence with network architecture parameters, once taking into account the effect of intrinsic growth rates. Specifically, we studied this association using the partial fitted values from a binomial regression (*47*) of the fraction of surviving species on nestedness (*N*), mean mutualistic strength (), and mutualistic trade-off (δ), while controlling for the deviations from the structural vectors of intrinsic growth rates (η_{(}_{A}_{)}, η_{(}_{P}_{)}). The full description of this binomial regression and the calculation of partial fitted values are provided in (*48*). These partial fitted values are the contribution of network architecture to the logit of the probability of species persistence, and in turn, these values are positively proportional to the size of the feasibility domain.

### Results

We analyzed each observed mutualistic network independently because network architecture is constrained to the properties of each mutualistic system (*11*). For a given pollination system located in the KwaZulu-Natal region of South Africa, the extent to which its network architecture modulates structural stability is shown in Fig. 5. Specifically, the partial fitted values are plotted as a function of network architecture. As shown in Fig. 5A, not all architectural combinations have the same structural stability. In particular, the architectures that maximize structural stability (reddish/darker regions) correspond to the following properties: (i) a maximal level of nestedness, (ii) a small (sublinear) mutualistic trade-off, and (iii) a high level of mutualistic strength within the constraint of any feasible solution being globally stable (*49*).

A similar pattern is present in all 23 observed mutualistic networks (*45*). For instance, using three different levels of interspecific competition (ρ = 0.2, 0.4, 0.6), we always find that structural stability is positively associated with nestedness and mutualistic strength (*45*). Similarly, structural stability is always associated with the mutualistic trade-off by a quadratic function, leading quite often to an optimal value for maximizing structural stability (*45*). These findings reveal that under the given characterization of interspecific competition, there is a general pattern of network architecture that increases the structural stability of mutualistic systems.

Yet, one question remains to be answered: Is the network architecture that we observe in nature close to the maximum feasibility domain of parameter space under which species coexist? To answer this question, we compared the observed network architecture with theoretical predictions. To extract the observed network architecture, we computed the observed nestedness from the observed binary interaction matrices (table S1) following (*42*). The observed mutualistic trade-off δ is estimated from the observed number of visits of pollinators or fruits consumed by seed-dispersers to flowering plants (*41*, *50*, *51*). The full details on how to compute the observed trade-off is provided in (*52*). Because there is no empirical data on the relationship between competition and mutualistic strength that could allow us to extract the observed mutualistic strength γ_{0}, our results on nestedness and mutualistic trade-off are calculated across different levels of mean mutualistic strength.

As shown in Fig. 5, B to D, the observed network (blue solid lines) of the mutualistic system located in the grassland asclepiads of South Africa actually appears to have an architecture close to the one that maximizes the feasibility domain under which species coexist (reddish/darker region). To formally quantify the degree to which each observed network architecture is maximizing the set of conditions under which species coexist, we compared the net effect of the observed network architecture on structural stability against the maximum possible net effect. The maximum net effect is calculated in three steps.

First, as outlined in the previous section, we computed the partial fitted values of the effect of alternative network architectures on species persistence (*48*). Second, we extracted the range of nestedness allowed by the network given the number of species and interactions in the system (*45*). Third, we computed the maximum net effect of network architecture on structural stability by finding the difference between the maximum and minimum partial fitted values within the allowed range of nestedness and mutualistic trade-off between δ ∈ [0, ..., 1.5]. All the observed mutualistic trade-offs have values between δ ∈ [0, ..., 1.5]. Last, the net effect of the observed network architecture on structural stability corresponds to the difference between the partial fitted values for the observed architecture and the minimum partial fitted values extracted in the third step described above.

Looking across different levels of mean mutualistic strength, in the majority of cases (18 out of 23, *P* = 0.004, binomial test) the observed network architectures induce more than half the value of the maximum net effect on structural stability (Fig. 6, red solid line). These findings reveal that observed network architectures tend to maximize the range of parameter space—structural stability—for species coexistence.

## Structural stability of systems with other interaction types

In this section, we explain how our structural stability framework can be applied to other types of interspecific interactions in complex ecological systems. We first explain how structural stability can be applied to competitive interactions. We proceed by discussing how this competitive approach can be used to study trophic interactions in food webs.

For a competition system with an arbitrary number of species, we can assume a standard set of dynamical equations given by , where α* _{i}* > 0 is the intrinsic growth rate, β

*> 0 is the competition interaction strength, and*

_{ij}*N*is the abundance of species

_{i}*i*. Recall that the Lyapunov diagonal stability of the interaction matrix β would imply the global stability of any feasible equilibrium point. However, in nonsymmetric competition matrices Lyapunov stability does not always imply Lyapunov diagonal stability (

*53*). This establishes that we should work with a restricted class of competition matrices such as the ones derived from the niche space of (

*54*). Indeed, it has been demonstrated that this class of competition matrices are Lyapunov diagonally stable and that this stability is independent of the number of species (

*55*). For a competition system with a symmetric interaction-strength matrix, the structural vector is equal to its leading eigenvector. For other appropriate classes of matrices, we can compute the structural vectors in the same way as we did with the effective competition matrices of our mutualistic model and numerically simulate the feasibility domain of the competition system. In general, following this approach we can verify that the lower the average interspecific competition, the higher is the feasibility domain, and in turn, the higher is the structural stability of the competition system.

In the case of predator-prey interactions in food webs, so far there is no analytical work demonstrating the conditions for a Lyapunov–diagonally stable system and how this is linked to its Lyapunov stability. Moreover, the computation of the structural vector of an antagonistic system is not a straightforward task. However, we may have a first insight about how the network architecture of antagonistic systems modulates their structural stability by transforming a two-trophic–level food web into a competition system among predators. Using this transformation, we are able to verify that the higher the compartmentalization of a food web, then the higher is its structural stability. There is no universal rule to study the structural stability of complex ecological systems. Each type of interaction poses their own challenges as a function of their specific population dynamics.

## Discussion

We have investigated the extent to which different network architectures of mutualistic systems can provide a wider range of conditions under which species coexist. This research question is completely different from the question of which network architectures are aligned to a fixed set of conditions. Previous numerical analyses based on arbitrary parameterizations were indirectly asking the latter, and previous studies based on local stability were not rigorously verifying the actual coexistence of species. Of course, if there is a good empirical or scientific reason to use a specific parameterization, then we should take advantage of this. However, because the set of conditions present in a community can be constantly changing because of stochasticity, adaptive mechanisms, or global environmental change, we believe that understanding which network architectures can increase the structural stability of a community becomes a relevant question. Indeed, this is a question much more aligned with the challenge of assessing the consequences of global environmental change—by definition, directional and large—than with the alternative framework of linear stability, which focuses on the responses of a steady state to infinitesimally small perturbations.

We advocate structural stability as an integrative approach to provide a general assessment of the implications of network architecture across ecological systems. Our findings show that many of the observed mutualistic network architectures tend to maximize the domain of parameter space under which species coexist. This means that in mutualistic systems, having both a nested network architecture and a small mutualistic trade-off is one of the most favorable structures for community persistence. Our predictions could be tested experimentally by exploring whether communities with an observed network architecture that maximizes structural stability stand higher values of perturbation. Similarly, our results open up new questions, such as what the reported associations between network architecture and structural stability tell us about the evolutionary processes and pressures occurring in ecological systems.

Although the framework of structural stability has not been as dominant in theoretical ecology as has the concept of local stability, it has a long tradition in other fields of research (*29*). For example, structural stability has been key in evolutionary developmental biology to articulate the view of evolution as the modification of a conserved developmental program (*27*, *28*). Thus, some morphological structures are much more common than others because they are compatible with a wider range of developmental conditions. This provided a more mechanistic understanding of the generation of form and shape through evolution (*56*) than that provided by a historical, functionalist view. We believe ecology can also benefit from this structuralist view. The analogous question here would assess whether the invariance of network architecture across diverse environmental and biotic conditions is due to the fact that such a network structure is the one increasing the likelihood of species coexistence in an ever-changing world.

## Supplementary Materials

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Any matrix
*B*is called Lyapunov-stable if the real parts of all its eigenvalues are positive, meaning that a feasible equilibrium at which all species have the same abundance is at least locally stable. A matrix*B*is called*D*-stable if*DB*is a Lyapunov-stable matrix for any strictly positive diagonal matrix*D*.*D*-stability is a stronger condition than Lyapunov stability in the sense that it grants the local stability of any feasible equilibrium (*15*). In addition, a matrix*B*is called Lyapunov diagonally stable if there exists a strictly positive diagonal matrix*D*so that*DB*+*B*is a Lyapunov-stable matrix. This notion of stability is even stronger than^{t}D*D*-stability in the sense that it grants not only the local stability of any feasible equilibrium point, but also its global stability (*14*). A Lyapunov–diagonally stable matrix has all of its principal minors positive. Also, its stable equilibrium (which may be only partially feasible; some species may have an abundance of zero) is specific (*57*). This means that there is no alternative stable state for a given parameterization of intrinsic growth rates (*55*). We have chosen the convention of having a minus sign in front of the interaction-strength matrix*B*. This implies that*B*has to be positive Lyapunov-stable (the real parts of all eigenvalues have to be strictly positive) so that the equilibrium point of species abundance is locally stable. This same convention applies for*D*-stability and Lyapunov–diagonal stability. - ↵
- ↵
- ↵All simulations were performed by integrating the system of ordinary differential equations by using the function ode45 of Matlab. Species are considered to have gone extinct when their abundance is lower than 100 times the machine precision. Although we present the results of our simulations with a fixed value of ρ = 0.2 and
*h*= 0.1, our results are robust to this choice as long as ρ < 1 (*45*). In our dynamical system, the units of the parameters are abundance for*N*, 1/time for α and*h*, and 1/(time · biomass) for β and γ. Abundance and time can be any unit (such as biomass and years), as long as they are constant for all species. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Materials and methods are available as supplementary materials on
*Science*Online. - ↵
- ↵
- ↵For a given network, we look at the association of the fraction of surviving species with the deviation from the structural vector of plants η
_{(}_{P}_{)}and animals η_{(}_{A}_{)}, nestedness*N*, mutualistic trade-off δ, and mean level of mutualistic strength using the following binomial generalized linear model (*47*): . Obviously, at a level of mutualism of zero (γ = 0), nestedness and mutualistic trade-off cannot influence the probability of a species to survive. We have to include an interaction between the mean level of mutualism and the nestedness and mututalistic trade-off. We have also included a quadratic term in order to take into account potential nonlinear effects of nestedness, mutualistic trade-off, and mean level of mutualistic strength. The effect of network architecture is confirmed by the significant likelihood ratio between the full model and a null model without such a network effect (*P*< 0.001) for all the observed empirical networks. The effect of network architecture on structural stability can be quantified by the partial fitted values defined as follows: , where are the fitted parameters corresponding to the terms to , respectively. - ↵
- ↵
- ↵
- ↵To estimate the mutualistic trade-off, δ, from an empirical point of view, we proceed as follows. First, there is available data on the frequency of interactions. Thus,
*q*is the observed number of visits of animal species_{ij}*j*on plant species*i*. This quantity has been proven to be the best surrogate of per capita effects of one species on another (γ_{ij}^{P}and γ_{ij}^{A}) (*41*). Second, the networks provide information on the number of species one species interacts with (its degree or generalization level;*k*_{i}^{A}and*k*_{i}^{P}). Following (*41*,*50*,*51*), the generalization level of a species has been found to be proportional to its abundance at equilibrium. Thus, the division of the total number of visits by the product of the degree of plants and animals can be assumed to be proportional to the interaction strengths. Mathematically, we obtain two equations, one for the effect of the animals on the plants and vice versa: (*q*)/(_{ij}*k*_{i}^{P}*k*_{j}^{A}) ∝ γ_{ij}^{P}and (*q*)/(_{ij}*k*_{i}^{A}*k*_{j}^{P}) ∝ γ_{ij}^{A}. Introducing the explicit dependence between interaction strength and trade-of—for example, γ_{ij}^{P}= γ_{0}/(*k*_{i}^{P})^{δ}and (*q*)/(_{ij}*k*_{i}^{P}*k*_{j}^{A}) ∝(γ_{0})/[(*k*_{i}^{P})^{δ}] —we obtain (*q*)/(_{ij}*k*_{i}^{P}*k*_{j}^{A}) ∝ (γ_{0})/[(*k*_{i}^{P})^{δ}] and (*q*)/(_{ij}*k*_{i}^{A}*k*_{j}^{P}) ∝ (γ_{0})/[(*k*_{i}^{A})^{δ}]. In order to estimate the value of δ, we can just take the logarithm on both sides of the previous equations for the data excluding the zeroes. Then, δ is simply given by the slope of the following linear regressions: and , where*a*and^{P}*a*are the intercepts for plants and animals, respectively. These two regressions are performed simultaneously by lumping together the data set. The intercept for the effect of animals on plants (^{A}*a*) may not be the same as the intercept for the effect of the plants on the animals (^{P}*a*).^{A} - ↵The competition matrix given by is Lyapunov-stable but not
*D*-stable. Because Lyapunov diagonal stability implies*D*-stability, this matrix is also not Lyapunov–diagonally stable. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**We thank U. Bastolla, L.-F. Bersier, V. Dakos, A. Ferrera, M. A. Fortuna, L. J. Gilarranz, S. Kéfi, J. Lever, B. Luque, A. M. Neutel, A. Pascual-García, D. B. Stouffer, and J. Tylianakis for insightful discussions. We thank S. Baigent for pointing out the counter-example in (*53*). Funding was provided by the European Research Council through an Advanced Grant (J.B.) and FP7-REGPOT-2010-1 program under project 264125 EcoGenes (R.P.R.). The data are publicly available at www.web-of-life.es.