Interspecific competition shapes the structural stability of mutualistic networks

the community is both feasible and stable — of ecological communities in which both mutualistic interactions between plants and pollinators and competitive interactions among plants and among pollinators are present. Using the structure of 50 real networks for mutualistic interactions, combined with analytical and numerical analyses, we show that the structure of the competitive network radically alters the necessary conditions for species coexistence in these communities. Our mathematical framework also allows to accurately characterize the structural stability of these systems. Moreover, we introduce a new metric that accurately links the network structures of competitive and mutualistic interactions to species coexistence. Our results highlight the joint role of the structures of different interaction types to understand the stability of ecological communities and facilitate the analysis of similar natural and artificial systems in which mutualism and competition coexist.


Introduction
Species rarely live in isolation, but constantly interact with other species with different interaction types, such as predation, competition and mutualism [38].Mutualism, in which different species interact for their mutual benefit, is ubiquitous in terrestrial ecosystems [4].Examples include, but are not limited to, plants receiving effective pollination or seed-dispersion by offering rewards of nutrients to their visiting animals, plants gaining resistance to insect herbivores by offering nutrients and shelter to fungi or ants, and leguminous plants obtaining nitrogen by rewarding nitrogen-fixing bacteria.
Interspecific competition, where species within the same guild compete for shared mutualistic partners, is one of the identified costs when progressing from two species mutualism to species-rich mutualism [30].This had already been reported in field experiments of mutualistic systems of plant and pollinators by Charles Robertson [31] in 1895, which was then followed by extensive studies in mutualisms of ants and plants, and of parrots and plants [12,10,15,13,30,35,7,28].Intra-guild competition among plants is specifically recognized when common pollinators frequently visit pollen-or nectar-rich species while reducing or avoiding the visitation to less-rewarding plants [27,3,39,18,17].Extremely attractive species can become dominant in the long run (e.g.well-identified representatives of Lythrum salicaria [8] and Impatiens glandulifera [11]), possibly threatening the persistence of lessrewarding species.
Competition among pollinators is likewise relevant for the functioning of plant-pollinator communities.In fact, for certain hummingbird species, interspecific competition may be as important as mutualistic interactions in shaping the evolution of species that coexist in particular geographical areas [19,9].More specifically, as a consequence of competition, hummingbird species may experience morphological changes, such as in bill length, which improve the pollination efficiency or to expand the diversity of pollinated flowers [19,9].A large amount of evidence suggests that morphological specialization is also an evolutionary strategy to avoid or reduce inter-and intra-specific competition in communities of bumblebees [40].Yet the effects of competition are perhaps more perceptible on shorter time scales.For instance, empirical studies reveal that competition for floral resources significantly alters the feeding performance and the harvest of nectar in pollinator communities when foreigner bees are inserted in a given environment [23,33,34]; the presence of invader bees is also likely to affect the availability of food and nest sites, which in turn may undermine the abundance of native species [16,34,41].The immediate rearrangement of mutualistic interactions by some pollinators after the intentional removal of competing species is another notable example of the key role that the shared use of resources plays in the dynamics of ecological communities [30].
Main theories of biodiversity, however, have largely ignored the diversity of interaction types that link species in nature and have instead focused on a small subset of well studied interactions, such as predation, competition and mutualism, each of them being typically studied in isolation from the others [24].Decades of studies on these interactions have shown that ecological networks have a specific architecture, which plays a key role for their dynamics and stability (e.g.[14,29,5,36]).Mutualistic networks -such as plant-pollinator networks, for example -have attracted increasing attention in the ecological literature in the last decades [22].These networks have been shown to be highly nested, with more specialist species interacting with a subset of the species that interact with more generalist species [5], which has been suggested to contribute to the maintenance of species diversity [37].
Despite the tremendous contribution of these previous studies to the understanding of the link between the structure and dynamics of ecological networks, the lack of studies explicitly incorporating the diversity of interaction types has hindered the advancement of our understanding of the factors that drive the number of species that can coexist in a given community, one of the oldest questions in ecology [46].
In previous studies, competitive interactions between plants sharing pollinators were modeled as an all-to-all connectivity pattern; that is, in a plant-pollinator scenario, all pollinators are considered to compete equally for plants, and all plants are assumed to compete for pollinators regardless of the heterogeneous organization of the mutualistic interactions.However, in real ecosystems, it seems unlikely that species would compete equally for shared partners without accounting for the structure of the mutualistic links.To which extent homogeneous competition can predict species coexistence and how variation in competition alters mechanisms that maintain biodiversity remains unknown.Given the empirical evidence of the structure of mutualistic networks and the lack of empirical knowledge about the structure of the associated competitive networks, it is of utmost importance to understand the effect of different competitive network structures on the community stability of mutualistic networks.
Here, we investigate the effect of different assumptions regarding the assignment of competitive links among plants and among pollinators (and thereby the structure of the competitive networks) on the stability of ecological communities including competitive and mutualistic interactions.Real network structures were used for the mutualistic part of the ecological communities.More specifically, we investigate the "feasible area", i.e. the set of conditions (parameters) under which all species coexist and have a positive abundance.We develop a framework that predicts accurately the boundaries of the feasible area.Furthermore, we find that different competitive network structures yield significantly different feasibility patterns, showing that the structure of competitive interactions does have strong implications for the species diversity of multilayer networks including competition and mutualism.

Population dynamics of competitive-mutualistic networks
To study the impact of interspecific competition on communities persistence, we begin by describing dynamics between plant and pollinator species.We consider a plant-pollinator system consisting of a set A of N A animal species that interact mutualistically with a set P of N P plant species, denoting the total biodiversity by N = N P + N A .The mutualistic interactions are fully encoded in a N P × N A bipartite matrix K, where K ij = 1 if plant species i is pollinated by pollinator species j, and 0 otherwise.Each plant (resp.animal) species is characterized by the abundance s P i (resp.s A i ), whose dynamics depend on the intrinsic growth rate α P i (resp.α A i ) and on the influence of competitive and

Mutualistic Network
P l a n t

Full mean field
A n i m a l P l a n t A n i m a l mutualistic interactions as follows:

Soft mean field Multilayer
where i = 1, ..., N P , and M P i = k∈A K ik s A k is the total abundance of pollinators interacting with plant i.The first term in the right-hand side of the above equations corresponds to the intrinsic growth of each species; the second and fourth are also identical in all models and correspond to the intra-species competition and mutualistic interactions, respectively; and β refers to the intensity of intra-specific competition.The intensities of inter-specific competition and mutualism are denoted by β 0 and γ 0 , respectively.Parameter h, known as the handling time, imposes a nonlinear saturation effect on mutualism.What distinguishes the models in Eqs. ( 13)-( 15) is the definition of the third term, which accounts for the intra-guild competition.We highlight schematically the differences between each competition scenario with Fig. 1.In the "full mean field" model, all species within a layer compete equally with each other, i.e., all plants compete with all other plants with the same intensity [see Fig. 1(a)], and all pollinators compete with other pollinators with the same intensity, irrespective of the mutualistic links existing between plants and pollinators (this is the approach used in previous studies; see, e.g., [6]).In the "soft mean-field" model, the homogeneous competition assumption is relaxed by placing a competitive link between plants i and j (A P ij = 1) only if they share at least one pollinator.The intra-guild competition links in this model are encoded in matrices A P and A A .Notice that, in this scenario, the competitive links are not weighted (only present or absent).Finally, in the "weighted competition" model, the intra-guild competitive links have the same structure as in the "soft mean-field" model but are now weighted by the abundance of the shared partner [see Fig. 1(c)].This weight is set via matrix W P , whose elements are given by W P ij = k∈A K ik K jk s A k [equations for the pollinator abundances s A i follow mutatis mutandis from Eqs. ( 13)-( 15)].Therefore, in the weighted scenario, the higher the abundance of mutualistic partners, the stronger the competition among plants which have common mutualistic connections.Notice also that the intra-guild competition terms in Eq. ( 15) are asymmetric, since the biomass of shared species is normalized by the total biomass of mutualistic partners j∈P,i =j W P ij /M P i .In other words, two plant species i and j perceive the competition with one another differently according to the importance of their shared pollinators in relation with their respective total abundance of pollinators.
In the Supplemental Material, we provide an exact and thorough bifurcation analysis of the toy network depicted in Fig. 1.In the next section we present an analytical calculation of the solutions of Eqs. ( 13)-( 15) for arbitrary networks.

Structural stability conditions
Our goal here is to derive an analytical calculation for the feasible equilibrium solution of the nonlinear population dynamics in Eqs. ( 13)- (15), i.e., the solution that corresponds to the maximum biodiversity (s P,A i > 0 ∀i).It is argued [32] that a specific parameterization can be inconclusive for empirical networks due to the strong dependence of species coexistence on parameterization.Accordingly, we investigate a range of parameter values termed feasible area under which all species coexist.Henceforth, we refer to the "feasible area" as the region in the space spanned by parameters β 0 and γ 0 in which all species have positive abundances at equilibrium.

Condition for the full and soft mean field competitions
We first address the solution for the full and soft mean-field models.By applying a linear approximation to the nonlinear mutualism term in Eqs. ( 13) and ( 14), we rewrite the population dynamics in a matrix form as where A P,A and M P,A are the matrices that set the competition and mutualism interactions, respectively; and s P,A are the vectors containing the individual abundances of plant and pollinators, respectively.In the full mean-field model, we have A P ij = 1 for i = j, while in the soft mean-field model A P ij = Θ k∈A K ik K jk , where Θ(•) is the heaviside function.The elements of matrices M P,A are obtained by applying the Taylor expansion for the mutualistic term in Eqs. ( 13)- (14), that is: expanded around a point (M 0 ) i close to a fixed point that is challenging to obtain without a prior knowledge.We sometimes omit the subscript (M 0 ) i when there is no ambiguity.After substituting both the competitive and mutualistic terms, a feasible equilibrium (s P i , s A i > 0) can thus be obtained by solving the following linear equation where the vector c = h . Without a prior knowledge on the fixed points of the system, the linearizion of the system near a fixed point appears to be challenging or even unfeasible.We approach this challenge by analyzing the interplay between the mutualistic interactions and the intra-guild competition, which separately lead to abundance gain and abundance loss at equilibrium.When the mutualistic strength is equal to the competition strength, the species abundance on average follows s i = α i β i .Assuming the average abundance s A k = s i for all the animal species pollinating plant i, we linearize the nonlinear population dynamics at M P 0 i = d P i s i for each plant i, where d P i = k K ik denoting the number of animals pollinating plant i.The fixed point for animal species M A 0 is approximated similarly.Equation ( 6) provides an approximated solution for the abundances of general networks in the full and soft-mean field competition.Notice that Eq. ( 6) provides the equilibrium solution of the system, but does not guarantee feasibility.In order to estimate the feasible area, one needs to solve Eq. ( 6) for different parameters seeking solutions satisfying s P,A i > 0. Equation 6 has the numerical advantage that it allows one to scan the parameter space of ecological networks and delineate the feasible area much more quickly than by evolving the original dynamics.Differently from the solutions assuming h = 0 in [32], Equation ( 6) is applicable to any real h ≥ 0, thus enabling the investigation of various mutualistic regimes.A point (β 0 , γ 0 ) is colored in blue if all species survive with a positive abundance in the stationary regime of the simulations of Eqs. ( 13)- (15) for that parameter choice.Parameters: α i = 1∀i, and β = 5.Diagrams in the upper panels have h = 0, while simulations in the lower panels are for h = 0.1.Grid size: 100 × 100.Solid lines are obtained by solving Eq. ( 6).
In Fig. 2(a,d) and (b,e) we compare the analytical predictions provided by Eq. ( 6) with direct simulations of the systems in Eqs. ( 13)-( 15) over the parameter space spanned by competition and mutualism strengths, β 0 and γ 0 , respectively.As it is seen, the analytical prediction delineates the boundaries of the feasible area with remarkable accuracy for the soft mean-field model, for both h = 0 and h = 0.1.In the full mean-field model, reasonable precision is achieved for h = 0, while for h = 0.1 the matching between numerical and theoretical boundaries is lost as competition strength increases.In the Supplemental Material, we show that the approximate solution of Eq. ( 6) is successful in predicting the feasible area for several real networks under the full mean-field and soft mean-field competition scenarios.

Condition for the weighted competition
To analyze the impact of heterogeneity in the competitive strengths among intra-guild species, we derive conditions under which species coexist for the dynamical model with weighted competition.The weighted competition is reduced to the homogeneous competitive strength (soft-mean field competition) only when any pair of intra-guild species depend on and share exactly the same set of mutualistic partners, corresponding to complete-like bipartite mutualistic networks.
In the dynamic model for weighted competition [Eq.( 15)], the inter-specific competition interactions are weighted by the relative importance of shared resources in a nonlinear form.To tackle the nonlinear inter-specific competition, we harness the microscopic perspective of intra-guild competition induced by a single mutualistic plant-pollinator interaction.When plant i is pollinated by an animal k (i.e., K ik = 1), the inter-specific competition between plant i and the other plants j that are pollinated by animal k reads j∈P,j =i K ik K T kj s P j .Summing over all the pollinators k that pollinate plant i yields the total inter-specific competition k j∈P,j =i K ik K T kj s P j s A k for that plant i.Armed with the view of a single mutualistic interaction, the nonlinear competition term in Eq. ( 15) can be rewritten as j∈P,j =i Analytical estimates can be obtained by reigning in the weighted competition terms.This can be accomplished by using the mediant inequality (see SI material): Equation ( 8) allows reducing the complexity of species abundance from two-guilds into a singleguild species abundance.After the complexity reduction, the weighted competition structure can be encoded in a competition matrix ÃP (resp.ÃA ) for plant (resp.animal) species, in analogy to the competition matrix, incorporated in the unweighted adjacency matrix A P , for the soft mean-field case.Additionally, Eq. ( 8) establishes lower and upper bounds for the competition term.
Seeking to find an accurate estimate for the competition term, it is appropriate to consider that plant i competes with plant j mediated by sharing the animal k, whose degree is, over all pollinators of plant i, the closest to the local average of pollinated plants.Specifically, the element of the competition matrix ÃP can be written as where ) is the number of animals (plants) with which plant (animal) i interacts, meaning that the weighted competition is approximated by an effective mutualistic partner k whose degree is the closest to the average competitive species per mutualistic interaction.
Combining Eq. ( 9) with the corresponding mutualistic term into an expression analogous to Eq. ( 6), the feasible solution for weighted competition model is eventually obtained by The analytical prediction of Eq. ( 10) for the weighted model is also well confirmed by numerical simulations on real mutualistic networks [Fig.2(c) and (f)].Besides checking the validity of our calculations, Fig. 2 also allows us to highlight the marked differences in the dynamics yielded by the three competition models.Notice, in particular, how strongly the shape of the feasible area changes from the full mean-field to the soft mean-field and then to the weighted competition case: the mere inclusion of heterogeneity in the competition term in Eq. ( 13) shifts the region of occurrence of feasible states from strong competition (full mean-field) to weak competition (soft mean-field, weighted).
Another noteworthy difference between soft-mean field and weighted competition is that the former allows a wider feasible region for low competition and high mutualism, while the latter is favored by modest values of both competition and mutualism.These results clearly demonstrate how crucial the structure of the intra-guild competition networks is for the species coexistence and the diversity of ecological communities.

Stability-complexity paradox
After investigating the impact of different choices of the structure of intra-guild competition networks on biodiversity, let us now study how the feasible area relates with the topological properties of the multilayer ecological networks.To do so, we consider a set of 50 real plant-pollinator networks retrieved from the Web of Life platform [43] and, for each network in the database, we evolve Eqs. ( 13)- (15) and calculate the corresponding feasible area for each intra-guild competition scenario.We look at different components of complexity independently, namely species diversity, the average interaction strength and connectance.Following May [46], we define the average interaction strength as the average of the off-diagonal elements of the Jacobian matrices of systems ( 13)- (15); and connectance is defined as the density of non-zero values of the Jacobian matrices (see Supplemental Material).As shown in Fig. 3, the feasible area correlates positively with connectance and interaction strength in the full mean-field, soft mean-field and weighted competition scenarios.In other words, the higher the number of interactions among species and the stronger their intensity, the more likely that the ecological community exhibits feasible states.This result actually points back to the long debated diversity-stability paradox initiated by May [46].In his work, May proved that the probability of facing stable states converges almost certainly to zero for sufficiently large communities.In mathematical terms, suppose that species i and j interact with probability C and via an interaction strength J ij , which is a random variable with mean E(J ij ) = 0 and with variance given by Var(J ij ) = σ 2 .Under such conditions, and setting the self-interaction terms as constants, J ii = −d ∀i, May proved that the dynamical system ds/dt = Js is almost surely stable if in particular with d = 1, i.e. the condition for which the leading eigenvalue of J is negative.Consequently, the increases in size, connectivity and interaction strength favor the dynamical destabilization of the system.This finding triggered the aforementioned paradox because it seems in contradiction with the high diversity of species observed in many natural communities [26].At first sight, our results seem to violate the stability-diversity paradox, since they show a positive correlation between feasible area and connectance, and between feasible area and interaction strength, suggesting that more connected -and thereby more complex-networks tend to have a higher feasibility area, meaning that such systems would tend to be more stable [Fig.3, panels (a) and (c)].However, this conclusion is reached by looking at the different aspects of complexity independently, whereas a closer inspection suggests that some of these aspects are related.For example, incorporating the third element of May's relation and inspecting the size of the networks, we realize that the networks with a high connectance are also those that are the smallest networks in the database (Fig. 3, size of the dots stands for number of species in the corresponding networks).What is more, by checking manually the networks with connectance values 0.12, one notices that they correspond to almost fully connected bipartite matrices and, hence, to almost fully connected Jacobian matrices as well.Thus, the networks with the highest possible values for the feasible area in Fig. 3 are in fact the networks with the less "complex" structure in the data set, in consonance with the notion that complexity, as quantified by this trade-off between system size and connectance, tends to destabilize ecological communities.
We now look more specifically at May's criteria combining the three network metrics and explore how the size of the feasibility area is related to that criteria.Our goal here is to check whether there is a clear relation between the feasible area and the expressions in Eqs. ( 13)- (15).However, it is noteworthy that May's condition [Eq.(11)] is related to the probability that the ecological system is stable for a particular set of parameters, whereas the feasible area results from a sum over different parameter combinations.Therefore, to appropriately evaluate how the dependencies of the feasible area on network properties relate with Eq. ( 11), we define the following quantity: where • corresponds to an average over the Jacobian matrix's elements, and • (β,γ) stands for the average over the parameters β and γ considered in Fig. 3.The first term in Eq. ( 25), J ii (β,γ) , is the average taken over the diagonal elements, since, contrarily to the random model considered by May, the diagonal elements of the Jacobian, J P,A ii , are not constant, but rather are heterogeneously distributed over the diagonal (see Supplemental Material); the term σ(J ij ) (β,γ) is the average standard deviation of the off-diagonal values of J.In practical terms, variable C May defined in Eq. ( 25) quantifies how distant a given network is from the critical point established by May's stability condition, meaning how stable it is according to that criteria.For each real network considered, we numerically calculate the Jacobian elements evaluated at the stationary points, which in turn are obtained by evolving Eqs ( 13)- (15) numerically.We visualize the dependence of the feasible area on C May in Fig. 4. Interestingly, the size of the feasibility area does not correlate with May's criteria for the full and soft mean field scenarios, but it is very strongly correlated to it in the case of the weighted competition scenario, despite the fact that May's criteria was formulated for much more idealized systems (Fig. 4).The patterns we see in Fig. 4 also agree with the scatter plots in Fig. 3; i.e., there is a noticeable correlation between network size and the value of the coefficient C May .In particular, we observe that large networks tend to exhibit low values for C May , while "less complex" networks (in terms of size and connectivity) have higher values of C May .
This raises the question of why the weighted competition model adheres with May's criteria better than the full and soft mean-field scenarios.The answer lies in the expression of the Jacobian terms of the different models in Eqs. ( 13)- (15).For the weighted competition scenario, the off-diagonal Jacobian elements J ij are proportional to the terms 1/s P,A i and 1/(s P,A i ) 2 (see the SI material).This makes the elements J ij to narrowly peak around an average value, thus making the standard deviation σ(J ij ) less significant than the average of the diagonal terms J ii (see SI material) and thereby creating a positive correlation between feasible area and C May .Comparing the models of heterogeneous competitive interactions (i.e.soft mean-field and weighted in Fig. 4), one notices that the introduction of weights in the competitive structure acts as a stabilizing factor in the dynamics, in the sense that networks tend to exhibit a wider feasible area in the weighted scenario.The interplay of interaction types among species, a well-defined structure, and interacting strength shows a markedly driving force in stabilizing the system that is highly predictable and consistent with May's criterion.Interestingly, the feasible area of the weighted competition scenario correlates with May's criterion as initially formulated by May for random matrices, but no significant relation was observed with other, more recent stability criteria formulated for random matrix models that incorporate features from predator-prey and mutualistic interactions [44] (see SI material).

Discussion
We investigated to which extent the incorporation of intra-guild competition alters the maintenance of biodiversity in mutualistic systems.Compared to a scenario where all species from a guild homogeneously compete with each other, as commonly assumed in the literature, heterogeneous competition leads to markedly different patterns for the feasibility area of plant-pollinator networks.Without sufficient empirical data about intra-guild competition in plant-pollinator communities, deriving the structure of the competitive links from the observed mutualistic interactions enables us to theoretically explore how the structure of different intra-guild competition networks affects the structural stability of mutualistic ecological communities.
Our results show that previously identified implications can be restricted to homogeneous competition and cannot be readily generalized to heterogeneous competition.Specifically, we found that feasibility patterns are dramatically modified when competitive interactions become heterogeneous.This finding suggests that a series of important conclusions regarding the dynamics of ecological networks might have been overlooked given that theoretical models have been traditionally studied under the assumption of homogeneous intra-guild competitive interactions.Therefore, getting information on the structure of competitive networks in mutualistic systems is key to better understand what constrains the assembly of mutualistic communities during the dynamical coevolutionary process, which could be an important driving force of coevolution [21].
Finally, by investigating the feasible area in terms of global network properties, we found that smaller and less connected networks exhibited larger regions sustaining maximum biodiversity.Interestingly, this result agrees with the long-standing May's stability-diversity paradox, which states that complex systems are more prone to be destabilized as their size, connectance and mean interaction strength increase.Indeed, we have verified that the more structurally stable networks turned out to be the "less complex" ones according to May's criterion [46].Our results therefore show that the complexity introduced in the model by the weighted competition scenario yields a phenomenology which is well predicted by a condition originally derived for random systems [Fig.4(c)], whereas it is not the case for the other competition scenarios.Importantly, the analysis on feasible solutions performed here is not limited to the specific equations studied here, but can also be extended to ecological networks with different types of interactions, such as facilitation [42], and predator-prey models [20].
funded by FAPESP (Grant No. 2013/07375-0).Y.M. acknowledges partial support from the Government of Aragon and FEDER funds, Spain through grant E36-20R (FENOL), by MINECO and FEDER funds (FIS2017-87519-P), and by Intesa Sanpaolo Innovation Center.The funders had no role in study design, data collection, and analysis, decision to publish, or preparation of the manuscript.

Supplementary Material 6 Exact solution for the structural stability of the minimal model
For the sake of clarity, let us briefly reintroduce the models discussed in the main text.We consider plant-pollinator systems composed of N A animal species, which interact with N P plant species.The matrix encoding the mutualistic is the N P × N A bipartite matrix K, whose elements are K ij = 1 if plant species i is pollinated by animal species j.The abundance of the ith plant is defined as s P i (s A i , analogously for the animal species), and its time-depend dynamics is governed by the following equations: where α P,A i are the intrinsic growth rate; parameters β and β 0 stand for the intra-and inter-specific competition strength, respectively; and γ 0 is the mutualism strength.Variable M P i is given by where k denotes the indexes belonging to the pollinators set.Matrix A P ij of the "Soft mean-field competition" scenario encodes the competitive connections within the plantguild: A P ij = 1 if plants i and j share at least one pollinator, and 0 otherwise.For the "Weighted competition scenario", the elements of matrix W P are defined as is the abundance of the kth pollinator.The equations for ds A k /dt are obtained by interchanging the labels P and A in the equations and by defining the terms A A , M A and W A accordingly.
In this section, we present the exact solution for the structural stability of the models in Eqs.( 13)-( 15) considering the toy network shown in Fig. 1 in the main text.In order to determine which regions of parameter space are actually accessible by a dynamical system one needs to consider both the stability and the feasibility of the equilibrium solution.The feasibility alone is insufficient as numerical simulations only allow the stable regions to be accessed.The feasibility conditions, that is, positive abundances for all species, should therefore be considered as prerequisite for any equilibrium solution, after which the stability of the equilibrium can be established.
For the minimal model illustrated in Fig. 1 in the main text, the governing equation for the population dynamics can be written in a matrix form as We first linearize the system ( 16) around the feasible equilibrium s * .Let us denote a small perturbation around the equilibrium as u = s − s * .Substituting s = u + s * to Eq. ( 16) yields where the entries of matrix K eff are given by where I is the identity matrix.By computing the eigenvalues of matrix K eff we can retrieve the stability of the equilibrium.The stability changes when one of the eigenvalues crosses zero.Such points can be found by setting the determinant of matrix K eff equal to zero.In order to find the stability region of the feasible equilibrium of Eq. ( 16), with all s i > 0, we perform a stability analysis.That is we do not consider the equilibrium with one of the species having zero (or negative) abundance.Taking matrix K of the full mean field as an example, feasible equilibrium is given by the following expressions with With the feasible equilibrium s * and the competition matrix K, we construct the matrix K eff .Setting the determinant of matrix K eff to zero, we find the expression for the curves with a vanishing eigenvalue.In the parameter space of β 0 and γ 0 , the stability status changes at each time we across those zero lines.From det K eff = 0, we obtain 15 − 2β 0 which are the boundary curves for stability regions.Analogously, we describe the exact solution for the feasible equilibrium and stability conditions for the soft mean-field and weighted competition models as shown in Figure 5.
For the case of nonzero Holling term, h = 0, similar feasibility and stability analysis can be established by applying the linear approximation of the mutualistic terms.Figure 6 illustrates the analytical results of feasibility and stability conditions for the full mean-field, soft mean-field and weighted competition scenarios, where panel (a) is for h = 0, and panel (b) is for h = 0.3.Figure 6 suggests an altered feasible area modulated by inter-specific competition for the mutualistic network shown in Figure 5.

Analytical prediction for the population dynamics of the weighted competition model
To accurately predict the population dynamics of the weighted competition scenario, we derive the weighted competition matrix.For an observed mutualistic network, and a given initial condition where Soft mean-field: Feasibility conditions Stability conditions Stability conditions Figure 5: Feasibility and stability conditions for soft-mean field and the weighted competition model for the minimal network.
one realization of population dynamics is sufficient, we derive the competition matrix as where ) is the number of animals (plants) with which plant (animal) i interacts.Term quantifies the local average of plant competitors mediated by sharing pollinators over number of pollinators d P i of plant i.Recall that in the weighted scenario the competition intensity between two plants (pollinators) is directly proportional to the relative abundance of shared pollinators (plants).Therefore, the more abundant the mutualistic partners, the fiercer the intra-guild competition among the members that have common mutualistic connections.
Figures 7 and 8 compare the simulation results of different networks with the analytical predictions.As can be seen, for h = 0 (Fig. 7), the matching between the numerical results and our calculations is remarkable, especially for the full and soft mean-field scenarios.Notice, in particular, how the detailed contour of the feasible area in the latter scenario is almost exactly captured by the analytical curves.As we can see in Fig. 7, the feasible area in the weighted competition case is slightly underestimated by the theory, but nonetheless the general trends of the boundaries are reproduced very closely by the solid curves.For the case h = 0.1 in Fig. 8 the analytical results accurately predict the feasible area for a large range of mutualistic strength γ 0 , with exceptions at a high value of mutualism due to the saturation effect.All in all, our analytical results work very well for a variety of real plant-pollinator networks.

Evaluation of analytical predictions on parameterization of population dynamics
To investigate the performance of our analytical predictions of the feasible area, in this section we present results from extensive simulations considering several parameter combinations.Specifically, we compare the analytical solutions derived in the previous section and in the main text for real plant-pollinator networks under different choices for the intraspecific competition β and the Holling term h.

Variation in the intraspecific competition β
Figures 9 and 10 show the numerical and analytical results when we increase or decrease the intensity of intraspecific competition β i for the observed mutualistic networks M-PL-048 and M-PL-016.For the population dynamics with soft mean-field competition, the analytical results accurately predict the feasible area for both the increase and decrease of intraspecific competition β.For the weighted competition scenario, the linear approximation captures the changing behavior of species coexistence.

Weak mutualism h 1
Compared to the majority of the theories in ecology dealing with either weak mutualism (h 1) or strong mutualism h = 0, our theoretical framework fills the gap and covers the full range (0 < h < 1).We proceed by presenting the feasible area in the extreme cases of h = 0 and h 1.For a strong mutualism h = 0, the dynamical model is reduced to Lotka-Volterra equations with type II functional responses.Strong mutualism without saturation is often unstable leading to indefinite and unbounded growth of species which is argued to be biologically unrealistic due to the environmental constraints like carrying capacity [47].However, strong mutualism at the same time overpowers the inter-specific competition resulting in a cut-off of the mutualism intensity γ 0 .
For a weak mutualism h 1, the mutualism term is saturated to a constant 1/h, i.e., lim When the mutualism is saturated, the dynamical model is reduced to a linear model, and the feasible equilibrium is obtained by solving where u is the all-one vector.For the weighted competition case, a linear approximation on the matrices ÃP and ÃA is applied.The mutualism intensity γ 0 is disentangled from affecting the feasible area (which is confirmed by the numerical results in Fig. 11).Moreover, saturated mutualism, in turn, imposes limitations on the intensity β 0 of competition.Because without the compensation of mutualism benefit, the intra-specific competition might be overly strong and eventually make species go extinct.

Structural stability and mutualistic network properties
In this section, we provide complementary results on the correlation between feasible area and network architectures for dynamical models of full mean field, soft mean field and weighted competition.We test two more network properties, namely, the variance of the degree sequence and the ratio of interdegree to intra-degree for both plants and animals.In addition, we test the maximum interspecific competition β 0 before losing any species and its relation to network architectures when the mutualism is saturated.
Here, we provide the correlation between feasible area and network architectures for all the considered dynamical models.In particular, we consider the regimes of strong (h = 0) and saturated mutualism (h = 0).For the strong mutualism regime (Fig. 12), global network properties, such as the number of species and connectance, correlate strongly and negatively with the feasible area, with goodness-of-fit of R 2 = 0.83 and R 2 = 0.81 for the soft mean-field competition.The maximum degree shows a less strong correlation with the feasible area, and neither does nestedness show a clear dependence on the same quantity.The same trends are observed for other competition scenarios in the strong mutualism regime (see Fig. 12) as well as for the dynamics with saturated mutualism (Fig. 13).
10 Analysis on the diversity-stability relation May [46] has established a stability-diversity relation for a large ecological system whose stability is characterized by dx dt = Jx, indicating that a large system is less stable.In May's assumption, each element of the Jacobian matrix J is equally likely to be positive or negative, having an absolute magnitude chosen from a random distribution with zero mean and standard deviation α.Matrix J can also be written in the form J = B − I. Based on Theorem 1 for the largest eigenvalue of a random matrix [45], the system is stable if σ where σ is the standard deviation of the random variable from which the off-diagonal elements of the Jacobian matrix J take value; N is the size of the network, and C denotes the network connectance [45].
where σ is the standard variation.
Our goal now is to examine how the stability condition in Eq. 23 relates with the feasible area of real plant-pollinator networks.However, as discussed in the main text, Eq. ( 23) is associated with the probability that the system is stable given a particular set of parameters; the feasible area, on the other hand, is defined for a set of parameters.Therefore, in order to establish a relation between May's condition and the feasible area, we define where • is an average over the Jacobian matrix's elements, and • (β 0 ,γ 0 ) corresponds to the average over certain ranges of parameters β 0 and γ 0 .The first term in Eq. ( 25), J ii (β 0 ,γ 0 ) , is the average taken over the diagonal elements, since, contrarily to the random model considered by May, the diagonal elements of the Jacobian, J P,A ii , are not constant, but rather are heterogeneously distributed over the diagonal (see Appendix); the term σ(J ij ) (β 0 ,γ 0 ) is the average standard deviation of the off-diagonal values of J.
The stability criterion by May is derived considering purely random interactions among the elements in a complex system, i.e., the Jacobian matrix is a random matrix without any constraints on its elements.More recently, Allesina and Tang (AT) [44] generalized May's stability analysis for random matrix models that incorporate aspects of predator-prey, mutualisc and competitive interactions [44].Seeking to verify how these more complex stability criteria relate with the feasible area of networks, we define where C AT,Pred corresponds to the stability condition of random predator-prey matrices [44], and C AT,Mutu is the analogous quantity for random matrix models that emulate the interplay between competitive and mutualistic relationships [44].Figure 14 shows the observed correlation of May's and Allesina and Tang's stability with the feasible area for full mean-field, soft mean-field and weighted competition scenarios.Interestingly, we do not observe a significant dependence of the feasible area on conditions ( 25)-( 27) for full and soft mean-field scenario.A pattern does emerge in the weighted scenario, for which we have a positive correlation between feasible area and all conditions; however, it is interesting to note that the strongest correlation occurs for C May .The surprise about this result resides in the fact that C May is the condition derived for the simplest random matrix model.At first sight, one would expect to observe a more significant dependence of the feasible area on C AT,Mutu , which accounts for competitive and mutualistic networks, but what we have is the opposite: the most complex model (weighted competition scenario) adheres best with the condition derived from May's stability criteria.As we argue in the main text, the explanation for the latter results lies in the expression of the Jacobian elements (see Appendix).In the weighted competition scenario, the off-diagonal Jacobian elements J ij depend on terms 1/s P,A i and 1/(s P,A i ) 2 .Since the abundance values are generally less than 1, the elements J ij end up being narrowly peaked around an average value, thus making the standard deviation σ(J ij ) to be less significant than the average off-diagonal terms J ii .To exemplify this phenomenon, in Figs. 15 and 16 we show the distribution of elements in Jacobian matrix for networks MPL-010 and MPL-048, respectively.
In order to get further insights into the dynamics of the three competition models, let us denote the right hand side of Eq. ( 15) as f (s P i ).To determine the stability of the equilibrium point, we analyze the Jacobian matrix J which can be written as where B can be expressed in the form of block matrix as .An analogous form can be obtained for B 22 and B 21 by replacing superscript P indicating plant species to superscript A representing pollinator species and vice versa.We derive the expression of submatrices B 11 and B 12 for three cases of full mean-filed, soft mean-field and weighted competition (see also Appendix A).
(i) Case of full mean-field competition: the submatrix B 11 = β 0 J N P , where J N P is the all one matrix.Therefore, all elements in Jacobian submatrix B 11 are linearly correlated.Each element in B 12 is determined mainly by mutualistic interactions (B 12 ) iv = All the interacting pollinators of plant i have the same value in Jacobian submatrix B 12 and thus linearly correlated.
(ii) Case of soft mean-field competition: the submatrix B 11 has an element of 1 whenever there is competition which has the same value for all competed species.Each element in B 12 is determined mainly by mutualistic interactions (B 12 ) iv = 2 .(iii) Case of weighted competition: the submatrix B 11 has elements computed by For each pollinator u of plant i, the value in Jacobian submatrix is varied, in contrast to the same value of 1 in the case of full mean field and soft mean field.Each element in submatrix B 12 is computed by The second term in the right hand side of Eq. ( 31) is introduced due to mutualistic interactions and shows an analogous pattern to the full mean field case and soft mean field case.However, the first term uniquely appears in the weighted competition scenario.In addition, the value is varied for different pollinators v, determined by the number of introduced plant competitions j K iv K vj (s * ) P j mediated by sharing a common pollinator v. Weighted competition reduces the correlation between elements in each row of the Jacobian submatrix B 12 and, therefore, shows a well agreement with May's stability criteria, which is built upon the assumption of independence among elements of the Jacobian matrix.In the Appendix A we provide the complete expressions for the elements of the Jacobian matrix.

A Expressions for the Jacobian matrix elements
The Jacobian matrices of the models considered in the main text can be expressed as follows In the sequel we write the expression for each competition scenario.A.1 Full mean-field competition model ∂ ṡP i ∂s P j ≡ J P ij = −s P i β 0 , i, j ∈ P.
∂ ṡi ∂s A.2 Soft mean-field competition model A.3 Weighted competition model  Figure 9: Analytical prediction (black curve) for the soft mean-field, considering different values for the inter-specific competition strength β A i = β P i = β.The system in Eq. ( 14) was numerically integrated with the Heun's method, considering total simulation time T = 2000 and time step dt = 0.01.The simulation result (shaded area) is obtained with parameters α P i = α A i = 1 and h = 0.1.

Figure 1 :
Figure 1: (a) Illustration of the minimal mutualistic network that distinguishes the (b) full meanfield, (c) soft mean-field and (d) weighted competition scenarios.For any other network with fewer nodes than shown in (a), the competition interactions of the models in (b)-(d) become identical.For the full mean-field, each layer is a complete, unweighted graph representing an entire inter-specific competition with the same magnitude.For the soft mean-field, each layer is an unweighted graph with connections (representing inter-specifies competition) only between pollinators (plants) who share plants (pollinators).For the weighted competition scenario, each layer is a weighted graph with weights representing the strength of inter-specific competition.

Figure 2 :
Figure 2: Feasible area patterns illustrated here for a real plant-pollinator network (MPL-16) from the Web of Life platform [43], in the (left) full mean-field, (center) soft mean-field and (right) weighted scenarios.A point (β 0 , γ 0 ) is colored in blue if all species survive with a positive abundance in the stationary regime of the simulations of Eqs.(13)-(15) for that parameter choice.Parameters: α i = 1∀i, and β = 5.Diagrams in the upper panels have h = 0, while simulations in the lower panels are for h = 0.1.Grid size: 100 × 100.Solid lines are obtained by solving Eq. (6).

Figure 3 :
Figure 3: Feasible area as a function of different network metrics for full mean field [(a) and (d)], soft mean-field [(b) and (e)], and weighted [(c) and (f)] competition scenarios.Each dot corresponds to a real mutualistic network from the Web of Life platform [43], with the size of dot proportional to the network size.Parameters common in all panels: α i = 1∀i, h = 0.1 and β = 5.Feasible area was calculated considering the same ranges of β 0 and γ 0 shown in Fig. 2.

Theorem 1 .
Let M be a random N × N real and symmetric matrix where elements m ij = m ji are independent random variables.Assume that these random variables possess a common mean E[m ij ] = 0 and common variance Var[m ij ] = σ 2 and E[m ii ] = µ.Then the largest eigenvalue is upper bounded by max |λ 1 | ≤ 2σ

Figure 6 :
Figure 6: Feasible domain of the minimal model for full mean-field, soft mean-field and the weighted competition.Panel (a) shows the theoretical results for h = 0, and (b) shows the approximation result for h = 0.3.Other parameters are taken as α P 1 = α A 2 = 1, β P = β A = 5.

Figure 7 :
Figure 7: Feasible area patterns for several networks from the Web of Life platform [43], in the (left) full mean-field, (center) soft mean-field and (right) weighted scenarios.Parameters h = 0, β = 5, and α P,A i = 1 ∀i.

Figure 8 :
Figure 8: Feasible area patterns for several networks from the Web of Life platform [43], in the (left) full mean-field, (center) soft mean-field and (right) weighted scenarios.Parameters h = 0, β = 5, and α P,A i = 1 ∀i.

Figure 10 :Figure 11 :
Figure 10: Analytical prediction (black curve) for the weighted competition scenario, considering different values for the inter-specific competition strength β A i = β P i = β.The systems in Eqs.(15) was numerically integrated with the Heun's method, considering total simulation time T = 2000 and time step dt = 0.01.The simulation result (shaded area) is obtained with parameters α P i = α A i = 1 and h = 0.1.
B = (B 11 ) N P ×N P (B 12 ) N P ×N A (B 21 ) N A ×N P (B 22 ) N A ×N A (29) Each element in B 11 is calculated by (B 11 ) iu = is the abundance at equilibrium.Each element in B 12 is calculated by (B 12 ) iv = *