The joint influence of competition and mutualism on the biodiversity of mutualistic ecosystems

In the past years, there have been many advances –but also many debates– around mutualistic communities, whose structural features appear to facilitate mutually beneficial interactions and increase biodiversity, under some given population dynamics. However, most approaches neglect the structure of inter-species competition by adopting a mean-field perspective that does not deal with competitive interactions properly. Here, we build up a multilayer network that naturally accounts for mutualism and competition and show, through a dynamical population model and numerical simulations, that there is an intricate relation between competition and mutualism. Specifically, the multilayer structure is coupled to a dynamical model in which the intra-guild competitive terms are weighted by the abundance of shared mutualistic relations. We find that mutualism does not have the same consequences on the evolution of specialist and generalist species, and that there is a non-trivial profile of biodiversity in the parameter space of competition and mutualism. Our findings emphasize how the simultaneous consideration of positive and negative interactions derived from the real networks is key to understand the delicate trade-off between topology and biodiversity in ecosystems and call for the need to incorporate more realistic interaction patterns when modeling the structural and dynamical stability of mutualistic systems.


Multilayer Networks
A graph (i.e., a single-layer network) is a tuple G = (V, E), where V stands for the set of nodes and E ⊂ V × V is the set of edges that connect pairs of nodes. The edges of the graph induce a binary relation on V that is called the adjacency relation of G. This relation can be represented through the adjacency matrix A, where A i,j indicates the number of links from node i to node j.
As a particular case, a bipartite graph is a graph whose nodes can be divided into two disjoint sets V 1 and V 2 (V = V 1 ∪ V 2 , V 1 ∩ V 2 = ∅) such that every edge connects a node in V 1 to one in V 2 . The adjacency matrix of a bipartite graph has the form of a block matrix: where O represents the n × n null matrix, O ij = 0.
In order to build a more complete network containing the interactions among species of different guilds as well as intra-guild interactions, we can further exploit the information contained in the bipartite network of mutualistic interactions and consider the projections of the mutualistic interactions onto the guilds of plants (P ) and animals (A). This set-up gives us a multilayer network where the guilds constitute the layers, the mutualistic interactions provide the interlayer links and their projections account for the intra-layer links, which in turn represent the number of pollinators (resp., plants) shared by the corresponding plants (resp., pollinators). Accordingly, the resulting multilayer network consists of two layers (V P , V A ), the set of inter-layer links (C ∈ V P × V A ), and two sets of intra-layer For the sake of clarity, here we adopt the notation commonly used in the biological literature. The mutualistic plant-pollinator interactions (i.e., the inter-layer connections) are given by a bipartite N P × N A matrix, K, with K ik = 1 if animal species k pollinates the plant species i, and K ik = 0 otherwise. The projection of this matrix onto the sets of plants is given by: where V P ij represents the number of pollinators shared by plant species i and j. In our dynamical equations, we consider the resources that are shared by species of the same guild. According to equation 2, the abundance of the pollinators shared by plant species i, j is given by: Equivalently, the respective formulas for pollinators are obtained by interchanging the P and A superscripts and vice versa, that is, Tables 1 and 2 show the main features (number of species, plants and animals), the location of the biotope and the reference of the ecological networks used to numerically solve the differential equations. Table 1 refers to plant-pollinator networks and table 2

Role of the heterogeneity in the competition parameters
Here we explore a simplified version of the model given by eq.(1) and (2) of the main text in order to investigate if the obtained results could be just a consequence of heterogeneity in the competition intensity. To do so, we consider a situation where the competition term between species i and j of the same layer reads:C ij = β ijṼij s j ,Ṽ ij is the binary projected matrix corresponding to each guild, as explained in the main text, and β ij are taken from an uniform distribution (0.05β 0 , 1.95β 0 ). It is worth noticing that this introduces a quenched disorder, at a difference with our model where the competition intensity is a dynamical variable. We show below that introducing heterogeneity in the competition term is not enough to reproduce the results observed in the right panel of Figure 2 of the main text. Figures 1, 2 and 3 represent the species persistence as a function of the interspecies competition and mutualistic parameters (resp., β 0 and γ 0 ). The rows correspond to different real mutualistic networks, plant-pollinator networks in Fig. 1-2 and seed-dispersal networks in Fig. 3. In all figures, left panels show the results of integrating the dynamical equations when the competitive term is homogeneous β ij = β 0 , for all interacting species of a given layer (some of which are also shown in the main text). In the central panels, which correspond to the heterogeneous competition parameter, the same qualitative behavior holds but the final biodiversities are lower than in the homogeneous case. Finally, the right panels correspond to the model defined by equations equations (1) and (2) of the main text, and therefore they represent the case in which the competition terms have been weighted by the shared abundance of pollinators (resp., plants). Again, some panels have been introduced in the main text. A direct comparison of the central and right panels clearly shows that the obtained results are not just a mere consequence of heterogeneity in the competition parameter, but they stem from the dynamical character of the competition term associated to the network structure.

Robustness of the model
Real mutualistic networks exhibit a wide range of size, connectance and nestedness. The connectance of a network is defined as the fraction of actual mutualistic links with respect to the total possible links in the network: N P × N A . Among the different metrics to measure the nestedness, here we use the NODF metric [15], which is bounded to the [0, 1] interval, where higher values indicate higher nestedness. In order to test the robustness of the proposed dynamics under the election of the network, we chose for this study real networks covering a wide range of size, connectance and nestedness. As shown in tables 1 and 2, the size of the selected networks varies from 27 to 719 species. Figure 4 represents the connectance as a function of the nestedness for all the plant-pollinator networks on the web-oflife.es database, highlighting the correlation between connectance and nestedness. The networks studied here are shown by red squares.
We have also checked that the results are not dependent on the value of the Holling constant that controls the saturation of the mutualist term. To do so we have integrated the dynamical equations for different values of h. Figure 5 represents the species persistence as a function of the inter-species competition and mutualistic terms (resp., β 0 and γ 0 ). The different rows correspond to the different mutualistic networks shown in the main text, while each column corresponds to a different value of the Holling constant h = 0.1, 0.5, 0.9. As it can be seen, the qualitative behavior remains unchanged under the variation of h, that is, regardless the choice of h, in the high-competition regime an increase in mutualism involves a decrease in persistence.

Additional intra-species competition term
The model described by equations (1) and (2) of the main text implies that the competition between two plants (respectively, pollinators) is proportional to the abundance of the shared pollinators (respectively, plants). Nevertheless, there are situations in which competition is not driven exclusively by pollination, as in the case of competition for abiotic resources like water, light or nutrients. In order to take into account this fact, here we explore the addition of an extra inter-species competition term to the equation (1) of the main text. As this extra competition term is not related to the projected networks, each species competes with all the others in the layer, with a uniform parameter β c .
According to this alternative formulation, the equation that describes the evolution of the relative abundances of a given plant i now takes the form: The rest of the terms remain invariant with respect to equation (1) of the main text: the first term of equation 6 represents the intrinsic growth of species i, the second term refers to the intra-specific competition term (saturation), the fourth term accounts for the inter-specific competition for pollinators and, finally, the fifth term gives the contribution of mutualism. Figure 6 shows the species persistence as a function of the inter-species competition and mutualistic terms (resp., β 0 and γ 0 ). The rows correspond to different real mutualistic networks as indicated. Left panels show the results of the approach described in this section, that is, when an extra inter-species competitive term not driven by pollinators is added to the equation for the evolution of plant species. In these simulations, the value of β P c has been chosen equal to the one used for the pollinators-driven competition β P 0 . Right panels of Fig. 6 correspond to the model defined by equations (1) and (2) of the main text, according to which the inter-specific competition is driven exclusively by the species of the other guild.
As can be seen by comparing the results of both approaches, the introduction of the additional term does not have a qualitative effect on the system's behavior, what means that the results are robust with respect to this modification. Similar results are obtained when the additional term described in equation 6 is introduced in the evolution of animal species or in both guilds.
5 Approximated analysis of the behaviour of the competition term.
Here we discuss an approximated version of the heuristic argument provided in the Methods section of the main text which provides a simple interpretation of the reasons why the structured competition term affects generalists and specialists in different ways. Let us consider the evolution of the abundance of a plant species i, s P i . The factor W P ij of equation 3 involves the abundance of each animal species of the other guild. These abundances are all different because for the animal guild, generalists and specialists also grow faster and slower with γ 0 , respectively. We will then consider an average relative abundance for animals s A . Then the competing term, for plants, reads: where V P ij is the projected matrix on the plant space. As V P ij ≤ min(k P i , k P j ) we can consider two situations: • i is a generalist so k P i ≥ k P j ∀j then, the maximum possible contribution to the competition term is : • i is a specialist, so k P i < k P j then the maximum contribution to the competition term is : C M AX spe ≈ β P 0 j∈P,j =i s P j which shows that the specialists are, on average, more affected by competition than generalists when mutualism increases. The same reasoning holds symmetrically for the animal's guild.  Figure 3: Biodiversity in seed-dispersal networks under different treatments of the competition term. Number N of species in the steady state as a function of the inter-species competition and mutualistic parameters (resp., β 0 and γ 0 ), for different seed-dispersal networks (rows). Left panels: results for the homogeneous competition model. All interacting pairs of species in the layer have the same interacting constant β ij = β 0 ∀(i, j). Central panels: results for the heterogeneous competition model. Interacting constants are drawn from a uniform distribution β ij ∈ [0.05β 0 , 1.95β 0 ]. Finally, the right panels stand for the model given by eq.
(1) and (2) of the main text. The maximum diversity corresponds to N = 209 (resp., 318) species for network M-SD-019 (resp., M-SD-022). Other values are: Dynamical competition coupling Figure 6: Additional term in the intra-species competition. Number N of species in the steady state as a function of the inter-species competition and mutualistic terms (resp., β 0 and γ 0 ), for different mutualistic networks (rows). Left panels show the results when a mean-field term for the intra-species competition is added to the model given by eqs. (1) and (2) of the main text. Right panels stand for the model here presented. The maximum diversity corresponds to N = 266 (resp., 719) species for network M-PL-044 (resp., M-SD-044). Other values are: β P i = β A k = 5.0, α P i ∈ (0.9, 1.1), α A k ∈ (0.9, 1.1), β P c = β P 0 and h P = h A = 0.1.