Understanding the mechanisms of ecosystem stability is a fundamental issue in ecology (May, 1972; Chesson, 2000; Levine and HilleRisLambers, 2009; Allesina and Levine, 2011). Many proposed mechanisms rely heavily on exogenous factors mitigating the effects of competition, such as niche differentiation (Chesson, 2000; Levine and HilleRisLambers, 2009; Chu and Adler, 2015) and intermediate disturbance hypothesis (Connell, 1978; Roxburgh et al., 2004; Liao et al., 2022). In contrast to these exogenous mechanisms, intransitive competition provides an endogenous mechanism promoting coexistence (Laird and Schamp, 2006; Allesina and Levine, 2011; Soliveres et al., 2015; Levine et al., 2017), and, in particular, the simple rock-paper-scissors game has become a typical case to explain ecosystem stability (Huisman and Weissing, 1999; Kerr et al., 2002; Reichenbach et al., 2007). The classic, cyclic game of rock-paper-scissors usually leads to species’ abundances neutrally cycling without converging to a stable equilibrium point, which is rarely observed in nature. To overcome this problem, many approaches have been proposed to explain the robust persistence of the cyclically competing system, such as higher-order interactions (Grilli et al., 2017) and spatially structured interactions (Durrett and Levin, 1997; Czárán et al., 2002; Rojas-Echenique and Allesina, 2011; Calleja-Solanas et al., 2021; Zhang et al., 2022). In particular, the latter mechanism of local structured interactions stabilizing coexistence has been documented experimentally (Kerr et al., 2002). Thus, it is widely believed that the inclusion of spatial structure, where the interactions and dispersal of individuals are local, can result in stable coexistence in rock-paper-scissors games. Yet, these conclusions are mostly drawn from lattice-based models, where each individual is assumed to only interact with its surrounding neighbors. This assumption is relatively restrictive, as species dispersal and interactions in nature might occur on complex networks with variation in patch connectivities (i.e., spatial heterogeneity in dispersal networks; Urban and Keitt, 2001; Fortuna et al., 2006; Dale and Fortin, 2010; Galpern et al., 2011; Grilli et al., 2015; Fortin et al., 2021; He et al., 2021; Li et al., 2021; Zhang et al., 2021).

There is abundant evidence that landscape structure, and other factors (e.g., patch quality; Liao et al., 2013), can result in anisotropic (i.e., directionally biased) dispersal behavior and therefore spatial heterogeneity in patch connectivities (Urban and Keitt, 2001; Fortuna et al., 2006; Dale and Fortin, 2010; Galpern et al., 2011; Grilli et al., 2015; Fortin et al., 2021; He et al., 2021). Since different patches in the landscape might be perceived differently by species (Hansbauer et al., 2010; Dondina et al., 2018), and the resulting dispersal network would display diverse patterns of patch connectivity (Yeaton and Bond, 1991; Bunn et al., 2000; Nicholson and Possingham, 2006; Fortuna et al., 2009; Bearup et al., 2013; Hirt et al., 2018; Germain et al., 2019). For instance, species dispersal between sub-reefs within the Great Barrier Reef has been characterized with scale-free networks (Kininmonth et al., 2010), while seed dispersal by birds, as opposed by winds, is better described by an irregular network than a spatially uniform network. In addition, Fortuna et al. (2006) identified a large spatial dispersal network of temporary ponds, which are used as breeding sites for amphibian species, following a power-law degree distribution. As such, there has been an increasing interest in exploring the effects of network heterogeneity on ecosystem stability using graph theory (Szabó et al., 2004; Szolnoki and Szabó, 2004; Masuda and Konno, 2006; Dale and Fortin, 2010; Schütt and Claussen, 2010; Galpern et al., 2011; Laird, 2014; Nagatani et al., 2018; Fortin et al., 2021; He et al., 2021; Zhang et al., 2021). Many studies have found that increasing network heterogeneity (i.e., increasing variation in patch connectivities) can promote stable coexistence in cyclic competition (Masuda and Konno, 2006; Schütt and Claussen, 2010; Nagatani et al., 2018).

Despite these advances, several lattice-based models (Laird and Schamp, 2008; Rojas-Echenique and Allesina, 2011; Zhang et al., 2022) have observed that local intransitive competition can reduce species coexistence compared to long-range competition, in stark contrast to current knowledge of local interactions stabilizing coexistence (Durrett and Levin, 1997; Huisman and Weissing, 1999; Czárán et al., 2002; Kerr et al., 2002; Calleja-Solanas et al., 2021). Rojas-Echenique and Allesina (2011) and Zhang et al. (2022) attributed the opposite outcomes to different interaction modes, which can induce distinct stabilizing mechanisms in lattice-structured models. Indeed, species in diverse natural ecosystems might interact in different ways. For example, there are two typical interaction modes often observed in nature: seedlings of tree species or propagules of grass species compete to fill gaps, while animal species fight directly for colony sites (Rojas-Echenique and Allesina, 2011; Calleja-Solanas et al., 2021; Zhang et al., 2022). However, the effects of these different interaction modes, in combination with spatial heterogeneity in dispersal networks mentioned above, have not been well integrated into our general understanding of how stable coexistence emerges in cyclic competition. In this study, we thus incorporate dispersal network heterogeneity into the classic rock-paper-scissors games, in order to make a comparative analysis of ecosystem stability between the two typical interaction modes.

We begin our analysis by inspecting the temporal evolution of species abundances (Figure 2). When species compete to fill an empty patch (Mode 1), increasing network heterogeneity (from lattice-structured regular to scale-free networks) increases the magnitude of population fluctuations and speeds up stochastic species extinctions. In contrast, if the superior competitor directly substitutes the inferior one (Mode 2), all three species can coexist, and the magnitude of population fluctuations around the equilibrium point (≈1/3) decreases as network heterogeneity increases. In lattice-structured regular networks, both interaction modes yield similar fluctuation size (Figure 2A vs. Figure 2E with different ranges of y-axis). This is more clear when looking at their dynamic trajectories (Figure 2I vs. Figure 2M). In these dynamic trajectories (Figures 2I–P), the abundances of the species at each generation represent a point in the 3-simplex, whose vertices correspond to a monospecific population. Over time, the point follows a trajectory which will reflect the state of the system. If the system is close to the equilibrium abundances, the trajectory eventually occupies a small area near it, whilst larger fluctuations cover larger areas. Similar to Figures 2A–H, the interaction mode 1 results in greater fluctuations and ultimately monoculture in more heterogeneous networks (Figures 2I–L), as opposed to the interaction mode 2 where ecosystems are more stable in networks with higher heterogeneity.

To test whether these outcomes are robust, we simulate 20 replicates for each case, by regenerating the dispersal network in each replicate (Figure 3). In Mode 1, three species can coexist with oscillations around the equilibrium point in the lattice-structured regular network (see Figure 2A), unlike other three heterogeneous networks where only one species ultimately dominates the system (Figure 3A). As such, the outcome in the regular network is absent in Figure 3A. In other three heterogeneous networks, we record the extinction time of the first species in each replicate, with more time required for species extinction implying more stable coexistence. As shown in Figure 3A, increasing network heterogeneity generally leads to species going extinct sooner. In Mode 2, all species can coexist with fluctuations around the equilibrium point in all dispersal networks. Thus, we use variation in fluctuations size (C.V. around the mean species abundances at steady state) to characterize ecosystem stability, with lower values of C.V. indicating more stable coexistence. As shown in Figure 3B, an increase in network heterogeneity generally results in smaller fluctuations in species abundances, yielding a more stable system. These outcomes further confirm our previous conclusions in Figure 2.

To better understand the mechanisms behind these dynamic behaviors, we display several snapshots of the spatial patterns of the three competitors in complex networks (Figure 4). In the lattice-structured regular network, both interaction modes produce similar conspecific clumping patterns (Figures 4A,B). To clearly illustrate the spatial organization in dispersal networks with contrasting heterogeneities, we take a system of N = 100 patches for example by displaying their snapshots at the 5th generation in Mode 1 and the 20th generation in Mode 2 (Figures 4C–J). In Mode 1 (Figures 4C–F), when three species have almost equal abundance in the regular network, increasing network heterogeneity which tends to increase network modularity, enlarges differences in species abundances, most likely accelerating stochastic extinctions. In Mode 2 (Figures 4G–J), if three species have unequal abundances in the regular network, dispersal network heterogeneity (i.e., network modularity) acts as a driving force to equalize their abundances, thereby stabilizing the ecosystem. We also observe that species can form self-organized conspecific clusters with the most connected patches at the core in these heterogeneous networks.

Under both interaction modes, we finally focus on how system size (i.e., increasing network size) affects fluctuations in species abundance in networks with contrasting heterogeneities (Figure 5). Intuitively, increasing network size promotes ecosystem stability, regardless of interaction mode and network heterogeneity. Specifically, stochastic extinctions in Mode 1 require more time to occur as network size increases, especially in networks with less heterogeneity (Figure 5A). This also demonstrates that increasing network heterogeneity greatly destabilizes the system. In Mode 2, the size of population fluctuation declines linearly with increasing network size (with log-log scale in Figure 5B), thereby promoting ecosystem stability. Furthermore, increasing network heterogeneity leads to smaller fluctuations, further confirming that network heterogeneity can stabilize the cyclically competing system.

Incorporating complex networks into the classic rock-paper-scissors game, we find that whether network heterogeneity can stabilize competitive dynamics depends on the interaction mode. Specifically, if species compete directly for colony sites, increasing spatial heterogeneity in dispersal networks stabilizes the cyclically competing ecosystems, further confirming previous theoretical arguments (Masuda and Konno, 2006; Schütt and Claussen, 2010; Nagatani et al., 2018). In contrast, when species compete to fill an empty site, an increase in network heterogeneity leads to stronger population fluctuations and thus increases the risk of stochastic extinctions. Interestingly, in the lattice-structured regular network, both interaction modes display similar coexistence patterns based on the same mechanism: local interactions allow species to survive by forming conspecific clusters (via self-organization) where interspecific competition only takes place at the borders between heterospecific clusters, thereby decreasing the effective interspecific competition (by reducing interspecific encounter rate) and therefore population fluctuations (Figures 4A,B).

There exist different mechanisms between the two interaction modes that induce contrasting effects of network heterogeneity on ecosystem stability. In our model, the stabilizing role of network heterogeneity is observed when competition occurs through the superior competitor displacing the inferior one in a randomly chosen pair of neighbors. This competition mode has been commonly used and is symmetric in a cyclical way. However, increasing variation in patch connectivities implies that a few patches are highly connected while most have only few connections. A patch with more connections is more likely to be selected as a neighbor, while another neighbor can be sampled from multiple patches that might accommodate any of the three species. As such, which species can occupy this highly connected patch will change more frequently in networks with higher heterogeneity. To some extent, this prohibits the growth of self-organized conspecific clusters with the most connected patches at the core. In fact, increasing network heterogeneity breaks the symmetry in species interactions, that is, those individuals in highly connected patches have a higher chance to participate in competition than other individuals in poorly connected patches. This results in asymmetric interactions that can produce negative frequency dependence (NFD) which is absent in symmetric interactions (Rojas-Echenique and Allesina, 2011). Specifically, the NFD can decrease the average fitness of individuals when they become more common (cf. Zhang et al., 2022). Thus, the NFD, which is more significant in more heterogeneous networks, can suppress the population growth of dominant species but promote the growth of rare species, thereby increasing the frequency of oscillations and stabilizing competitive dynamics.

Network heterogeneity is destabilizing for the interaction mode where species compete for an empty site. This is because, in this mode, each species has the same probability to be selected to die, regardless of patch connectivities, so NFD cannot emerge. As such, conspecific clusters, with highly connected patches at the core, can grow (via self-organization) without restriction especially in networks with higher heterogeneity. Within these conspecific clusters, a competition event does not contribute to variation in species abundance, as local competition occurs between conspecific individuals. Thus, changes in species abundances can only take place along the borders between heterospecific clusters. However, the length of these borders increases more rapidly in more heterogeneous networks (i.e., a rapid increase in the number of “active” individuals that can change the state at the borders), promoting interspecific encounter rate and therefore interspecific competition. In addition, these heterospecific clusters can be treated as compartments dominated by different species, but there are lots of links connecting these compartments. If a superior competitor invades a compartment dominated by an inferior one, then it can colonize this compartment via the highly connected core. Ultimately, this results in larger population fluctuations and makes stochastic extinction events more likely.

This study incorporates complex networks into the rock-paper-scissors game under two commonly observed modes of interaction. In the lattice-structured regular network, both interaction modes yield similar system stability through forming conspecific clusters to reduce interspecific competition (Figure 4A vs. Figure 4B). Interestingly, we find that the inclusion of network heterogeneity can induce contrasting coexistence patterns between the two interaction modes, due to different mechanisms as explained above. This strongly suggests that particular attention should be devoted to testing, theoretically and experimentally, which mode of interaction is more appropriate for modeling a given competing system. For example, whereas the death of a tree or grass creates a gap in plant communities (i.e., Mode 1), animal species more typically fight for territory (i.e., Mode 2). Thus, these communities should be modeled with different interaction modes and so can be expected to be affected differently by the structure of the landscape they inhabit. Overall, even if these outcomes are obtained using an extremely simplified model, our findings can help identify different mechanisms to explain the role of different interaction modes in stabilizing competitive dynamics in complex networks with contrasting heterogeneities.

The original contributions presented in this study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

GG, ZZ, and HZ built the model, performed simulations, and analyzed the results. JL conceived the study and wrote the manuscript. DB contributed substantially to revisions. All authors contributed to the article and approved the submitted version.