Balance of interactions determines
optimal survival in multi-species communities

Anshul Choudhary and Sudeshna Sinha

Department of Physical Sciences, Indian Institute of Science Education and Research (IISER) Mohali, Knowledge City, SAS Nagar, Sector 81, Manauli PO 140 306, Punjab, India.

* E-mail of corresponding author: sudeshna@iisermohali.ac.in

## Abstract

We consider a multi-species community modelled as a complex network of populations, where the links are given by a random asymmetric connectivity matrix , with fraction of zero entries, where reflects the over-all connectivity of the system. The non-zero elements of are drawn from a gaussian distribution with mean and standard deviation . The signs of the elements reflect the nature of density-dependent interactions, such as predatory-prey, mutualism or competition, and their magnitudes reflect the strength of the interaction. In this study we try to uncover the broad features of the interspecies interactions that determine the global robustness of this network, as indicated by the average number of active nodes (i.e. non-extinct species) in the network, and the total population, reflecting the biomass yield. We find that the network transitions from a completely extinct system to one where all nodes are active, as the mean interaction strength goes from negative to positive, with the transition getting sharper for increasing and decreasing . We also find that the total population, displays distinct non-monotonic scaling behaviour with respect to the product , implying that survival is dependent not merely on the number of links, but rather on the combination of the sparseness of the connectivity matrix and the net interaction strength. Interestingly, in an intermediate window of positive , the total population is maximal, indicating that too little or too much positive interactions is detrimental to survival. Rather, the total population levels are optimal when the network has intermediate net positive connection strengths. At the local level we observe marked qualitative changes in dynamical patterns, ranging from anti-phase clusters of period cycles and chaotic bands, to fixed points, under the variation of mean of the interaction strengths. We also study the correlation between synchronization and survival, and find that synchronization does not necessarily lead to extinction. Lastly, we propose an effective low dimensional map to capture the behavior of the entire network, and this provides a broad understanding of the interplay of the local dynamical patterns and the global robustness trends in the network.

## Introduction

Complex networks provides a common framework to address a vast range of phenomena in large interactive systems [1]. The use of network theory in studying the stability and dynamics of model ecosystems started with the landmark paper of Robert May[2] and the success and effectiveness of such inquiry can be gauged from the fact that even today most studies in theoretical ecology heavily rely on the network framework [3]. Our current understanding of the stability of an ecological network hinges around two key aspects: interaction topology and nature of interactions. Now, there are wide ranging situations where one either does not have sufficient information on the exact underlying topology, or one finds that the web of interactions essentially appears to be a random network. In such cases, the interactions are modelled by random connectivity matrices, and the broad nature of interactions is our only guiding principle in analyzing the dynamics and survival properties of the complex system.

In this study, we are going to explore the effect of the balance of different kinds of interactions in a multi-scpecies community on the collective dynamical behaviour of the network. Our focus will be on the global robustness of the system, as exemplified by the total population of all species [4]. It is evident that the total population relects the global state of the network effectively, while being sensitive to the underlying dynamics at the local species level as well. So here we will explore how diversity in interactions influence the emergent dynamics, and the relation of these dynamical patterns to survival of populations, lending yet another perspective to the stability-diversity debate.

Recently, Mougi & Kondoh[5] studied the interesting effects of diversity in interaction types on the stability of an ecological community and they found that diversity is a key element in determining stability and biodiversity. However their results are based on linear stability analysis for small perturbations about local equilibria, and they do not give the relationship between survival and the emergent dynamical patterns. In this context, our study provides a complementary exploration of the global survival features of such systems [6] and also relates it to dynamical behavior of the constituent populations.

## Model

The model we consider here is inspired by the earlier theoretical studies conducted by Robert May [2, 7]. However, we would like to mention here that unlike most studies regarding stability [8, 5] which assumes the existence of time-independent population densities when system reaches steady state, we consider the more general condition where the attracting state can have complex temporal behavior, rather than a fixed point solution [6]. The principal motivation for this approach to the question of stability is its wider relevance and broader applicability. Further, rather than local stability about an equilibrium, we will focus on a different set of global quantifiers of robustness and survival in the complex network.

Specifically, in this work we consider a prototypical map, the Ricker (Exponential) Map, modelling population growth of species with non-overlapping generations:

(1) |

with interpreted as an intrinsic growth rate.

We then consider the evolution of such interacting populations given as:

(2) |

where , and the connectivity or community matrix, represents how species are mutually interacting. Further we consider that if , where is the minimum population density below which population cannot sustain on their own and therefore becomes extinct, namely the Allee Effect [9]. So we have a system (Eqn. 2) with well mixed populations at the nodes that display chaotic dynamics, and are extinction prone due to the population threshold when uncoupled. So clearly, a persistent community can only be sustained through suitable interactions among the species.

The connectivity matrix is a random asymmetric matrix, with fraction of zero entries, where reflects the over-all connectivity of the system. The non-zero elements are drawn from a gaussian distribution, , with mean and standard deviation .

The signs of the elements of the connectivity matrix reflects the nature of density-dependent interactions. In general, neutral interactions are reflected by zero matrix elements. Interactions that reduce the population at a node, for instance through parasitism, grazing, and predation, will be reflected by a negative sign, while interactions that benefit a species, for instance by provding refuge from physical stress, predation or competition, will bear a positive sign.

So then, when we have mutualism or symbiosm, both and are positive. When we have competition or antagonistic interactions, both and are negative. When the effect of one species on the other is positive, but neutral the other way round, we have commensalism, reflected by and . Similarly, amensalism is reflected by and . General predator-prey interactions are captured by and having different signs.

In general, positive interactions or facilitative interactions between species that benefit the growth of the species, give rise to more positive mean . Namely, high positive implies the dominance of mutualism in the ecosystem where as would imply balance of different kinds of interactions in the system. The standard deviations, on the other hand controls the degree of variability in the strength of these interactions. Finally, connectedness, is another important factor that tell us how many species are interacting with each other. Namely, it reflects the fraction of neutral interactions in the network. The main aim of this study is to understand the relation between broad features of the interaction matrix and the collective dynamics of the system, and then go on to link this to the local and global survival in the system.

## Methods

At the outset, we present our tools and describe the measures for analyzing the survival properties of the system. To gauge the robustness of the system, we first calculate the number of active nodes, namely the number of non-extinct species with non-zero population, after transience. This quantity is then averaged over a period of time and further averaged over different initial conditions. We denote this averaged number of non-extinct nodes by , and it is defined as:

(3) | ||||

(4) |

The next important measure of global survival is the total population , of the system, reflecting the biomass yield in a multi-species community. This quantity, averaged over a period of time and over different initial realizations, is denoted by , and mathematically expressed as:

(5) |

Synchronization Order Parameter :

In order to probe collective patterns in the network, we studied the level of synchronization that emerges in the system. To quantify the degree of synchronization we have employed two different order parameters.

(i) : Here we measure synchronization error as the mean square deviation of the local state of the nodes, averaged over time (after transience) and over different realizations [23, 24, 25], mathematically expressed as :

(6) |

When it goes to zero, this measure reflects complete synchronization in the system.

(ii) : This is a phase order parameter that reflects the degree of variation in the phases of the local dynamics at the nodes. Specifically, it is a measure of the fraction of nodes in the largest phase cluster, averaged over time and over different network realizations . When , it implies that the entire system is phase synchronized (though not necessarily in complete synchronization).

For the specific case of the local dynamics being a -cycle or being in a period chaotic band (which is observed in this system over a large parameter range) then is the supremum of the quantity

where is if lies in the specified band, and otherwise.

So these measures provide complementary information about the synchrony in phase and amplitude of the dynamics of the local constituents of the network.

System Parameters:

In this work parameter in Eq. 1, namely the local map is in the chaotic regime, and the threshold value . All results reported here are robust with respect to small variations around these values. The survival and synchronization measures were calculated by averaging over random initial conditions, i.e. in the equations above. The system sizes ranged from , with connectedness and standard deviation in the connectivity matrix . Further, to explore the effect of the mean interaction strength on the dynamics, which is a focus of our work here, we investigated the range: . Note that several earlier studies have been confined to the balanced situation, where . In the sections below, we present the principal observations from our extensive simulations over this wide-ranging window of parameters.

## Results and Analysis

### Survival in the Network

We first calculate the average number of active nodes (namely the average number of non-extinct species) , as a function of the mean interaction strength of the connectivity matrix . As evident from the results displayed in Fig. 1, the average number of active nodes in the network rises sharply as a function of mean interaction strength around . When the mean interaction strength is quite negative, the number of active nodes goes to zero, i.e. the entire system is driven to extinction. For positive all nodes in the network are non-zero, i.e. no species goes extinct. The connectedness and the variability of interaction strengths then does not affect the number of active nodes in the network when the network is far from balanced, namely considerably positive mean interaction strengths yield , while considerably negative interactions results in , irrespective of and . The transition from complete extinction to a completely active network is sharper for systems with low variability in interaction strengths (i.e. low ), and for systems with higher connectedness (i.e. high ).

To gain further quantitative understanding of the nature of this transition, we explore the scaling behaviour near the transition, and discover that the average number of active sites scales with respect to as:

(7) |

Further the number of active sites scales with respect to as:

(8) |

Here, and are appropriate critical exponents for scaling functions and respectively. A good data collapse, shown in Fig. 1 (insets), is obtained for , with and . These scaling relations suggest that a transition from complete extinction to a fully active ecosystem occurs around , namely around the state of balanced interaction strengths or completely neutral interactions [10]. We also performed finite size scaling with respect to system size , and found the simple scaling: , implying that the active fraction is independent of .

The next important measure of global survival is the average total population of the system, reflecting the biomass yield in a multi-species community. The variation of , as a function of mean of the interaction strengths and connectedness of the interaction matrix , is shown in Figs. 2-3. It is clear that for , i.e, when there are no interactions, the local extinctions accumulate, eventually leading to mass extinction.

When interactions are present, different global scenarios emerge with respect to varying mean and connectedness of matrix . For fixed , the total population increases sharply with increasing , namely with increasing net positive interactions, around . So we find that for networks close to the balanced situation, we have enhanced population densities indicating greater survival, for increasing net positive interactions.

Further, there exists an interval of mean , around , where the average total population always increases with the increase in the number of interactions among species, namely increasing . This would imply that connectivity always enhances survival of the system here. However, when mean is smaller, or larger, than the above interval, one finds that at intermediate values of connectedness, the population is the largest(cf. Fig.2). Namely, a mix of neutral interactions along-side other interactions is most conducive to enhanced population yeild.

There is a critical negative mean, (where is a function of ) for which the local species experience severe loss of population leading to global extinction. There is also a critical positive mean, , where , such that for mean the nodes experience unbounded and explosive growth, destabilizing the whole network. We consider in our study.

Interestingly, we uncovered a scaling pattern between the total population and characteristics of the connectivity matrix. The data collapse of the population onto a single non-monotonic curve in Fig. 3B reveals a scaling relation between total population and the product of the mean interaction strength, and connectedness, of the network. This implies that the most relevant quantity in understanding the global behaviour of the network is , rather than and alone. So clearly, survival is dependent not merely on the number of links, but on the combination of the sparseness of the connectivity matrix and the net interaction strength. For instance, fewer interactions (i.e. low ) tends to decrease the population, but this effect may be compensated by more positive interactions, i.e. higher . More importantly, the existence of an intermediate window of positive where the total population is maximum indicates that too little or too much positive interactions is detrimental to survival. Infact survival is optimal when the network has intermediate net positive connection strengths. So counter intuitively, if positive interactions such as mutualistic or symbiotic links dominate other kinds of interactions too much, its effect ceases to be beneficial, causing the total population to reduce.

### Local Dynamics

Now we attempt to correlate these global features to local species-level dynamics. Namely, we attempt to correlate the survival and global stability of the ecological network to dynamical patterns emerging in the network as a result of interactions.

From the bifurcation diagrams displayed in Fig. 4, one can clearly discern the presence of coherent collective dynamics in the system. This coherence breaks down as one approches , as evident in Fig. 4, with the network displaying unsynchronized spatio-temporal chaos.

In order to gauge the degree of synchronization among the nodes quantitatively we calculate the synchronization order parameter . Our attempt now will be to find the correlation between synchronization and survival. This is an important question, as synchronization has often been seen as increasing risks of extinction. Fig. 5 exhibits this synchronization error, along side the number of active patches (i.e. nodes with non-zero population), the total population and collective dynamics of the whole network. It is clear that synchronization does not necessarily lead to extinction. In fact for positive mean interactive strengths, even when the entire system is completely synchronized, all nodes are active. The rationale for the above observation is that synchronization is a threat only when the synchronized dynamics covers a large range of population densities, such as in synchronized chaos, which typically is ergodic over state space. Here on the other hand, the synchronized dynamics is confined to the “safe zone” and the attractor trajectory does not enter the extinction region [11]. So the synchronized patches survive.

We further investigate the nature of the time series of the local patches to discern cluster formation, and the phase relation between the clusters. We find that when the mean of the interaction strengths has a low positive value, the populations are attracted to a period cycle, and the sytem divides into two anti-phase clusters (cf. Fig. 6). Namely, alternately in time, one set of nodes in the network have low population densities, while the other set has high population densities. This behaviour is reminiscent of the field experiment conducted by Scheffer et al [12] which showed the existence of self-perpetuating stable states alternating between blue-green alage and green algae.

We also studied the phase clusters emerging in the system, by calculating a phase order parameter , which gives the fraction of species whose dynamics are in-phase in the network. This quantity is in a large range of positive mean interaction strengths (), indicating that here the network always splits into two approximately equal clusters. In each cluster the nodes are in-phase with respect to each one another, and anti-phase with respect to nodes in the other cluster. However note that the degree of complete synchronization, which depends on both phase and amplitude, will be dependent on (as evident from Fig. 5). So changing the mean interaction strength changes the nature of the dynamics without destroying this phase relationship and two phase synchronized clusters of varying amplitudes emerge in a large range of .

### Effective map for nodal dynamics

To gain further understanding of the dynamical patterns, we construct an effective map to mimick the essential dynamics of the nodes. Our approach is to split the interactive part in Eqn. 2 into an average part and a term capturing the fluctuations. Here the mean interaction strength, which is the dominant contribution, is , as there are a fraction of non-zero interaction strengths drawn from a distribution with mean . So as first approximation, neglecting fluctuations, we can model the local dynamics as:

(9) |

when , and otherwise. Such an effective map is an accurate representation of the population dynamics when there is a high degree of coherence in the system.

To further gauge if the bifurcation diagram obtained numerically in Fig. 4(d) can be understood qualitatively using this effective map picture, we argue that the deviations from the effective dynamics can be modelled by random fluctuations about the mean interaction strength . So we study the dynamics given by Eqn. 9 under the influence of multiplicative noise as well. The results are shown in Fig. 7. It can be clearly seen that the effective dynamical map, under random noise, qualitatively captures the collective dynamics of the multi-species communities.

Analysis: We can also straight-forwardly analyse the effective map dynamics given by Eqn. 9 to find the windows where the positive steady state is stable. Note that this non-trivial fixed point, which is a function of , can be obtained as a solution of:

(10) |

and its stability is determined by the condition . Clearly as , the fixed point becomes unboundedly large. This also explains the presence of the critical positive mean , with , in the system. From Fig. 8A it is clear that the parameter yielding fixed populations in the system is very close that found in the effective map (cf. Fig. 7).

Further we employed another, more accurate, approach to gauge the stability of the synchronized steady state. In this extension, the stability analysis takes into account the entire network by considering the extremal eigenvalues of the connectivity matrix. This yields stability conditions on the fixed point (which is a solution of Eqn. 10) given by:

(11) |

where, and are the average maximum and minimum eigenvalues of the random gaussian matrix, .

From Fig. 8B, one observes that the region where synchronized steady state is stable, namely where lies within the two bounds, corresponds quite closely to . This matches the value observed in simulations very closely (cf. Fig.4), and thus provides an accurate description of the effective collective behavior of the system.

Lastly, consider the scenario that leads the dynamics of the nodes to extinction, namely to . This will happen if , which then will map to . Notice that for populations very close to zero, . So from Eqn. 9, . This implies that the subsequent iterate can become negative if and only if is negative, as is non-negative and if . This suggests why extinctions are seen to arise for .

### Effect of the degree of variability in inter-species interactions

Having gained understanding of the collective dynamics of the system in terms of the dynamics of a single effective stochastic map, we now try to understand the effect of the standard deviation of the connectivity matrix on the dynamics of the multi-species community. It seems reasonable to argue that the strength of the noise term in the effect dynamical map is directly related to , namely we can associate the stochasticity in the effective map to variability in the inter-species interaction strength across the network. Thus we investigate the changes in the bifurcation structure for two different values of , which represents different degrees of spatial variability in the network. From Fig. 9, one can observe that with increasing , the bifurcation diagram gets more noisy. This indicates that one can incorporate the spatial variability in interaction structure easily in the effective map.

Also, as discussed earlier, for populations close to zero, Eqn. 9 effectively gives , and these are driven to extinction if . This implies that the sign of the subsequent iterate for a system close to zero is very sensitive to large fluctuations in the distribution of interaction strengths. Namely, large variability around the mean value in Eqn. 9 can push the system into the extinction zone, or out of it, when , and are close to zero. This accounts for the spread in values around , in the presence of large in Fig. 1.

## Discussions

In summary, we have analyzed the survival properties in ecological networks. In particular, we considered a complex network of populations where the links are given by a random asymmetric connectivity matrix , with fraction of zero entries, where reflects the over-all connectivity of the system. The non-zero elements are drawn from a gaussian distribution with mean and standard deviation . The signs of the elements of the connectivity matrix reflect the nature of density-dependent interactions, such as predatory-prey, mutualism or competition, and their magnitude reflect the strength of the interaction. Unlike many earlier studies, we investigate the full range of mean interaction strengths, and do not confine ourselves to the balanced situation which assumes .

Also note that one can potentially draw a parallel between our model and the system of metapopulations with density dependent dispersal [13]. Namely, our system can also be interpreted as a network of metapopulation patches [14], or “a population of populations” [15]. In particular, it can describe a system comprising many spatially discrete sub-populations connected by migrations where inter-patch dispersal is both high enough to ensure demographic connectivity among patches, yet low enough to maintain some degree of independence in local population dynamics. The connectivity matrix in this scenario reflects density dependent dispersal and migration, as is commonly seen in vertebrate and invertebrate populations [16, 17, 18, 19, 20, 21].

A problem of vital importance here is understanding how broad features, such as the connectedness and net positive interaction strength, modulates the emergent dynamics in such a network. First, in order to gauge the global stability of the system, we calculate the average number of active nodes, namely number of non-extinct species, in the network. We find that the network transitions from a completely extinct system to one where all nodes are active, as the mean interaction strength goes from negative to positive. This transition, marked distinctly by scaling relations, gets sharper with increasing and decreasing . This result has much relevance, as realistic ecosystems are unlikely to have a perfect balance of interactions. So understanding the implications of imbalance in interaction types and strengths in the network (namely ) is important.

Another significant observation is that the total population, reflecting the biomass production in a multi-species community, displays distinct non-monotonic scaling behaviour with respect to the product , implying that survival is dependent not merely on the number of links, but rather on the combination of the sparseness of the connectivity matrix and the net interaction strength. Interestingly, in an intermediate window of positive , the total population is maximal, indicating that too little or too much positive interactions is detrimental to survival. Infact survival is optimal when the network has intermediate net positive connection strengths. Counter-intuitively then, if positive interactions such as mutualistic or symbiotic links are too dominant, its effect ceases to be beneficial and in fact results in reduction of the total population.

At the local level we observed marked qualitative changes in dynamical patterns, ranging from fixed points to spatioteporal chaos, under variation of mean of the interaction strengths. Specifically we found anti-phase clusters of period cycles and the presence of period- chaotic bands, in certain windows of mean . This behaviour is reminiscent of the field experiment conducted by Scheffer et al [12] which showed the existence of self-perpetuating stable states alternating between blue-green alage and green algae. We also studied the correlation between synchronization and survival, and find that synchronization does not necessarily lead to extinction. Lastly, we proposed an effective low dimensional map to capture the behavior of the entire network, and this provides a broad understanding of the interplay of the local dynamical patterns and global stability of the network.

## References

- 1. M. Newman, Networks: An Introduction (Oxford University Press, Oxford, UK, 2010).
- 2. R. May, Nature, 238, 413-414(1972).
- 3. See representative work: B Blasius, A Huppert, L Stone, Nature 399 (1999), 354; T Gross, B Blasius, J. of R. Soc. Interface 5 (2008), 259; M. Baurmann, T. Gross, U. Feudel, J. of Theor. Bio. 245 (2007) 220.
- 4. A. R. Ives & S. R. Carpenter, Science, 317 58-62(2007).
- 5. A. Mougi & M. Kondoh, Science, 337, 349-351(2012).
- 6. S. Sinha and S. Sinha, Phys. Rev. E, 71, (2005) 020902; S. Sinha and S. Sinha, Phys. Rev. E, 74, (2006) 066117.
- 7. R. May, Ecology, 54, 638-641(1973).
- 8. S. Allesina & S.Tang, Nature, 483, 205-208(2012)..
- 9. PA Stephens, WJ Sutherland & RP Freckleton, Oikos(JSTOR), 87, 185 (1999)
- 10. N.K. Kamal and S. Sinha, Chaos, Solitons and Fractals 44 (2011) pp. 71-78
- 11. A. Choudhary, V. Kohar, S. Sinha, Eur. Phys. J. B 87 (2014), 202
- 12. M. Scheffer, S. Rinaldi, A. Gragnani, L. R. Mur, & E. H. van Nes Ecology, 78, 272-282(1997).
- 13. I Hanski, Nature, 396, 41 (1998).
- 14. R Levins, Extinction. In Some mathematical questions in biology (ed. M. Gerstenhaber), pp. 77–107. Providence, RI: The American Mathematical Society (1970).
- 15. R Levins, Bulletin of the ESA, 15, 237 (1969).
- 16. B Kerr, MA Riley, MW Feldman & BJM Bohannan, Nature, 418, 171 (2002).
- 17. B Kerr, C Neuhauser, BJM Bohannan & AM Dean, Nature, 442, 75 (2006).
- 18. S Dey & A Joshi, Science, 312, 434 (2006).
- 19. JH Brown & A Kodric-Brown, Ecology, 58, 445 (1977).
- 20. T Reichenbach, M Mobilia & E Frey, Nature, 448, 1046 (2007)
- 21. B Cazelles, S Bottani & L Stone, Proc. R. Soc., B, 268, 2595 (2001)
- 22. Stachowicz J.J., Mutualism, facilitation, and the structure of ecological communities, Bioscience (2001) 51 235
- 23. Amritkar, R. E., & Rangarajan, G. (2006) Phys. Rev. Letts., 96, 258102.
- 24. Nishikawa, T., & Motter, A. E. (2010). PNAS 107, 10342.
- 25. Sorrentino, F., & Ott, E. (2008) Phys. Rev. Letts., 100, 114101.