Evolutionary Game Dynamics with Environmental Feedback in a Network with Two Communities

Recent developments of eco-evolutionary models have shown that evolving feedbacks between behavioral strategies and the environment of game interactions, leading to changes in the underlying payoff matrix, can impact the underlying population dynamics in various manners. We propose and analyze an eco-evolutionary game dynamics model on a network with two communities such that players interact with other players in the same community and those in the opposite community at different rates. In our model, we consider two-person matrix games with pairwise interactions occurring on individual edges and assume that the environmental state depends on edges rather than on nodes or being globally shared in the population. We analytically determine the equilibria and their stability under a symmetric population structure assumption, and we also numerically study the replicator dynamics of the general model. The model shows rich dynamical behavior, such as multiple transcritical bifurcations, multistability, and anti-synchronous oscillations. Our work offers insights into understanding how the presence of community structure impacts the eco-evolutionary dynamics within and between niches.


Introduction
Evolutionary game theory is the study of population changes driven by competition among different strategies.A recent adjustment of evolutionary game models with the aim of better representing the natural world is the inclusion of strategy-dependent feedback, specifically, environmental feedback [1].This type of game is called an eco-evolutionary game.This type of strategy-dependent feedback can be seen in many complex systems, such as ecological metacommunities [2], collectives of insect individuals [3,4], microbial populations [5][6][7], and human social and reproductive structures [8,9].A major question with models of ecoevolutionary game dynamics is conditions under which cooperation in a population can thrive when the payoff matrix, which we regard as the environment, is influenced by the action of players.Extensions of the original eco-evolutionary game dynamics models include the addition of finite carrying capacity [10], renewable and decaying resources [11][12][13], imitation and aspiration dynamics [14], mutation of players [15], reciprocity dynamics [16], and extension to public goods games [12,[17][18][19].The models can also be extended in terms of additional types of dynamic feedback, such as non-constant enhancement or degradation rates of the environmental variable, which depends on the payoff of players [20], and global and local environment fluctuations [18].
Given that players of the game are embedded in structured populations in reality, evolutionary game models have been extended to the case of various networks [21][22][23][24].Similarly, players involved in an eco-evolutionary game may be better interpreted to inhabit on nodes of a network.Therefore, eco-evolutionary games have been extended to the case of networks.For example, in eco-evolutionary games on regular graphs, it was found that a higher degree of the node creates oscillatory behavior in the population and that a lower degree promotes spread of cooperation [25,26].Spatial networks are also commonly used for exploring how environmental feedback promotes cooperation [27][28][29][30][31][32][33][34][35][36][37].Lastly, through the use of bimatrix payoffs, which are equivalent to the complete bipartite graph as population structure in the case of symmetric payoff matrices, periodic orbits in the state space have been proven to exist [38][39][40][41].
However, there are some vital gaps missing in the prior research on eco-evolutionary games on networks.First, in complete bipartite graphs [38][39][40][41], the players do not interact within each community.This assumption is suitable for modeling situations in which the population of players is divided into two different roles but otherwise not in general.Second, in most of the previous studies, the environmental state is assumed to be either a globally shared variable [25,26,[38][39][40][41] or local to each node (i.e., player) [27][28][29][30][31][32][33][34][35][36][37].However, it may be more realistic to assume that the environment is shared across some, but not all, players [42,43].For example, a meta-community in ecological systems may be an appropriately sized unit for considering an environmental variable [2,[44][45][46].Other eco-evolutionary game models assume network structure and assign a local environmental variable to each edge between a pair of players [35,37].
In the present study, we extend a previously proposed model of eco-evolutionary dynamics [1] to the case of networks with equally sized two communities.Unlike the complete bipartite graph models proposed in [38,39,41] where the players in each community only interact with those in the other community, we assume that players not only interact with those in the other community but also with those in the same community.Next, we assume that the state of the environment depends on the type of edge in the network, similarly to [35,37].We crucially assign one environmental variable to each type of edge, i.e., the edges within the first community, those within the second community, and those connecting the two communities.
Figure 1: Schematic of the two-community network.A filled circle represents a player.Two players from the same community interact at rate 1 − δ.Two players from the opposite communities interact at rate δ.Without loss of generality, we normalize the rate parameter 0 < δ < 1.We only show some edges for visualization purposes.
In this manner, we model the situation in which two players forming an edge may improve or deteriorate their shared environment, which is assumed to be on the edge.We do not distinguish between edges of the same type because of the symmetric population structure assumed.Unlike the previous studies similarly assuming edge-dependent environmental states [35,37], our two-community network model, which is a minimal network model, allows analytical investigations.
Our paper is organized as follows.In Section 2, we describe our model in detail and focus on eco-evolutionary dynamics with two network communities.In Sections 3 and 4, respectively, we present our stability analysis of the simplified replicator dynamics resulting from different symmetry assumptions.In Section 5, we numerically investigate the rich dynamical behavior of the general model.Finally, we discuss contributions of the current work along with an outlook for future work.

Model
Consider an eco-evolutionary game in a population composed of two communities.Each player chooses either of the two actions, i.e., cooperation or defection.We assume that there are N players in total and N/2 players in each community.We assume that the entire population is infinite (i.e., N → ∞) and that the players interact with each other player within the same community at rate 1 − δ and with each player in the other community at rate δ > 0. See Fig. 1 for a schematic.
We consider replicator dynamics for a population on the two-community network with feedback-evolving games.Crucially, we assume that the state of the environment depends on the type of edge in the network.We denote by n 1 ∈ [0, 1] the state of the environment in community 1, representing the edges within community 1, by n 2 ∈ [0, 1] the state of the environment in community 2, and by n 12 ∈ [0, 1] the state of the environment used when a player in community 1 and one in community 2 interact.The environment-dependent payoff matrices for community 1, 2, and in between are assumed to be given by where n is either n 1 , n 2 , or n 12 .We assume that, if n = 0, then cooperation is the unique Nash equilibrium, i.e., R 0 > T 0 and S 0 > P 0 .If n = 1, then defection is the unique Nash equilibrium, i.e., R 1 < T 1 and S 1 < P 1 .We label the prior inequalities as Let us define q 1 and q 2 as the two-dimensional payoff vector for a player in community 1 and 2, respectively.The first entry of the vector is the payoff for a cooperator.The second entry of the vector is the payoff for a defector.Define x and y as the fraction of cooperators in community 1 and 2, respectively.The fraction of defectors in community 1 and 2 is 1 − x and 1 − y, respectively.We obtain where x = x 1 − x ⊤ , y = y 1 − y ⊤ , and ⊤ denotes the transposition.The first term on the right-hand side of Eqs. ( 3) and ( 4) is the payoff obtained by playing with the other players in the same community.The second term is the payoff obtained by playing with the players in the opposite community.
We assume that the competition between cooperation and defection occurs only within each community because players inhabiting different communities may perceive the different environments due to the different state of the environment.Then, the replicator dynamics are given by ẋ where q 11 , q 12 , q 21 , and q 22 are defined by with i ∈ {1, 2}.We give the dynamics of the environmental state of each type of edge by where z is the fraction of cooperators in the entire population, i.e., z ≡ (x+y)/2, and θ 1 > 0, θ 2 > 0, and θ 12 > 0 are the ratio of enhancement to degradation of the environmental variable for the respective edge type.For example, if θ 1 is large, then enhancement of the environment in community 1 occurs at a relatively small fraction of cooperators, x.
3 Three-dimensional system with θ 1 ̸ = θ 12 In this section, we assume that θ 1 = θ 2 , and that the initial condition satisfies x = y and n 1 = n 2 .Then, x = y and n 1 = n 2 hold true for any t > 0. We further assume that θ 1 = θ 2 ̸ = θ 12 .In this case, the original five-dimensional dynamical system is reduced to the three-dimensional dynamical system given by ẋ We analyze the equilibria and dynamics of this three-dimensional dynamical system.The Jacobian of this dynamical system is given by where and q 11 and q 12 are given by Eq. ( 7).

Corner equilibria
We denote by x * the equilibrium of x and similar for the other dynamical variables.By setting x * , n * 1 , and n * 12 to 0 or 1, specifying the corners of the unit cube defined by 0 ≤ x, n 1 , n 12 ≤ 1, we obtain 8 corner equilibria.We show in Appendix A that each corner equilibrium is a saddle.

Face equilibria
In this section, we seek equilibria on the face of the unit cube, i.e., those in which just one of x * , n * 1 , or n * 12 is either 0 or 1 and the other two are between 0 and 1.We call these equilibria face equilibria.Similarly to the case of the edge equilibria, if we let x * = 0 or 1, then we obtain a corner equilibrium.Therefore, we assume that 0 < x * < 1.By setting just one of n * 1 or n * 12 to 0 or 1, we obtain the four face equilibria shown in Table 1.

Movement of stable equilibria as δ varies
The results in sections 3.1-3.4indicate that, for given θ 1 and θ 12 (̸ = θ 1 ) values, there are three equilibria, two of which are face equilibria and one is an edge equilibrium.Just one of these three equilibria is stable for a given value of δ.Specifically, when θ 1 < θ 12 , a face equilibrium is stable when 0 < δ < δ c,1 , an edge equilibrium is stable when δ c,1 < δ < δ c,2 , and another face equilibrium is stable when δ c,2 < δ < 1; see Fig. 3(a).As δ varies, the position of the stable equilibrium continuously moves, including through δ = δ c,1 and δ = δ c,2 .The dynamical system undergoes a transcritical bifurcation at δ = δ c,1 , with which the face equilibrium and the edge equilibrium exchange the stability.Another similar transcritical bifurcation occurs at δ = δ c,2 .See Figs.4(a) and 4(b) for visualization.When θ 1 > θ 12 , a different set of three equilibria, which reside on the opposite side of the unit-cube state space, are stable for a respective range of δ, as shown in Fig. 3(b).Similarly to the case of θ 1 < θ 12 , these equilibria undergo transcritical bifurcations at δ = δ c,3 and δ c,4 .
We point out that, as the transcritical bifurcation is approached as δ gradually increases 4 Three-dimensional system with θ 1 = θ 12 In this section, as in section 3, we assume that θ 1 = θ 2 and that the initial condition satisfies x = y and n 1 = n 2 .Then, x = y and n 1 = n 2 hold true for any t > 0. We now further assume that θ 1 = θ 2 = θ 12 .

Corner equilibria
By setting x * , n * 1 , and n * 12 to 0 or 1, we obtain eight corner equilibria.Similar to the case of θ 1 ̸ = θ 12 (see section 3.1), each corner equilibrium is a saddle.See Appendix D for the proof.

Interior equilibria
In this section, we look for equilibria in the interior of the unit cube, i.e., those satisfying 0 < x * , n * 1 , n * 12 < 1.By setting ṅ1 = 0 and ṅ12 = 0 in Eqs. ( 14) and ( 15), respectively, with θ 1 = θ 12 , and imposing n * 1 , n * 12 / ∈ {0, 1}, we obtain By substituting Eq. ( 42) in Eq. ( 13) and imposing ẋ = 0, we obtain Any point on this line is an equilibrium.We call Eq. ( 43) the line of equilibria and denote it by L; it is the equilibrium manifold.We show in Appendix E that L is neutrally stable along the direction of L and that the other two eigenvalues, λ 2 and λ 3 , have negative real part if Eqs. ( 2) and ( 19) hold true.In this case, line L attracts trajectories near L.
To demonstrate L, we numerically simulate trajectories with θ 1 = 5 and δ = 0.5, for which λ 2,3 = −0.019± 0.810i.We show trajectories of the dynamics starting from two initial conditions in Fig. 6.The figure indicates that the solution spirals into L as expected.

Edge equilibria
Let us examine possible edge equilibria.It should be noted that ρ 1 = ρ 12 when θ 1 = θ 12 ; we recall that ρ 1 and ρ 12 are defined in Eqs. ( 24) and ( 25), respectively.We find that there are just two edge equilibria when θ 1 = θ 12 , which are the same as those found for the case θ 1 ̸ = θ 12 in section 3.3.These two edge equilibria occur where line L intersects the edge specified by n * 1 = 0, n * 12 = 1 or that specified by n * 1 = 1, n * 12 = 0. We show in Appendix F that the edge equilibrium (x , 0, 1 is marginally stable with two zero eigenvalues and one negative eigenvalue if and only if Eq. ( 19) holds true and Figure 6: System's behavior near the equilibrium manifold L. Shown are trajectories of the three dimensional system given by Eqs. ( 13), ( 14), and ( 15) when θ 1 = θ 12 = 5 for two initial conditions.The green line indicates L, the line of equilibria given by Eq. ( 43).We use the payoff matrices given by Eq. ( 38), initial conditions (x, n 1 , n 12 ) = (0.5, 0.4, 0.1), shown in blue, and (0.1, 0.9, 0.9), shown in orange, and set δ = 0.5.
When δ ̸ = δ c,1 , the Jacobian has two positive eigenvalues and one negative eigenvalue.Similarly, the edge equilibrium , 1, 0 is marginally stable if and only if Eq. ( 19) holds true and When δ ̸ = δ c,3 , the Jacobian has two positive eigenvalues and one negative eigenvalue.
To understand the location of the face equilibria depending on the value of δ, we examine the movement of line L on the (n 1 , n 12 ) plane as we vary δ.The two intersections of L with the boundary of the square defined by 0 ≤ n 1 , n 12 ≤ 1, combined with x * = 1 1+θ 1 , give the two face equilibria.When the intersection is at a corner of the square, it is an edge equilibrium.We show L as a function of δ in Fig. 7 for the payoff matrices given by Eq. (38). Figure 7 indicates that the two edge equilibria are realized at different δ values, which is consistent with the results shown in section 4.3.The figure also indicates that L passes through a particular point regardless of the δ value.By setting both the coefficient of δ and the constant term to 0 in Eq. ( 43), we obtain this point as follows: Figure 7 also indicates that, when δ is small, n * 1 is highly variable between 0 and 1, but the range of n * 12 is small.When δ is large, the converse is true.This result is natural because a larger δ implies that more interaction between players occur between the two communities than in the same community.
As δ varies, our three-dimensional dynamical system undergoes two bifurcations at δ = δ c,1 and δ = δ c,3 .When 0 < δ < δ c,1 , the face equilibrium with n * 12 = 1 is stable except along the direction of L (therefore, the Jacobian has two negative eigenvalues and one 0 eigenvalue), and the edge equilibrium given by (x * , n * 1 , n * 12 ) = P 0 −S 0 −δ(P 0 −P 1 −S 0 +S 1 ) R 0 −T 0 −S 0 +P 0 −δγ , 0, 1 and the face equilibrium with n * 1 = 0 are saddles (when disregarding the 0 eigenvalue along the direction of the line of equilibria; same in the following text).When δ = δ c,1 , the dynamical system undergoes a transcritical bifurcation and the stability of the two face equilibria switches.At δ = δ c,1 , the edge equilibrium has two 0 eigenvalues and one negative eigenvalue.These three equilibria collide at δ = δ c,1 , which we depict in Fig. 4(c) and (d).When δ c,1 < δ < 1, the face equilibrium with n * 12 = 1 and the edge equilibrium given by , which owes to Eq. ( 46).
, 0, 1 are saddles, and the face equilibrium with n * 1 = 0 is stable.There are three other equilibria located at the other end of L intersecting a face or edge of the state space, i.e., the unit cube.The structure of the bifurcation occurring at δ = δ c,3 , involving this second triplet of equilibria, which are composed of two face equilibria (one with n * 1 = 1 and the other with n * 12 = 0) and one edge equilibrium given by (x * , n * 1 , n * 12 ) = P 1 −S 1 +δ(P 0 −P 1 −S 0 +S 1 ) , 1, 0 , is qualitatively the same.Similarly to when θ 1 ̸ = θ 12 , as δ gradually increases from 0 to approach the first transcritical bifurcation, the two eigenvalues except the 0 eigenvalue are first complex conjugates with negative real parts and then change to real negative values.Figure 5(c) shows the dependence of the real part of the two eigenvalues on δ around δ = δ c,1 .Therefore, when L intersects the unit cube at a point not close to an edge, trajectories on the face spiral into the stable face equilibria, which is consistent with the numerical results shown in Fig. 6.When the stable face equilibrium approaches an edge of the unit cube, it becomes a sink, enabling the transcritical bifurcation on the edge.The dependence of the Jacobian eigenvalues of the three equilibria near δ = δ c,3 is qualitatively the same as that near δ = δ c,1 (see Fig. 5(d)).

Five-dimensional system
In this section, we analyze the five-dimensional dynamical system given by Eqs. ( 8)- (12) without assuming symmetry between the two communities.We exhaustively examine its equilibria as follows.First, we search for all possible combinations of x, y, n 1 , n 2 , and n 12 by classifying the value of each variable to be either 0, 1, or between 0 and 1.Because three options are available for each variable, there are 3 5 = 243 possible combinations.Second, we find that the 2 5 = 32 corners of the state space given by x, y, n 1 , n 2 , n 12 ∈ {0, 1} are equilibria, more specifically, saddles.Third, out of the remaining 211 combinations, we have found that 60 combinations are equilibria; the other 151 combinations are not.We show these equilibria in Appendix G.By analyzing the Jacobian of the 60 equilibria with the assistance of Mathematica, we find that 21 of them are stable under some conditions (see Appendix G).
Figure 8(c) shows an oscillatory trajectory for θ 1 = 3, θ 2 = 5, θ 12 = 8, δ = 0.31, and initial condition (x, y, n 1 , n 2 , n 12 ) = (0.1, 0.5, 0.1, 0.9, 0.5).The inset of the figure, showing the time courses of x and y, indicates anti-synchronization behavior during the oscillatory dynamics.We point out the environmental state between the two communities is bountiful (i.e., n 12 ≈ 1) and almost constant despite the anti-synchronous dynamics between x and y.When one increases δ to δ = 0.4, with all the other parameter values being the same as those used in Fig. 8(c), the oscillations become apparently aperiodic while keeping anti- synchronous behavior between x and y (see Fig. 8(d)).We observe n 12 ≈ 1 and n 2 ≈ 0 during this apparently aperiodic dynamics.It should be noted that n 1 is similarly aperiodic.

Discussion
We extended a previously proposed model of eco-evolutionary dynamics [1] to the case of networks with two equally sized communities.In the three-dimensional dynamical system given by Eqs. ( 13), (14), and (15), which assumes symmetry between the two communities, a further assumption that n 1 = n 12 lends the model the same as the original well-mixed population model [1], and the requirement for the stability of equilibria, i.e., Eq. ( 19), is the same as that derived in [1] as well.
Under the generic condition n 1 ̸ = n 12 , our stability requirement for the equilibria again contained that of [1], i.e., Eq. ( 19).However, the stability of the equilibria in our model also requires conditions on the edge weight between two communities, i.e., δ, and on environment recovery rates, i.e., θ 1 (= θ 2 ) and θ 12 .When θ 1 = θ 12 , the line of equilibria, L, only requires Eq. ( 19) for stability, but the position of L depends on θ 1 and δ.This result implies that the network has no effect on the stability requirements when θ 1 = θ 12 .In contrast, when θ 1 ̸ = θ 12 , the network and the environment recovery rates affect the stability of the system.As a remark, it was mathematically found [15] that the eco-evolutionary dynamical system proposed in [1] has no limit cycles.This mathematical result corroborates with the theoretical results in [1], in which it was proven that the oscillations converge to a heteroclinic cycle, and our numerical results; because we have analytically shown that there is no internal unstable equilibrium, it is unlikely that our system has a limit cycle.
There exists another commonly explored family of dynamic payoff matrices dependent on environmental feedback, given by where T > R and P > S [1,25,26,[38][39][40].With Eq. ( 47), we retain mutual cooperation as a Nash equilibrium when n = 0 and mutual defection when n = 1.In addition, this payoff matrix causes Eq. ( 19) to be satisfied with equality.By using this payoff matrix and holding the assumption that θ 1 = θ 12 , it is straightforward to analytically obtain a neutrally stable interior line of equilibria, which implies closed periodic orbits in the interior of the state space, corroborating the results in [1].When θ 1 ̸ = θ 12 , our system with Eq. ( 47) in fact shows a closed periodic orbit on a face of the hypercubic state space.Therefore, we claim that the closed periodic orbits found in the previous studies with Eq. ( 47) are at least partially due to the symmetry in the payoff matrix given by Eq. ( 47).In the absence of such a symmetry, our results suggest that convergence to stable equilibria is a norm regardless of the population structure.When we removed the assumption of symmetry between the two communities by allowing θ 1 ̸ = θ 2 , we obtained a rich repertoire of stable equilibria, some of which coexist to realize multistability, especially when δ is large.Multistability was also found in other eco-evolutionary models [10,11], but these models are ecological extensions of [1] and are not network-based models as our model is.Bistability was also found in a spatial eco-evolutionary model [30], but for the trivial equilibria (i.e., bistability between an equilibrium with no cooperators in a replete environment and an equilibrium only with cooperators in a rich environment) and under the snowdrift game.In contrast to these previous studies showing multistability in eco-evolutionary game dynamics, our model is a direct network extension of the original model proposed in [1] and without additional ecological assumptions.The present results suggest that multistability may be commonly found in the same eco-evolutionary model on various networks.We also found anti-synchronization behavior during oscillatory population dynamics.This type of behavior was found in a prior complete bipartite graph model [40], but for the division of labor game rather than the typical prisoner's dilemma game.When our stability requirements are not satisfied, our system may converge to a heteroclinic cycle.Further exploring different types of oscillatory behavior in networked eco-evolutionary game dynamics may be interesting.
We emphasize that our model substantially varies from the previously proposed model composed of two interacting subpopulations, or precisely, complete bipartite graphs [38,39,41].Their model does not allow interaction between players in the same subpopulation, whereas our model does.Furthermore, these previous studies adopted the dynamic payoff matrix given by Eq. ( 47), which led to closed periodic orbits, as we discussed above.In [38,39], such cyclic orbits do not accompany anti-synchronous oscillation of the fraction of cooperation in the two subpopulations.Instead, the cyclic behavior originates from interplay of the fraction of one of the two subpopulations and the environmental variable.On the other hand, the orbits obtained in [41] show largely in-phase synchronous oscillation between the two subpopulations.The model in [38,39] was extended in [40] to include a different form of A(n) and different influences of strategies in two subpopulations on the environment.The inclusion of these parameters produces periodic orbits as did the models proposed in [38,39].In contrast, our model showed anti-phase oscillations in terms of the fraction of cooperators in the two communities (i.e., x and y) and multistability.Therefore, even within the family of two-subpopulation networks, which is one of the simplest network model, qualitatively different dynamical behavior may arise depending on the assumption on the environmental dynamics.
Prior extensions of the eco-evolutionary game models to larger networks include those to spatial lattices and regular graphs.The spatial extensions have been to the case of square lattices [27][28][29][30][31][32][33][34][35][36][37].A lattice model of eco-evolutionary game dynamics assuming local environmental variables, meaning that each node (i.e., player) has its own dynamical environmental state, resulted in spatiotemporal patterns, including clustering, flickering, and wave-like pat-terns [31].Enhanced cooperation due to the environmental feedback was also found in eco-evolutionary models on square lattices [27-29, 32, 34-37].Another type of network that has been studied with eco-evolutionary feedback is regular graphs, in which all nodes have degree k.Through the use of pair approximation, the extension of the original model [1] to regular graphs (therefore using the payoff matrix given by Eq. ( 47)) has clarified that an increased k induces the internal stable equilibrium to become neutrally stable, producing periodic orbits [25,26].These models are substantially different from ours not only in the network structure but also in that their model assumes that the environment is global to all nodes.Assigning an environmental state n ij to each edge (i, j), as has been done for square lattices in previous studies [35,37] and for a two-community network in the present study, in the case of regular graphs and general networks may be an interesting generalization.
In addition to the extension of the network structure, edge-dependent environmental state variable, and weighted networks, which we discussed above, there are further possible extensions of the present model as future work.First, in well-mixed populations, incorporation of intrinsic environmental dynamics, such as resource growth and decay, results in multistability and limit cycles [11], which one can explore for networks.Second, the incorporation of dynamic recovery and degradation rates for the environmental state, which are boosted by cooperators' and defectors' payoffs [20], leads to the same stability requirement as that in [1], i.e., Eq. ( 19).One can extend the present model to the case of dynamic rates of environment recovery and degradation by letting, e.g., θ 1 depend on x and n 1 .Third, the use of finite carrying capacity in an environment, which excludes any periodic orbits and enables bistability in the original model [10], should be possible.Fourth, the incorporation of aspiration dynamics, with which players update their strategies based on whether or not they are satisfied with their current payoff [14] is another possible direction of research.Lastly, although we studied the prisoner's dilemma, as other eco-evolutionary game dynamics models, our model can be studied for other games such as the prisoner's dilemma with voluntary participation [40,47], coordination game [1,10,31,41], anti-coordination game [1,10,31], and division-of-labor game [40].
In conclusion, we have studied an eco-evolutionary game dynamics model with two distinct network communities.We find that the interaction rates both within and between these communities significantly impact on the resulting dynamical behavior and the determination of possible equilibrium classes (i.e., interior, face, edge, and corner) of the system.In addition to numerical investigation of the full model, we have performed comprehensive stability analysis of the simplified system under symmetry conditions.Our work highlights the importance of community structures in impacting eco-evolutionary dynamics across different ecological niches.
A Corner equilibria of the three-dimensional system with θ 1 ̸ = θ 12 By evaluating Eq. ( 16) at each corner equilibrium, we obtain By Eq. ( 2), we obtain that S 0 −P 0 > 0, T 0 −R 0 < 0, S 1 −P 1 < 0, and T 1 −R 1 > 0. Therefore, each of these Jacobians has at least one positive eigenvalue and one negative eigenvalue, and each corner equilibrium is a saddle.
For positive δ values satisfying Eq. ( 65) to exist, it must hold true that R 0 − T 0 > S 0 − P 0 .Eigenvalues J 22 and J (2) 33 are negative if and only if and respectively.Lastly, eigenvalue J (2) holds true.However, we find that there is no δ value that simultaneously satisfies Eqs. ( 65), , or one that simultaneously satisfies Eqs. ( 65), ( Case 2: γ > 0 and µ 1 < 0 If γ > 0, then µ 1 < 0 if and only if For positive δ values satisfying Eq. ( 70) to exist and be less than 1, it must hold true that 22 and J 33 are negative if and only if and respectively.For a δ value satisfying Eqs. ( 71) and (72) to exist, it must hold true that θ 12 > θ 1 .
For positive δ values satisfying Eq. ( 86) to exist and be less than 1, it must hold true that R 0 − T 0 > S 0 − P 0 .
Eigenvalues J 22 and J (3) 33 are negative if and only if and respectively.Lastly, J 11 is negative if either holds true.However, we find that there is no δ value that simultaneously satisfies Eqs. ( 86), Case 2: γ > 0 and µ 2 < 0 If γ is positive, then µ 2 < 0 if and only if For positive δ values satisfying Eq. ( 91) to exist, it must hold true that R 1 − T 1 < S 1 − P 1 .Eigenvalues J 22 and J 33 are negative if and only if and respectively.For a δ value satisfying Eqs. ( 92) and (93) to exist, it must hold true that θ 12 < θ 1 .
Case 4: γ < 0 and µ 2 < 0 In this section, we assume that γ < 0 and µ 2 < 0, which requires Eq. (86).For positive δ values satisfying Eq. ( 86) to exist and be less than 1, it must hold true that R 0 −T 0 > S 0 −P 0 .The derivation of the stability of this case follows the same derivation as Case 2, and we find that the equilibrium is stable if and only if Eqs. ( 19), (26), and ( 27) hold true.
C Three face equilibria of the three-dimensional system with θ 1 ̸ = θ 12 In this section, we derive the stability conditions for three face equilibria of the threedimensional system with θ 1 ̸ = θ 12 .
C.1 Equilibrium , the Jacobian is reduced to where The characteristic equation is given by (4) 31 Eigenvalue is negative if and only if θ 12 < θ 1 (i.e., Eq. ( 27)).The real part of the other two eigenvalues is negative if and only if −J 31 > 0, which holds true if and only if δ > δ c,4 (i.e., Eq. ( 39)).
, the Jacobian is reduced to where The characteristic equation is given by 13 is negative.Therefore, −J 31 > 0, which holds true if and only if δ > δ c,2 (i.e., Eq. ( 41)).
D Corner equilibria of the three-dimensional system with θ 1 = θ 12 By evaluating Eq. ( 16) at each corner equilibrium, we obtain the same Jacobians as those in Appendix A with θ 1 = θ 12 .The stability analysis of these corner equilibria is the same as that in Appendix A, and we find that each Jacobian has at least one positive eigenvalue and one negative eigenvalue.Therefore, each of these corner equilibria is a saddle.
E Interior equilibria of the three-dimensional system with θ 1 = θ 12 In this section, we derive the stability requirements for the line of interior equilibria, L. By setting θ 1 = θ 12 and x * = 1 1+θ 1 , we obtain the (2, 2) and (3, 3) entries of the Jacobian given by Eq. ( 16) as follows: Therefore, we obtain By substituting x * = 1 1+θ 1 , we obtain ∂g ∂n Therefore, we obtain Likewise, using and Now, let us calculate the quantity for x By substituting Eq. ( 43) in Eq. (127), we obtain Next, using x * = 1 1+θ 1 and Eq. ( 43), we find Using Eqs. ( 128) and (129), we obtain Using Eqs. ( 124), (126), and (130), we find that the Jacobian at any point of L is given by where The characteristic equation is given by Eigenvalue λ 1 = 0 reflects the fact that the line of equilibria, L, is neutrally stable along the direction of L. The other two eigenvalues, λ 2 and λ 3 , are given by The real part of λ 2 and λ 3 is negative if and only if α > 0 and β > 0. Because the denominator of Eq. ( 130) is positive, then α > 0 if and only if Eq. ( 19) holds true.Now we seek the conditions under which β > 0. Because n * 1 and n * 12 are positive, we obtain σ 4 > 0 and σ 5 > 0. Therefore, a sufficient condition for β > 0 is that both σ 2 and σ 3 are negative.Using the assumptions in Eq. ( 2), we find that the numerators of σ 2 and σ 3 are always negative.Thus, under Eqs.( 2) and (19), we obtain α > 0 and β > 0 such that the real parts of λ 2 and λ 3 are negative.
For positive δ values satisfying Eq. ( 145) to exist, it must hold true that R 0 − T 0 > S 0 − P 0 .
Eigenvalues J 22 and J 33 are non-positive if and only if and respectively, which implies that Lastly, eigenvalue J 11 is negative if either δ < min holds true.However, we find that there is no δ value that simultaneously satisfies Eqs. ( 145), (148), and (149), or one that simultaneously satisfies Eqs.(145), (148), and (150).
Case 2: γ > 0 and µ 3 < 0 If γ > 0, then µ 3 < 0 if and only if For positive δ values satisfying Eq. ( 151) to exist and be less than 1, it must hold true that R 1 − T 1 < S 1 − P 1 .Eigenvalues J 22 and J 33 are non-positive if and only if and respectively, which implies that Eigenvalues J 22 and J 33 are non-positive if and only if respectively, which implies that Lastly, J 11 is negative if either holds true.However, we find that there is no δ value that simultaneously satisfies Eqs. ( 168), (171), and (172) or Eqs.(168), (171), and (173).
Case 2: γ > 0 and µ 4 < 0 If γ is positive, then µ 4 < 0 if and only if For positive δ values satisfying Eq. ( 174) to exist, it must hold true that R 1 − T 1 < S 1 − P 1 .Eigenvalues J 22 and J 33 are non-positive if and only if and respectively, which implies that Because we have assumed that µ 4 < 0, the numerator of Eq. (162) has to be positive for eigenvalue J and We find that there is no δ value that satisfies Eqs. ( 174), ( 177), (178), and (179).
Case 4: γ < 0 and µ 4 < 0 In this section, we assume that γ < 0 and µ 4 < 0, which requires Eq. (168).For positive δ values satisfying Eq. (168) to exist and be less than 1, it must hold true that R 0 −T 0 > S 0 −P 0 .The derivation of the stability of this case follows the same derivation as Case 2, and we find that the equilibrium is marginally stable if and only if Eqs. ( 19) and ( 45) hold true.

G Equilibria of the five-dimensional system
We show the 60 equilibria of the five-dimensional system with their stability requirements in Table 2.
Table 2: Equilbria of the five-dimensional system with their stability requirements.Symbol "a" represents an equilibrium value between 0 and 1.When the analytical expression is too complicated, we show the numerical values to the third significant digit.

Figure 4 :
Figure 4: Visualization of the transcritical bifurcations as δ varies.We use the payoff matrices given by Eq. (38).The solid and dashed lines indicate stable and unstable equilibria, respectively, both disregarding the 0 eigenvalues along the direction of L in the case of θ 1 = θ 12 .(a) Movement of three equilibria in the full state space as δ varies when θ 1 = 5 and θ 12 = 8.A transcritical bifurcation occurs involving the face equilibrium on n 12 = 1 and the edge equilibrium (x * , 0, 1), where x * = 1/6, at δ = 7/22.The second transcritical bifurcation occurs involving the face equilibrium on n 1 = 0 and the edge equilibirium (x * , 0, 1), where x * = 1/9, at δ = 10/31.(b) Positions of all the same three edge and face equilibria as a function of δ.The θ 1 and θ 12 values are the same as those used in (a).In (b), the three curves do not meet at a single point, as shown in the inset, which is a magnification of the main panel.(c) Same as (a) but when θ 1 = θ 12 = 5.A transcritical bifurcation occurs involving the face equilibrium on n 12 = 1 and that on n 1 = 0 at 1 6 , 0, 1 when δ = 7/22.Edge equilibrium (x * , 0, 1) also collides with the two face equilibria at this value of δ.(d) Same as (b) but when θ 1 = θ 12 = 5.There is another triplet of equilibria in addition to the triplet of equilibria shown in (c).For this second set of triplet of equilibria, a transcritical bifurcation occurs involving the face equilibrium on n 12 = 0 and that on n 1 = 1, and edge equilibrium (x * , 1, 0) collides with the bifurcation point, at δ = 15/22.Note that x * is not constant along the trajectories in (b), whereas it is in (d).