Interspecific competition shapes the structural stability of mutualistic networks
Xiangrong Wang, Thomas Peron, Johan L. A. Dubbeldam, Sonia Kèfi, Yamir Moreno
IInterspecific competition shapes the structural stability of mutualisticnetworks
Xiangrong Wang † , Thomas Peron † , Johan L. A. Dubbeldam , Sonia K´efi , and YamirMoreno Institute of Future Networks, Southern University of Science and Technology, People’sRepublic of China Research Center of Networks and Communications, Peng Cheng Laboratory, People’sRepublic of China Institute of Mathematics and Computer Science, University of S˜ao Paulo, 13566-590 S˜aoCarlos, SP, Brazil Faculty of Electrical Engineering, Mathematics and Computer Science, P.O. Box 5031, 2600GA Delft, The Netherlands ISEM, CNRS, Univ. Montpellier, IRD, EPHE, Montpellier, France Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA Institute for Biocomputation and Physics of Complex Systems (BIFI), University ofZaragoza, Zaragoza, Spain Department of Theoretical Physics, Faculty of Sciences University of Zaragoza, Zaragoza,Spain Institute for Scientific Interchange Foundation, Turin, Italy † These authors contributed equally to this work.26 January 2021
Abstract
Mutualistic networks have attracted increasing attention in the ecological literature in the lastdecades as they play a key role in the maintenance of biodiversity. Here, we develop an analyticalframework to study the structural stability of these networks including both mutualistic and com-petitive interactions. Analytical and numerical analyses show that the structure of the competitivenetwork fundamentally alters the necessary conditions for species coexistence in communities. Us-ing 50 real mutualistic networks, we show that when the relative importance of shared partners isincorporated via weighted competition, the feasibility area in the parameter space is highly corre-lated with May’s stability criteria and can be predicted by a functional relationship between thenumber of species, the network connectance and the average interaction strength in the community.Our work reopens a decade-long debate about the complexity-stability relationship in ecologicalcommunities, and highlights the role of the relative structures of different interaction types. a r X i v : . [ phy s i c s . s o c - ph ] F e b Introduction
Species rarely live in isolation, but constantly interact with other species with different interactiontypes, such as predation, competition and mutualism [38]. Mutualism, in which different speciesinteract for their mutual benefit, is ubiquitous in terrestrial ecosystems [4]. Examples include, but arenot limited to, plants receiving effective pollination or seed-dispersion by offering rewards of nutrientsto their visiting animals, plants gaining resistance to insect herbivores by offering nutrients and shelterto fungi or ants, and leguminous plants obtaining nitrogen by rewarding nitrogen-fixing bacteria.Interspecific competition, where species within the same guild compete for shared mutualisticpartners, is one of the identified costs when progressing from two species mutualism to species-richmutualism [30]. This had already been reported in field experiments of mutualistic systems of plantand pollinators by Charles Robertson [31] in 1895, which was then followed by extensive studies inmutualisms of ants and plants, and of parrots and plants [12, 10, 15, 13, 30, 35, 7, 28]. Intra-guildcompetition among plants is specifically recognized when common pollinators frequently visit pollen- ornectar-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 representativesof
Lythrum salicaria [8] and
Impatiens glandulifera [11]), possibly threatening the persistence of less-rewarding species.Competition among pollinators is likewise relevant for the functioning of plant-pollinator com-munities. In fact, for certain hummingbird species, interspecific competition may be as importantas mutualistic interactions in shaping the evolution of species that coexist in particular geographicalareas [19, 9]. More specifically, as a consequence of competition, hummingbird species may experiencemorphological changes, such as in bill length, which improve the pollination efficiency or to expandthe diversity of pollinated flowers [19, 9]. A large amount of evidence suggests that morphologicalspecialization is also an evolutionary strategy to avoid or reduce inter- and intra-specific competitionin communities of bumblebees [40]. Yet the effects of competition are perhaps more perceptible onshorter time scales. For instance, empirical studies reveal that competition for floral resources sig-nificantly alters the feeding performance and the harvest of nectar in pollinator communities whenforeigner bees are inserted in a given environment [23, 33, 34]; the presence of invader bees is also likelyto affect the availability of food and nest sites, which in turn may undermine the abundance of nativespecies [16, 34, 41]. The immediate rearrangement of mutualistic interactions by some pollinatorsafter the intentional removal of competing species is another notable example of the key role that theshared use of resources plays in the dynamics of ecological communities [30].Main theories of biodiversity, however, have largely ignored the diversity of interaction types thatlink species in nature and have instead focused on a small subset of well studied interactions, suchas predation, competition and mutualism, each of them being typically studied in isolation from theothers [24]. Decades of studies on these interactions have shown that ecological networks have a specificarchitecture, which plays a key role for their dynamics and stability (e.g. [14, 29, 5, 36]). Mutualisticnetworks – such as plant-pollinator networks, for example — have attracted increasing attention inthe 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 generalistspecies [5], which has been suggested to contribute to the maintenance of species diversity [37].2espite the tremendous contribution of these previous studies to the understanding of the linkbetween the structure and dynamics of ecological networks, the lack of studies explicitly incorporatingthe diversity of interaction types has hindered the advancement of our understanding of the factorsthat drive the number of species that can coexist in a given community, one of the oldest questions inecology [46].In previous studies, competitive interactions between plants sharing pollinators were modeled as anall-to-all connectivity pattern; that is, in a plant-pollinator scenario, all pollinators are considered tocompete equally for plants, and all plants are assumed to compete for pollinators regardless of the het-erogeneous organization of the mutualistic interactions. However, in real ecosystems, it seems unlikelythat species would compete equally for shared partners without accounting for the structure of themutualistic links. To which extent homogeneous competition can predict species coexistence and howvariation in competition alters mechanisms that maintain biodiversity remains unknown. Given theempirical evidence of the structure of mutualistic networks and the lack of empirical knowledge aboutthe structure of the associated competitive networks, it is of utmost importance to understand theeffect of different competitive network structures on the community stability of mutualistic networks.Here, we investigate the effect of different assumptions regarding the assignment of competitivelinks among plants and among pollinators (and thereby the structure of the competitive networks)on the stability of ecological communities including competitive and mutualistic interactions. Realnetwork 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 coexistand have a positive abundance. We develop a framework that predicts accurately the boundaries ofthe feasible area. Furthermore, we find that different competitive network structures yield significantlydifferent feasibility patterns, showing that the structure of competitive interactions does have strongimplications for the species diversity of multilayer networks including competition and mutualism.
To study the impact of interspecific competition on communities persistence, we begin by describingdynamics between plant and pollinator species. We consider a plant-pollinator system consisting of aset A of N A animal species that interact mutualistically with a set P of N P plant species, denotingthe 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 0otherwise. Each plant (resp. animal) species is characterized by the abundance s Pi (resp. s Ai ), whosedynamics depend on the intrinsic growth rate α Pi (resp. α Ai ) and on the influence of competitive and3 utualistic 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 Soft mean field
Multilayer P l a n t A n i m a l a) b)c) d) Figure 1: (a) Illustration of the minimal mutualistic network that distinguishes the (b) full mean-field, (c) soft mean-field and (d) weighted competition scenarios. For any other network with fewernodes than shown in (a), the competition interactions of the models in (b)-(d) become identical. Forthe full mean-field, each layer is a complete, unweighted graph representing an entire inter-specificcompetition with the same magnitude. For the soft mean-field, each layer is an unweighted graphwith connections (representing inter-specifies competition) only between pollinators (plants) who shareplants (pollinators). For the weighted competition scenario, each layer is a weighted graph with weightsrepresenting the strength of inter-specific competition.mutualistic interactions as follows:1 s Pi ds Pi dt = α Pi − βs Pi − β N P (cid:88) j (cid:54) = i s Pj + γ M Pi hγ M Pi (Full mean-field competition), (1)1 s Pi ds Pi dt = α Pi − βs Pi − β N P (cid:88) j =1 A Pij s Pj + γ M Pi hγ M Pi (Soft mean-field competition), (2)1 s Pi ds Pi dt = α Pi − βs Pi − β N P (cid:88) j =1 (cid:32) W Pij M Pi (cid:33) s Pj + γ M Pi hγ M Pi (Weighted competition), (3)4here i = 1 , ..., N P , and M Pi = (cid:80) k ∈A K ik s Ak is the total abundance of pollinators interacting withplant i . The first term in the right-hand side of the above equations corresponds to the intrinsicgrowth of each species; the second and fourth are also identical in all models and correspond to theintra-species competition and mutualistic interactions, respectively; and β refers to the intensity ofintra-specific competition. The intensities of inter-specific competition and mutualism are denoted by β and γ , respectively. Parameter h , known as the handling time, imposes a nonlinear saturationeffect 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 betweeneach competition scenario with Fig. 1. In the “full mean field” model, all species within a layercompete 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, irrespectiveof the mutualistic links existing between plants and pollinators (this is the approach used in previousstudies; see, e.g., [6]). In the “soft mean-field” model, the homogeneous competition assumption isrelaxed by placing a competitive link between plants i and j ( A Pij = 1) only if they share at least onepollinator. The intra-guild competition links in this model are encoded in matrices A P and A A . Noticethat, 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 Pij = (cid:80) k ∈A K ik K jk s Ak [equations forthe pollinator abundances s Ai follow mutatis mutandis from Eqs. (13)-(15)]. Therefore, in the weightedscenario, the higher the abundance of mutualistic partners, the stronger the competition among plantswhich have common mutualistic connections. Notice also that the intra-guild competition terms inEq. (15) are asymmetric, since the biomass of shared species is normalized by the total biomass ofmutualistic partners (cid:16)(cid:80) j ∈P ,i (cid:54) = j W Pij /M Pi (cid:17) . In other words, two plant species i and j perceive thecompetition with one another differently according to the importance of their shared pollinators inrelation with their respective total abundance of pollinators.In the Supplemental Material, we provide an exact and thorough bifurcation analysis of the toynetwork depicted in Fig. 1. In the next section we present an analytical calculation of the solutions ofEqs. (13)-(15) for arbitrary networks. Our goal here is to derive an analytical calculation for the feasible equilibrium solution of the nonlinearpopulation dynamics in Eqs. (13)-(15), i.e., the solution that corresponds to the maximum biodiver-sity ( s P,Ai > ∀ i ). It is argued [32] that a specific parameterization can be inconclusive for empiricalnetworks due to the strong dependence of species coexistence on parameterization. Accordingly, weinvestigate a range of parameter values termed feasible area under which all species coexist. Hence-forth, we refer to the “feasible area” as the region in the space spanned by parameters β and γ inwhich all species have positive abundances at equilibrium.5 ondition 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 approxima-tion to the nonlinear mutualism term in Eqs. (13) and (14), we rewrite the population dynamics in amatrix form as (cid:34) ds P dtds A dt (cid:35) = diag (cid:32)(cid:34) s P s A (cid:35)(cid:33) (cid:32)(cid:34) α P α A (cid:35) − (cid:32) βI + β (cid:34) A P A A (cid:35) − γ (cid:34) (cid:102) M P (cid:102) M A (cid:35)(cid:33) (cid:34) s P s A (cid:35)(cid:33) , (4)where A P,A and (cid:102) M P,A are the matrices that set the competition and mutualism interactions, re-spectively; 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 Pij = 1 for i (cid:54) = j , while in the soft mean-fieldmodel A Pij = Θ (cid:0)(cid:80) k ∈A K ik K jk (cid:1) , where Θ( · ) is the heaviside function. The elements of matrices (cid:102) M P,A are obtained by applying the Taylor expansion for the mutualistic term in Eqs. (13)-(14), that is: γ M i hγ M i = γ ( M ) i hγ ( M ) i + (cid:18) γ M i hγ M i (cid:19) (cid:48) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M i =( M ) i ( M i − ( M ) i ) , (5)expanded around a point ( M ) i close to a fixed point that is challenging to obtain without a priorknowledge. We sometimes omit the subscript ( M ) i when there is no ambiguity. After substitutingboth the competitive and mutualistic terms, a feasible equilibrium ( s Pi , s Ai >
0) can thus be obtainedby solving the following linear equation (cid:34) α P α A (cid:35) = βI + β (cid:34) A P A A (cid:35) − (cid:18) γ ( hγ ( M P ) i ) (cid:19) K diag (cid:18) γ ( hγ ( M A ) i (cid:105) ) (cid:19) K T (cid:34) s P s A (cid:35) − c, (6)where the vector c = h (cid:18)(cid:16) γ M P hγ M P (cid:17) , (cid:16) γ M A hγ M A (cid:17) (cid:19) T . Without a prior knowledge on the fixedpoints of the system, the linearizion of the system near a fixed point appears to be challenging oreven unfeasible. We approach this challenge by analyzing the interplay between the mutualistic in-teractions and the intra-guild competition, which separately lead to abundance gain and abundanceloss at equilibrium. When the mutualistic strength is equal to the competition strength, the speciesabundance on average follows (cid:104) s i (cid:105) = α i β i . Assuming the average abundance (cid:104) s Ak (cid:105) = (cid:104) s i (cid:105) for all theanimal species pollinating plant i , we linearize the nonlinear population dynamics at (cid:0) M P (cid:1) i = d Pi (cid:104) s i (cid:105) for each plant i , where d Pi = (cid:80) k K ik denoting the number of animals pollinating plant i . The fixedpoint for animal species M A is approximated similarly.Equation (6) provides an approximated solution for the abundances of general networks in the fulland 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) fordifferent parameters seeking solutions satisfying s P,Ai >
0. Equation 6 has the numerical advantagethat it allows one to scan the parameter space of ecological networks and delineate the feasible areamuch 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 variousmutualistic regimes. 6 .000.250.500.751.00 M u t ua li s m Feasible area ( a ) Analytical prediction
Competition M u t ua li s m Competition
Competition ( b ) ( c )( d ) ( e ) ( f ) h=0h=0.1 Full mean-field competition Soft mean-field competition Weighted competition
Figure 2: Feasible area patterns illustrated here for a real plant-pollinator network (MPL-16) from theWeb of Life platform [43], in the (left) full mean-field, (center) soft mean-field and (right) weightedscenarios. A point ( β , γ ) is colored in blue if all species survive with a positive abundance in thestationary 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 × β and γ , respectively. As it is seen, the analytical prediction delineates theboundaries of the feasible area with remarkable accuracy for the soft mean-field model, for both h = 0and h = 0 .
1. In the full mean-field model, reasonable precision is achieved for h = 0, while for h = 0 . Condition for the weighted competition
To analyze the impact of heterogeneity in the competitive strengths among intra-guild species, we de-rive conditions under which species coexist for the dynamical model with weighted competition. Theweighted competition is reduced to the homogeneous competitive strength (soft-mean field competi-tion) only when any pair of intra-guild species depend on and share exactly the same set of mutualisticpartners, corresponding to complete-like bipartite mutualistic networks.In the dynamic model for weighted competition [Eq. (15)], the inter-specific competition interac-tions are weighted by the relative importance of shared resources in a nonlinear form. To tackle the7onlinear inter-specific competition, we harness the microscopic perspective of intra-guild competitioninduced by a single mutualistic plant-pollinator interaction. When plant i is pollinated by an ani-mal k (i.e., K ik = 1), the inter-specific competition between plant i and the other plants j that arepollinated by animal k reads (cid:80) j ∈P ,j (cid:54) = i K ik K Tkj s Pj . Summing over all the pollinators k that pollinateplant i yields the total inter-specific competition (cid:80) k (cid:16)(cid:80) j ∈P ,j (cid:54) = i K ik K Tkj s Pj (cid:17) s Ak for that plant i . Armedwith the view of a single mutualistic interaction, the nonlinear competition term in Eq. (15) can berewritten as (cid:88) j ∈P ,j (cid:54) = i W Pij M Pi s Pj = (cid:80) k (cid:16)(cid:80) j ∈P ,j (cid:54) = i K ik K Tkj s Pj (cid:17) s Ak (cid:80) k K ik s Ak (7)Analytical estimates can be obtained by reigning in the weighted competition terms. This can beaccomplished by using the mediant inequality (see SI material):min k (cid:88) j ∈P ,j (cid:54) = i K ik K Tkj s Pj ≤ (cid:80) k (cid:16)(cid:80) j ∈P ,j (cid:54) = i K ik K Tkj s Pj (cid:17) s Ak (cid:80) k K ik s Ak ≤ max k (cid:88) j ∈P ,j (cid:54) = i K ik K Tkj s Pj (8)Equation (8) allows reducing the complexity of species abundance from two-guilds into a single-guild species abundance. After the complexity reduction, the weighted competition structure can beencoded in a competition matrix ˜ A P (resp. ˜ A A ) for plant (resp. animal) species, in analogy to thecompetition 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 thatplant i competes with plant j mediated by sharing the animal k , whose degree is, over all pollinators ofplant i , the closest to the local average of pollinated plants. Specifically, the element of the competitionmatrix ˜ A P can be written as˜ A Pij = K ik K Tkj (cid:54) = 0 for k ∈ arg min k ∈ A (cid:12)(cid:12)(cid:12) d Ak − (cid:106) (cid:80) s K is d As d Pi (cid:107)(cid:12)(cid:12)(cid:12) , d Pi ( d Ai ) is the number of animals (plants) with which plant (animal) i interacts, meaning thatthe weighted competition is approximated by an effective mutualistic partner k whose degree is theclosest 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 (cid:34) s P s A (cid:35) = βI + β (cid:34) ˜ A P
00 ˜ A A (cid:35) − (cid:18) γ ( hγ ( M P ) i ) (cid:19) K diag (cid:18) γ ( hγ ( M A ) i (cid:105) ) (cid:19) K T − (cid:32)(cid:34) α P α A (cid:35) + c (cid:33) (10)The analytical prediction of Eq. (10) for the weighted model is also well confirmed by numericalsimulations on real mutualistic networks [Fig. 2(c) and (f)]. Besides checking the validity of ourcalculations, Fig. 2 also allows us to highlight the marked differences in the dynamics yielded by thethree competition models. Notice, in particular, how strongly the shape of the feasible area changesfrom the full mean-field to the soft mean-field and then to the weighted competition case: the mere8 .0 0.2 0.4 0.6 0.80.10.20.30.40.5 F ea s i b l e a r ea Network size
Connectance Connectance Connectance
Full mean-field competition Soft mean-field competition Weighted competition a ) b ) c ) Interaction strength F ea s i b l e a r ea Interaction strength
Interaction strength d ) e ) f ) Figure 3: Feasible area as a function of different network metrics for full mean field [(a) and (d)], softmean-field [(b) and (e)], and weighted [(c) and (f)] competition scenarios. Each dot corresponds toa real mutualistic network from the Web of Life platform [43], with the size of dot proportional tothe network size. Parameters common in all panels: α i = 1 ∀ i , h = 0 . β = 5. Feasible area wascalculated considering the same ranges of β and γ shown in Fig. 2.inclusion of heterogeneity in the competition term in Eq. (13) shifts the region of occurrence of feasiblestates 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 formerallows a wider feasible region for low competition and high mutualism, while the latter is favored bymodest values of both competition and mutualism. These results clearly demonstrate how crucial thestructure of the intra-guild competition networks is for the species coexistence and the diversity ofecological communities. After investigating the impact of different choices of the structure of intra-guild competition networkson biodiversity, let us now study how the feasible area relates with the topological properties of themultilayer ecological networks. To do so, we consider a set of 50 real plant-pollinator networks retrievedfrom 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 averageinteraction strength and connectance. Following May [46], we define the average interaction strength asthe average of the off-diagonal elements of the Jacobian matrices of systems (13)-(15); and connectance9s defined as the density of non-zero values of the Jacobian matrices (see Supplemental Material). Asshown in Fig. 3, the feasible area correlates positively with connectance and interaction strength inthe full mean-field, soft mean-field and weighted competition scenarios. In other words, the higherthe number of interactions among species and the stronger their intensity, the more likely that theecological community exhibits feasible states. This result actually points back to the long debateddiversity-stability paradox initiated by May [46]. In his work, May proved that the probability of facingstable states converges almost certainly to zero for sufficiently large communities. In mathematicalterms, 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 ) = σ . Undersuch conditions, and setting the self-interaction terms as constants, J ii = − d ∀ i , May proved that thedynamical system d s /dt = Js is almost surely stable if √ N C < dσ , (11)in particular with d = 1, i.e. the condition for which the leading eigenvalue of J is negative. Conse-quently, the increases in size, connectivity and interaction strength favor the dynamical destabilizationof the system. This finding triggered the aforementioned paradox because it seems in contradictionwith 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 positivecorrelation 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 feasibilityarea, 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 acloser inspection suggests that some of these aspects are related. For example, incorporating the thirdelement of May’s relation and inspecting the size of the networks, we realize that the networks with ahigh connectance are also those that are the smallest networks in the database (Fig. 3, size of the dotsstands for number of species in the corresponding networks). What is more, by checking manually thenetworks with connectance values (cid:38) .
12, one notices that they correspond to almost fully connectedbipartite matrices and, hence, to almost fully connected Jacobian matrices as well. Thus, the networkswith 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 bythis 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 explorehow the size of the feasibility area is related to that criteria. Our goal here is to check whether thereis a clear relation between the feasible area and the expressions in Eqs. (13)-(15). However, it isnoteworthy that May’s condition [Eq. (11)] is related to the probability that the ecological system isstable for a particular set of parameters, whereas the feasible area results from a sum over differentparameter combinations. Therefore, to appropriately evaluate how the dependencies of the feasiblearea on network properties relate with Eq. (11), we define the following quantity: C May = (cid:104)(cid:104) J ii (cid:105)(cid:105) ( β,γ ) − (cid:104) σ ( J ij ) (cid:105) ( β,γ ) √ N C, (12)where (cid:104)·(cid:105) corresponds to an average over the Jacobian matrix’s elements, and (cid:104)·(cid:105) ( β,γ ) stands for theaverage over the parameters β and γ considered in Fig. 3. The first term in Eq. (25), (cid:104)(cid:104) J ii (cid:105)(cid:105) ( β,γ ) , is the10
10 15 20 250.10.20.30.40.5 F ea s i b l e a r ea
10 20 300.10.20.30.40.5 1 2 3 4 5 60.10.20.30.40.5
Network size
Full mean-field competition Soft mean-field competition Weighted competition a ) b ) c ) Figure 4: Comparison between feasible area and May’s stability condition [Eq. (25)] for several realplant-pollinator networks considering (a) full mean-field, (b) soft mean-field and (c) weighted intra-guild competition schemes. Parameters common in all panels: α i = 1 ∀ i , h = 0 . β = 5.average taken over the diagonal elements, since, contrarily to the random model considered by May, thediagonal elements of the Jacobian, J P , A ii , are not constant, but rather are heterogeneously distributedover the diagonal (see Supplemental Material); the term (cid:104) σ ( J ij ) (cid:105) ( β,γ ) is the average standard deviationof the off-diagonal values of J . In practical terms, variable C May defined in Eq. (25) quantifies howdistant a given network is from the critical point established by May’s stability condition, meaning howstable it is according to that criteria. For each real network considered, we numerically calculate theJacobian 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, thesize 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 thefact that May’s criteria was formulated for much more idealized systems (Fig. 4). The patterns wesee in Fig. 4 also agree with the scatter plots in Fig. 3; i.e., there is a noticeable correlation betweennetwork size and the value of the coefficient C May . In particular, we observe that large networks tendto exhibit low values for C May , while “less complex” networks (in terms of size and connectivity) havehigher values of C May .This raises the question of why the weighted competition model adheres with May’s criteria betterthan the full and soft mean-field scenarios. The answer lies in the expression of the Jacobian terms ofthe different models in Eqs. (13)-(15). For the weighted competition scenario, the off-diagonal Jacobianelements J ij are proportional to the terms 1 /s P,Ai and 1 / ( s P,Ai ) (see the SI material). This makes theelements J ij to narrowly peak around an average value, thus making the standard deviation σ ( J ij ) lesssignificant than the average of the diagonal terms J ii (see SI material) and thereby creating a positivecorrelation between feasible area and C May . Comparing the models of heterogeneous competitiveinteractions (i.e. soft mean-field and weighted in Fig. 4), one notices that the introduction of weightsin the competitive structure acts as a stabilizing factor in the dynamics, in the sense that networkstend to exhibit a wider feasible area in the weighted scenario. The interplay of interaction types amongspecies, a well-defined structure, and interacting strength shows a markedly driving force in stabilizingthe system that is highly predictable and consistent with May’s criterion. Interestingly, the feasible11rea of the weighted competition scenario correlates with May’s criterion as initially formulated by Mayfor random matrices, but no significant relation was observed with other, more recent stability criteriaformulated for random matrix models that incorporate features from predator-prey and mutualisticinteractions [44] (see SI material).
We investigated to which extent the incorporation of intra-guild competition alters the maintenanceof biodiversity in mutualistic systems. Compared to a scenario where all species from a guild homoge-neously compete with each other, as commonly assumed in the literature, heterogeneous competitionleads to markedly different patterns for the feasibility area of plant-pollinator networks. Withoutsufficient empirical data about intra-guild competition in plant-pollinator communities, deriving thestructure of the competitive links from the observed mutualistic interactions enables us to theoreticallyexplore how the structure of different intra-guild competition networks affects the structural stabilityof mutualistic ecological communities.Our results show that previously identified implications can be restricted to homogeneous com-petition and cannot be readily generalized to heterogeneous competition. Specifically, we found thatfeasibility patterns are dramatically modified when competitive interactions become heterogeneous.This finding suggests that a series of important conclusions regarding the dynamics of ecological net-works might have been overlooked given that theoretical models have been traditionally studied underthe assumption of homogeneous intra-guild competitive interactions. Therefore, getting informationon the structure of competitive networks in mutualistic systems is key to better understand what con-strains the assembly of mutualistic communities during the dynamical coevolutionary process, whichcould be an important driving force of coevolution [21].Finally, by investigating the feasible area in terms of global network properties, we found thatsmaller and less connected networks exhibited larger regions sustaining maximum biodiversity. Inter-estingly, this result agrees with the long-standing May’s stability-diversity paradox, which states thatcomplex systems are more prone to be destabilized as their size, connectance and mean interactionstrength increase. Indeed, we have verified that the more structurally stable networks turned out tobe the “less complex” ones according to May’s criterion [46]. Our results therefore show that the com-plexity introduced in the model by the weighted competition scenario yields a phenomenology whichis well predicted by a condition originally derived for random systems [Fig. 4(c)], whereas it is not thecase for the other competition scenarios. Importantly, the analysis on feasible solutions performed hereis not limited to the specific equations studied here, but can also be extended to ecological networkswith different types of interactions, such as facilitation [42], and predator-prey models [20].
Acknowledgments
X.W. acknowledges the project 62003156 supported by NSFC and “PCL Future Greater-Bay AreaNetwork Facilities for Large-scale Experiments and Applications (LZC0019)”. T.P. acknowledgesFAPESP (Grants No. 2016/23827-6 and No. 2018/15589-3). This research was carried out usingthe computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI)12unded by FAPESP (Grant No. 2013/07375-0). Y.M. acknowledges partial support from the Govern-ment of Aragon and FEDER funds, Spain through grant E36-20R (FENOL), by MINECO and FEDERfunds (FIS2017-87519-P), and by Intesa Sanpaolo Innovation Center. The funders had no role in studydesign, data collection, and analysis, decision to publish, or preparation of the manuscript.13 upplementary Material6 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 considerplant-pollinator systems composed of N A animal species, which interact with N P plant species. Thematrix encoding the mutualistic is the N P × N A bipartite matrix K , whose elements are K ij = 1if plant species i is pollinated by animal species j . The abundance of the i th plant is defined as s Pi ( s Ai , analogously for the animal species), and its time-depend dynamics is governed by the followingequations:1 s Pi ds Pi dt = α Pi − βs Pi − β N P (cid:88) j (cid:54) = i s Pj + γ M Pi hγ M Pi (Full mean-field competition), (13)1 s Pi ds Pi dt = α Pi − βs Pi − β N P (cid:88) j =1 A Pij s Pj + γ M Pi hγ M Pi (Soft mean-field competition), (14)1 s Pi ds Pi dt = α Pi − βs Pi − β N P (cid:88) j =1 (cid:32) W Pij M Pi (cid:33) s Pj + γ M Pi hγ M Pi (Weighted competition), (15)where α P,Ai are the intrinsic growth rate; parameters β and β stand for the intra- and inter-specificcompetition strength, respectively; and γ is the mutualism strength. Variable M Pi is given by M Pi = (cid:80) k ∈ A K ik s Ak , where k denotes the indexes belonging to the pollinators set. Matrix A Pij ofthe “Soft mean-field competition” scenario encodes the competitive connections within the plant-guild: A Pij = 1 if plants i and j share at least one pollinator, and 0 otherwise. For the “Weightedcompetition scenario”, the elements of matrix W P are defined as W Pij = (cid:80) k ∈ A K ik K jk s Ak , where s Ak is the abundance of the k th pollinator. The equations for ds Ak /dt are obtained by interchanging thelabels 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 whichregions of parameter space are actually accessible by a dynamical system one needs to consider boththe stability and the feasibility of the equilibrium solution. The feasibility alone is insufficient asnumerical 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 equilibriumsolution, 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 thepopulation dynamics can be written in a matrix form as dsdt = (cid:0) α T s − s T Ks (cid:1) , (16)We first linearize the system (16) around the feasible equilibrium s ∗ . Let us denote a small perturbationaround the equilibrium as u = s − s ∗ . Substituting s = u + s ∗ to Eq. (16) yields dudt = K eff u, (17)14here the entries of matrix K eff are given by K eff ij = αI − K ij s ∗ i , (18)where I is the identity matrix. By computing the eigenvalues of matrix K eff we can retrieve thestability of the equilibrium. The stability changes when one of the eigenvalues crosses zero. Suchpoints 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, weperform a stability analysis. That is we do not consider the equilibrium with one of the specieshaving zero (or negative) abundance. Taking matrix K of the full mean field as an example, feasibleequilibrium is given by the following expressions s ∗ = s ∗ = p q , s ∗ = s ∗ = p q , s ∗ = p q (19)with p = 25 + 5 β − β + 10 γ − β γ p = 25 − β + 5 γ − β γ − γ p = 25 − β + 10 γ + γ q = 125 + 50 β − β − β − γ + 2 β γ With the feasible equilibrium s ∗ and the competition matrix K , we construct the matrix K eff . Settingthe determinant of matrix K eff to zero, we find the expression for the curves with a vanishing eigen-value. In the parameter space of β and γ , the stability status changes at each time we across thosezero lines. From det (cid:0) K eff (cid:1) = 0, we obtain γ = ± (5 − β ) γ = ± (cid:112)
125 + 50 β − β − β √ − β which are the boundary curves for stability regions. Analogously, we describe the exact solution for thefeasible equilibrium and stability conditions for the soft mean-field and weighted competition modelsas shown in Figure 5.For the case of nonzero Holling term, h (cid:54) = 0, similar feasibility and stability analysis can beestablished by applying the linear approximation of the mutualistic terms. Figure 6 illustrates theanalytical results of feasibility and stability conditions for the full mean-field, soft mean-field andweighted competition scenarios, where panel (a) is for h = 0, and panel (b) is for h = 0 .
3. Figure 6suggests an altered feasible area modulated by inter-specific competition for the mutualistic networkshown in Figure 5.
To accurately predict the population dynamics of the weighted competition scenario, we derive theweighted competition matrix. For an observed mutualistic network, and a given initial condition where15 oft mean-field: Feasibility conditions s ∗ , s ∗ = p q , s ∗ , s ∗ = p q , s ∗ = p qp = 25 − β + 10 γ − β γ p = 25 − β + 5 γ − β γ − γ p = 25 − β − β + 10 γ − β γ + γ q = 125 + 25 β − β − β − γ + 4 β γ Stability conditions γ = ±√ (cid:112) − β γ = ± (cid:112)
125 + 25 β − β − β √ − β Weighted: Feasibility conditions s ∗ , s ∗ = p q , s ∗ , s ∗ = p q , s ∗ = p qp = (5 − β ) (5 + β + 2 γ ) p = 25 − β + 5 γ − β γ − γ p = 25 − β + 10 γ − β γ + γ q = (5 − β ) (cid:0)
25 + 10 β − β − γ (cid:1) Stability conditions γ = ±√ (cid:112) − β γ = ± β + 5 √ A Pij = K ik K Tkj (cid:54) = 0 for k ∈ arg min k ∈ A (cid:12)(cid:12)(cid:12) d Ak − (cid:106) (cid:80) s K is d As d Pi (cid:107)(cid:12)(cid:12)(cid:12) , d Pi ( d Ai ) is the number of animals (plants) with which plant (animal) i interacts. Term (cid:106) (cid:80) k K ik d Ak d Pi (cid:107) quantifies the local average of plant competitors mediated by sharing pollinators overnumber of pollinators d Pi of plant i . Recall that in the weighted scenario the competition intensitybetween two plants (pollinators) is directly proportional to the relative abundance of shared polli-nators (plants). Therefore, the more abundant the mutualistic partners, the fiercer the intra-guildcompetition 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 isremarkable, especially for the full and soft mean-field scenarios. Notice, in particular, how the detailedcontour 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 underestimatedby the theory, but nonetheless the general trends of the boundaries are reproduced very closely by thesolid curves. For the case h = 0 . γ , with exceptions at a high value of mutualism due to thesaturation effect. All in all, our analytical results work very well for a variety of real plant-pollinatornetworks. 16 Evaluation of analytical predictions on parameterization of pop-ulation dynamics
To investigate the performance of our analytical predictions of the feasible area, in this section wepresent 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 realplant-pollinator networks under different choices for the intraspecific competition β and the Hollingterm h . β Figures 9 and 10 show the numerical and analytical results when we increase or decrease the intensityof intraspecific competition β i for the observed mutualistic networks M-PL-048 and
M-PL-016 . Forthe population dynamics with soft mean-field competition, the analytical results accurately predictthe feasible area for both the increase and decrease of intraspecific competition β . For the weightedcompetition scenario, the linear approximation captures well the changing behavior of species coexis-tence. h (cid:29) Compared to the majority of the theories in ecology dealing with either weak mutualism ( h (cid:29)
1) orstrong mutualism h = 0, our theoretical framework fills the gap and covers the full range (0 < h < h = 0 and h (cid:29) h = 0, the dynamical model is reduced to Lotka–Volterra equationswith type II functional responses. Strong mutualism without saturation is often unstable leading toindefinite and unbounded growth of species which is argued to be biologically unrealistic due to theenvironmental constraints like carrying capacity [47]. However, strong mutualism at the same timeoverpowers the inter-specific competition resulting in a cut-off of the mutualism intensity γ .For a weak mutualism h (cid:29)
1, the mutualism term is saturated to a constant 1 /h , i.e.,lim h →∞ γ M Pi hγ M Pi ≈ h (21)When the mutualism is saturated, the dynamical model is reduced to a linear model, and the feasibleequilibrium is obtained by solving (cid:34) s P s A (cid:35) = (cid:32) βI + β (cid:34) A P A A (cid:35)(cid:33) − (cid:34) α P + h uα A + h u (cid:35) (22)where u is the all-one vector. For the weighted competition case, a linear approximation on thematrices ˜ A P and ˜ A A is applied. The mutualism intensity γ is disentangled from affecting the feasiblearea (which is confirmed by the numerical results in Fig. 11). Moreover, saturated mutualism, inturn, imposes limitations on the intensity β of competition. Because without the compensation ofmutualism benefit, the intra-specific competition might be overly strong and eventually make speciesgo extinct. 17 Structural stability and mutualistic network properties
In this section, we provide complementary results on the correlation between feasible area and networkarchitectures for dynamical models of full mean field, soft mean field and weighted competition. Wetest two more network properties, namely, the variance of the degree sequence and the ratio of inter-degree to intra-degree for both plants and animals. In addition, we test the maximum interspecificcompetition β before losing any species and its relation to network architectures when the mutualismis saturated.Here, we provide the correlation between feasible area and network architectures for all the con-sidered dynamical models. In particular, we consider the regimes of strong ( h = 0) and saturatedmutualism ( h (cid:54) = 0). For the strong mutualism regime (Fig. 12), global network properties, such asthe number of species and connectance, correlate strongly and negatively with the feasible area, withgoodness-of-fit of R = 0 .
83 and R = 0 .
81 for the soft mean-field competition. The maximum degreeshows a less strong correlation with the feasible area, and neither does nestedness show a clear de-pendence on the same quantity. The same trends are observed for other competition scenarios in thestrong 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 ischaracterized by dxdt = J x , indicating that a large system is less stable. In May’s assumption, eachelement of the Jacobian matrix J is equally likely to be positive or negative, having an absolutemagnitude 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 randommatrix [45], the system is stable if σ √ N C − ≤ σ is the standard deviation of the random variable from which the off-diagonal elements of theJacobian matrix J take value; N is the size of the network, and C denotes the network connectance [45]. 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 ] = σ and E [ m ii ] = µ . Then the largest eigenvalue is upperbounded by max | λ | ≤ σ √ N + O ( N / log N ) (24) where σ is the standard variation. Our goal now is to examine how the stability condition in Eq. 23 relates with the feasible area ofreal plant-pollinator networks. However, as discussed in the main text, Eq. (23) is associated withthe probability that the system is stable given a particular set of parameters; the feasible area, onthe other hand, is defined for a set of parameters. Therefore, in order to establish a relation betweenMay’s condition and the feasible area, we define C May = (cid:104)(cid:104) J ii (cid:105)(cid:105) ( β ,γ ) − (cid:104) σ ( J ij ) (cid:105) ( β ,γ ) √ N C, (25)18here (cid:104)·(cid:105) is an average over the Jacobian matrix’s elements, and (cid:104)·(cid:105) ( β ,γ ) corresponds to the averageover certain ranges of parameters β and γ . The first term in Eq. (25), (cid:104)(cid:104) J ii (cid:105)(cid:105) ( β ,γ ) , is the averagetaken over the diagonal elements, since, contrarily to the random model considered by May, thediagonal elements of the Jacobian, J P , A ii , are not constant, but rather are heterogeneously distributedover the diagonal (see Appendix); the term (cid:104) σ ( J ij ) (cid:105) ( β ,γ ) is the average standard deviation of theoff-diagonal values of J .The stability criterion by May is derived considering purely random interactions among the ele-ments in a complex system, i.e., the Jacobian matrix is a random matrix without any constraints on itselements. More recently, Allesina and Tang (AT) [44] generalized May’s stability analysis for randommatrix 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 C AT,Pred = (cid:104)(cid:104) J ii (cid:105)(cid:105) ( β,γ ) − (cid:104) σ ( J ij ) (cid:105) ( β,γ ) √ N C (cid:32) − (cid:10) E ( | J ij | ) (cid:11) ( β,γ ) (cid:104) σ ( J ij ) (cid:105) β,γ ) (cid:33) (26) C AT,Mutu = (cid:104)(cid:104) J ii (cid:105)(cid:105) ( β,γ ) − (cid:104) σ ( J ij ) (cid:105) ( β,γ ) √ N C (cid:32) (cid:10) E ( | J ij | ) (cid:11) ( β,γ ) (cid:104) σ ( J ij ) (cid:105) β,γ ) (cid:33) , (27)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 betweencompetitive and mutualistic relationships [44]. Figure 14 shows the observed correlation of May’s andAllesina and Tang’s stability with the feasible area for full mean-field, soft mean-field and weightedcompetition scenarios. Interestingly, we do not observe a significant dependence of the feasible areaon conditions (25)-(27) for full and soft mean-field scenario. A pattern does emerge in the weightedscenario, for which we have a positive correlation between feasible area and all conditions; however, itis interesting to note that the strongest correlation occurs for C May . The surprise about this resultresides in the fact that C May is the condition derived for the simplest random matrix model. At firstsight, 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: themost complex model (weighted competition scenario) adheres best with the condition derived fromMay’s stability criteria. As we argue in the main text, the explanation for the latter results lies inthe expression of the Jacobian elements (see Appendix). In the weighted competition scenario, theoff-diagonal Jacobian elements J ij depend on terms 1 /s P,Ai and 1 / ( s P,Ai ) . Since the abundance valuesare generally less than 1, the elements J ij end up being narrowly peaked around an average value, thusmaking 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 Jacobianmatrix for networks MPL-010 and MPL-048, respectively.In order to get further insights into the dynamics of the three competition models, let us denotethe right hand side of Eq. (15) as f ( s Pi ). To determine the stability of the equilibrium point, weanalyze the Jacobian matrix J which can be written as J = B − βI (28)19here B can be expressed in the form of block matrix as B = (cid:34) ( B ) N P × N P ( B ) N P × N A ( B ) N A × N P ( B ) N A × N A (cid:35) (29)Each element in B is calculated by ( B ) iu = ∂f (( s ∗ ) Pi ) ∂ ( s ∗ ) Pu , where s ∗ is the abundance at equilibrium.Each element in B is calculated by ( B ) iv = ∂f (( s ∗ ) Pi ) ∂ ( s ∗ ) Av . An analogous form can be obtained for B and B by replacing superscript P indicating plant species to superscript A representing pollinatorspecies and vice versa. We derive the expression of submatrices B and B for three cases of fullmean-filed, soft mean-field and weighted competition (see also Appendix A).(i) Case of full mean-field competition: the submatrix B = β J N P , where J N P is the all onematrix. Therefore, all elements in Jacobian submatrix B are linearly correlated. Each element in B is determined mainly by mutualistic interactions ( B ) iv = γ K iv ( hγ (cid:80) k K ik ( s ∗ ) Ak ) . All the interactingpollinators of plant i have the same value in Jacobian submatrix B and thus linearly correlated.(ii) Case of soft mean-field competition: the submatrix B has an element of 1 whenever there iscompetition which has the same value for all competed species. Each element in B is determinedmainly by mutualistic interactions ( B ) iv = γ K iv ( hγ (cid:80) k K ik ( s ∗ ) Ak ) .(iii) Case of weighted competition: the submatrix B has elements computed by( B ) iu = β (cid:80) k K ik K Tku ( s ∗ ) Ak (cid:80) k K ik ( s ∗ ) Ak (30)For each pollinator u of plant i , the value in Jacobian submatrix is varied, in contrast to the same valueof 1 in the case of full mean field and soft mean field. Each element in submatrix B is computed by( B ) iv = β (cid:80) j K iv K vj ( s ∗ ) Pj (cid:80) k (cid:54) = v K ik ( s ∗ ) Ak (cid:16)(cid:80) k K ik ( s ∗ ) Ak (cid:17) + γ K iv (cid:16) hγ (cid:80) k K ik ( s ∗ ) Ak (cid:17) (31)The second term in the right hand side of Eq. (31) is introduced due to mutualistic interactions andshows an analogous pattern to the full mean field case and soft mean field case. However, the first termuniquely appears in the weighted competition scenario. In addition, the value is varied for differentpollinators v , determined by the number of introduced plant competitions (cid:80) j K iv K vj ( s ∗ ) Pj mediatedby sharing a common pollinator v . Weighted competition reduces the correlation between elements ineach row of the Jacobian submatrix B and, therefore, shows a well agreement with May’s stabilitycriteria, 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 J = (cid:32) J P J P A J AP J A (cid:33) . (32)In the sequel we write the expression for each competition scenario.20 Competition M u t ua li s m Full mean-fieldSoft mean-fieldWeighted
Competition a ) b ) Figure 6: Feasible domain of the minimal model for full mean-field, soft mean-field and the weightedcompetition. Panel (a) shows the theoretical results for h = 0, and (b) shows the approximation resultfor h = 0 .
3. Other parameters are taken as α P = α A = 1, β P = β A = 5. A.1 Full mean-field competition model ∂ ˙ s Pi ∂s Pi ≡ J Pii = α Pi − βs Pi − β N P (cid:88) j (cid:54) = i s Pj + γ M Pi hγ M Pi , i ∈ P. (33) ∂ ˙ s Pi ∂s Pj ≡ J Pij = − s Pi β , i, j ∈ P. (34) ∂ ˙ s i ∂s Ak ≡ J P Aik = s Pi γ K ik (cid:0) hγ M Pi (cid:1) , i ∈ P, k ∈ A. (35) A.2 Soft mean-field competition model J Pii = α Pi − βs Pi − β N P (cid:88) j (cid:54) = i A Pij s Pj + γ M Pi hγ M Pi , i ∈ P. (36) J Pij = − s Pi β A Pij , i, j ∈ P. (37) J P Aik = s Pi γ K ik (cid:0) hγ M Pi (cid:1) , i ∈ P, k ∈ A. (38) A.3 Weighted competition model J Pii = α Pi − βs Pi − β (cid:80) j ∈ P,i (cid:54) = j W Pij s Pj M Pi + γ M Pi hγ M Pi , i ∈ P. (39) J Pij = − β s Pi M Pi W Pij , i, j ∈ P. (40)21igure 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,Ai = 1 ∀ i . 22igure 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,Ai = 1 ∀ i . 23 .00 0.25 0.50 0.75 1.000.00.20.40.60.81.0 M u t ua li s m MPL 16
Competition M u t ua li s m F ea s i b l e a r ea Analytical prediction
MPL 48
Competition a ) c ) b ) d ) Figure 9: Analytical prediction (black curve) for the soft mean-field, considering different values for theinter-specific competition strength β Ai = β Pi = β . The system in Eq. (14) was numerically integratedwith the Heun’s method, considering total simulation time T = 2000 and time step dt = 0 .
01. Thesimulation result (shaded area) is obtained with parameters α Pi = α Ai = 1 and h = 0 . J P Aik = − s Pi β K ik (cid:34) M Ak M Pi − s Pi M Pi − (cid:80) j ∈ P,i (cid:54) = j W Pij s Pj (cid:0) M Pi (cid:1) (cid:35) + s Pi γ K ik (cid:0) hγ M Pi (cid:1) , i ∈ P, k ∈ A. (41)24 .00 0.25 0.50 0.75 1.000.00.20.40.60.81.0 M u t ua li s m MPL 16
Competition M u t ua li s m Feasible area
MPL 48
Competition
Analytical prediction a ) c ) b ) d ) Figure 10: Analytical prediction (black curve) for the weighted competition scenario, consideringdifferent values for the inter-specific competition strength β Ai = β Pi = β . The systems in Eqs. (15)was numerically integrated with the Heun’s method, considering total simulation time T = 2000 andtime step dt = 0 .
01. The simulation result (shaded area) is obtained with parameters α Pi = α Ai = 1and h = 0 .
1. 25 .00 0.25 0.50 0.75 1.000.00.20.40.60.81.0 M u t ua li s m MPL-16
Analytical prediction
Competition M u t ua li s m MPL-48
MPL-16
Competition F ea s i b l e a r ea MPL-48
Soft mean-field competition Weighted competition a ) c ) b ) d ) Figure 11: Performance of analytical predictions (black lines) in the weak mutualism regime ( h (cid:29) α Pi = α Ai = 1, β = 5, and h = 10 . The systems in Eqs. (14) and (15) werenumerically integrated with the Heun’s method, considering total simulation time T = 2000 and timestep dt = 0 .
01. 26 Network size − 1 F ea s i b l e a r ea = xy − 0.74 + 0.7R = 0.66 Maximum degree − 1 F ea s i b l e a r ea = xy − 0.25 − 0.33R = 0.19 Connectance − 1 F ea s i b l e a r ea = xy − 0.45 + 0.45R = 0.76 − 1 Nestedness − 1 F ea s i b l e a r ea = 0.02 Network size − 1 = xy − 1.0 + 1.24R = 0.83 Maximum degree − 2 − 1 = xy − 0.21 − 0.39R = 0.09 Connectance − 1 = xy − 0.56 + 0.77R = 0.81 − 1 Nestedness − 2 − 1 = 0.13 Network size − 1 = xy − 0.89 + 1.03R = 0.61 Maximum degree − 2 − 1 = xy − 0.23 − 0.34R = 0.11 Connectance − 1 = xy − 0.52 + 0.64R = 0.7 − 1 Nestedness − 2 − 1 = 0.02 Full mean-field competition Soft mean-field competition Weighted competition
Figure 12: Feasible area and network architecture. Simulations are performed on 50 real-world mu-tualistic networks. Parameters to compute the feasible area are α = 1, β = 5, h = 0, β ∈ [0 ,
1] and γ ∈ [0 , Network size − 1 F ea s i b l e a r ea x Maximum degree − 1 F ea s i b l e a r ea = xy − 0.2 − 0.29R = 0.17 Connectance − 1 F ea s i b l e a r ea = xy − 0.39 + 0.4R = 0.78 − 1 Nestedness − 1 F ea s i b l e a r ea = 0.03 Network size − 1 Maximum degree − 1 = xy − 0.14 − 0.4R = 0.08 Connectance − 1 = xy − 0.37 + 0.39R = 0.82 − 1 Nestedness − 1 = 0.2 Network size − 1 Maximum degree − 1 = xy − 0.19 − 0.26R = 0.14 Connectance − 1 = xy − 0.36 + 0.36R = 0.68 − 1 Nestedness − 1 = 0.04 Full mean-field competition Soft mean-field competition Weighted competition R = 0.73=y − 0.66 + 0.67 R = 0.91= xy − 0.69 + 0.74 R = 0.66= xy − 0.64 + 0.68 Figure 13: Feasible area and network architecture. Simulations are performed on 50 real-world mu-tualistic networks. Parameters to compute the feasible area are α = 1, β = 5, h = 0 . β ∈ [0 ,
1] and γ ∈ [0 ,
10 15 20 25 F ea s i b l e a r ea = xy − 0.01 + 0.34R = 0.04
10 20 30 = xy − 0.01 + 0.42R = 0.23 = 0.59
10 15 20 25 F ea s i b l e a r ea = xy − 0.0 + 0.29R = 0.0
10 20 30 = xy − 0.01 + 0.4R = 0.09 −20 0 20 R = 0.27 F ea s i b l e a r ea = xy − 0.01 + 0.35R = 0.08
10 20 30 = xy − 0.01 + 0.41R = 0.29 = 0.46 a ) b ) c ) d ) e ) f ) g ) h ) i ) Full mean-field competition Soft mean-field competition Weighted competition
Figure 14: Correlation of feasible area with (a)-(c) May’s stability criterion [Eq. (25)], (d)-(f) stabilitycriterion by Allesina and Tang [44] for Predator-Prey models [Eq. (26)], and (g)-(i) criterion byAllesina and Tang for random matrices with competitive and mutualistic interactions [Eq. (27)]. Dotscorrespond to the real plant-pollinator networks with indexes between 01 and 50 retrieved from theWeb of Life database [43]. Size of the dots is proportional to the network size. The dynamics ofall networks was numerically integrated with the Heun’s method, considering total simulation time T = 1000 and a time step dt = 0 .
01. The feasible area was calculated over a β × γ grid with 100 × β , γ ∈ [0 , h = 0 . α i = 1 ∀ i , β = 5. Solid lines correspond tothe linear least-square regression, and R is the correlation coefficient.29 P r obab ili t y den s i t y −0.02 0.00 0.02 0.04050100150 0.0 0.1 0.2 0.3 0.40204060− 0.01 0.00 0.01 0.02 Jacobian element P r obab ili t y den s i t y − 0.02 0.00 0.02 0.04 Jacobian element
Jacobian element
Full mean-field competition Soft mean-field competition Weighted competition
Figure 15: Distribution of the elements of the Jacobian matrix for network MPL-10 for (upper row) h = 0 and (lower row) h = 0 .
1. Parameters α = 1, β = 5, β = 0 . γ = 0 . −0.01 0.00 0.01 0.0202505007501000 P r obab ili t y den s i t y −0.02 0.00 0.02 0.04 0.06050100150 0.0 0.1 0.20255075100−0.01 0.00 0.01 0.02 Jacobian element P r obab ili t y den s i t y −0.02 0.00 0.02 0.04 Jacobian element
Jacobian element
Full mean-field competition Soft mean-field competition Weighted competition
Figure 16: Distribution of the elements of the Jacobian matrix for network MPL-48 for (upper row) h = 0 and (lower row) h = 0 .
1. Parameters α = 1, β = 5, β = 0 . γ = 0 . eferences [1] .[2] S. Allesina and S. Tang. Stability criteria for complex ecosystems. Nature , 483(7388):205, 2012.[3] E. T. Aschehoug, R. Brooker, D. Z. Atwater, J. L. Maron, and R. M. Callaway. The mechanismsand consequences of interspecific competition among plants.
Annual Review of Ecology, Evolution,and Systematics , 47:263–281, 2016.[4] J. Bascompte and P. Jordano.
Mutualistic networks , volume 70. Princeton University Press, 2013.[5] J. Bascompte, P. Jordano, C. J. Meli´an, and J. M. Olesen. The nested assembly of plant–animalmutualistic networks.
Proceedings of the National Academy of Sciences , 100(16):9383–9387, 2003.[6] U. Bastolla, M. A. Fortuna, A. Pascual-Garcia, A. Ferrera, B. Luque, and J. Bascompte. Thearchitecture of mutualistic networks minimizes competition and increases biodiversity.
Nature ,458(7241):1018, 2009.[7] G. Blanco, F. Hiraldo, A. Rojas, F. V. D´enes, and J. L. Tella. Parrots as key multilinkers inecosystem structure and functioning.
Ecology and evolution , 5(18):4141–4160, 2015.[8] B. J. Brown, R. J. Mitchell, and S. A. Graham. Competition for pollination between an invasivespecies (purple loosestrife) and a native congener.
Ecology , 83(8):2328–2336, 2002.[9] J. H. Brown and M. A. Bowers. Community organization in hummingbirds: relationships betweenmorphology and ecology.
The Auk , 102(2):251–269, 1985.[10] D. R. Campbell and A. F. Motten. The mechanism of competition for pollination between twoforest herbs.
Ecology , 66(2):554–563, 1985.[11] L. Chittka and S. Sch¨urkens. Successful invasion of a floral market.
Nature , 411(6838):653, 2001.[12] J. H. Connell. On the prevalence and relative importance of interspecific competition: evidencefrom field experiments.
The American Naturalist , 122(5):661–696, 1983.[13] D. W. Davidson and S. R. Morton. Competition for dispersal in ant-dispersed plants.
Science ,213(4513):1259–1261, 1981.[14] P. C. de Ruiter, A.-M. Neutel, and J. C. Moore. Energetics, patterns of interaction strengths,and stability in real ecosystems.
Science , 269(5228):1257–1260, 1995.[15] P. Feinsinger. Effects of plant species on each other’s pollination: Is community structure influ-enced?
Trends in Ecology & Evolution , 2(5):123–126, 1987.[16] H. Fort. Quantitative predictions of pollinators’ abundances from qualitative data on their inter-actions with plants and evidences of emergent neutrality.
Oikos , 123(12):1469–1478, 2014.[17] E. C. Fricke, J. J. Tewksbury, E. M. Wandrag, and H. S. Rogers. Mutualistic strategies minimizecoextinction in plant–disperser networks.
Proc. R. Soc. B , 284(1854):20162302, 2017.3118] J. Ghazoul. Floral diversity and the facilitation of pollination.
Journal of ecology , 94(2):295–304,2006.[19] C. H. Graham, J. L. Parra, C. Rahbek, and J. A. McGuire. Phylogenetic structure in tropicalhummingbird communities.
Proceedings of the National Academy of Sciences , 106(Supplement2):19673–19678, 2009.[20] T. Gross, L. Rudolf, S. A. Levin, and U. Dieckmann. Generalized models reveal stabilizing factorsin food webs.
Science , 325(5941):747–750, 2009.[21] P. R. Guimar˜aes Jr, M. M. Pires, P. Jordano, J. Bascompte, and J. N. Thompson. Indirect effectsdrive coevolution in mutualistic networks.
Nature , 550(7677):511, 2017.[22] T. C. Ings, J. M. Montoya, J. Bascompte, N. Bl¨uthgen, L. Brown, C. F. Dormann, F. Edwards,D. Figueroa, U. Jacob, J. I. Jones, et al. Ecological networks–beyond food webs.
Journal ofAnimal Ecology , 78(1):253–269, 2009.[23] L. K. Johnson and S. P. Hubbell. Aggression and competition among stingless bees: field studies.
Ecology , 55(1):120–127, 1974.[24] S. K´efi, E. L. Berlow, E. A. Wieters, S. A. Navarrete, O. L. Petchey, S. A. Wood, A. Boit,L. N. Joppa, K. D. Lafferty, R. J. Williams, et al. More than a meal. . . integrating non-feedinginteractions into food webs.
Ecology letters , 15(4):291–300, 2012.[25] R. M. May. Will a large complex system be stable?
Nature , 238(5364):413, 1972.[26] K. S. McCann. The diversity–stability debate.
Nature , 405(6783):228–233, 2000.[27] R. J. Mitchell, R. J. Flanagan, B. J. Brown, N. M. Waser, and J. D. Karron. New frontiers incompetition for pollination.
Annals of botany , 103(9):1403–1413, 2009.[28] A. Montesinos-Navarro, F. Hiraldo, J. L. Tella, and G. Blanco. Network structure embracingmutualism–antagonism continuums increases community robustness.
Nature ecology & evolution ,1(11):1661, 2017.[29] A.-M. Neutel, J. A. Heesterbeek, and P. C. de Ruiter. Stability in real food webs: weak links inlong loops.
Science , 296(5570):1120–1123, 2002.[30] T. M. Palmer, M. L. Stanton, and T. P. Young. Competition and coexistence: exploring mech-anisms that restrict and maintain diversity within mutualist guilds. the american naturalist ,162(S4):S63–S79, 2003.[31] C. Robertson. The philosophy of flower seasons, and the phaenological relations of the ento-mophilous flora and the anthophilous insect fauna.
The American Naturalist , 29(338):97–117,1895.[32] R. P. Rohr, S. Saavedra, and J. Bascompte. On the structural stability of mutualistic systems.
Science , 345(6195):1253497, 2014. 3233] D. W. Roubik. Foraging behavior of competing africanized honeybees and stingless bees.
Ecology ,61(4):836–845, 1980.[34] D. W. Roubik, J. E. Moreno, C. Vergara, and D. Wittmann. Sporadic food competition withthe african honey bee: projected impact on neotropical social bees.
Journal of Tropical Ecology ,pages 97–111, 1986.[35] M. L. Stanton. Interacting guilds: moving beyond the pairwise perspective on mutualisms.
TheAmerican Naturalist , 162(S4):S10–S23, 2003.[36] D. B. Stouffer and J. Bascompte. Compartmentalization increases food-web persistence.
Proceed-ings of the National Academy of Sciences , 108(9):3648–3652, 2011.[37] E. Th´ebault and C. Fontaine. Stability of ecological communities and the architecture of mutu-alistic and trophic networks.
Science , 329(5993):853–856, 2010.[38] C. R. Townsend, M. Begon, J. L. Harper, et al.
Essentials of ecology.
Number Ed. 2. BlackwellScience, 2003.[39] C. Wagg, J. Jansa, M. Stadler, B. Schmid, and M. G. A. Van Der Heijden. Mycorrhizal fungalidentity and diversity relaxes plant–plant competition.
Ecology , 92(6):1303–1313, 2011.[40] P. Wilmer.
Pollination and floral ecology . Princeton University Press, 2011.[41] W. Wilms and B. Wiechers. Floral resource partitioning between native melipona bees and theintroduced africanized honey bee in the brazilian atlantic rain forest.
Apidologie , 28(6):339–355,1997.[42] F. Z´el´e, S. Magalh˜aes, S. K´efi, and A. B. Duncan. Ecology and evolution of facilitation amongsymbionts.
Nature communications , 9(1):1–12, 2018.[43] .[44] S. Allesina and S. Tang. Stability criteria for complex ecosystems.
Nature , 483(7388):205, 2012.[45] Z. F¨uredi and J. Koml´os. The eigenvalues of random symmetric matrices.
Combinatorica ,1(3):233–241, 1981.[46] R. M. May. Will a large complex system be stable?
Nature , 238(5364):413, 1972.[47] D. H. Wright. A simple, stable model of mutualism incorporating handling time.