Coevolutionary dynamics of feedback-evolving games in structured populations

When either of the conditions is satisfied: $T-R>{\rm max}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$ or $T-R<{\rm min}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$, the system has five equilibria and exhibits two distinct dynamics. These dynamics are primarily determined by whether the interior equilibrium is a center or a saddle point. The specific conclusions are provided in the following theorems.

The detailed proof of Theorem 1 can be found in A.2. This theorem indicates that if the condition $T-R>{\rm max}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$ is satisfied, the frequency of cooperation and the environmental state will exhibit persistent oscillations over time, regardless of the system’s initial state. This indicates that the system will consistently exhibit self-recovery capabilities.

In the following, a proposition concerning the existence of a heteroclinic cycle is presented under the conditions of Theorem 1.

The proof of Proposition 1 is provided in A.2. Based on Theorem 1 and Proposition 1, it can be observed that all solutions of the system (2) located inside the heteroclinic cycle are periodic, excluding the interior equilibrium $\left(\frac{1}{1+\theta},\frac{1}{2}\right)$.

However, periodic solutions may not always exist when there are five equilibria. When $T-R<{\rm min}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$, a phenomenon occurs in which the cooperation or defection disappears, while the environment tends towards a state of saturation or depletion, respectively. The corresponding conclusions are provided in the theorem below.

The proof details of Theorem 2 are given in A.2. When $T-R<{\rm min}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$, the system (2) exhibits bistability, with two stable equilibria represented by $(0,0)$ and $(1,1)$, respectively. Since the saddle point $(\frac{1}{1+\theta},\frac{1}{2})$ is the sole interior equilibrium, the basins of attraction of $(0,0)$ and $(1,1)$ are bounded by the stable manifold of $(\frac{1}{1+\theta},\frac{1}{2})$. Accordingly, the global dynamics correspond to that for any initial condition $(x_{0},m_{0})$ where $x_{0}\in(0,1)$ and $m_{0}\in(0,1)$ and the point does not lie on the stable manifolds of the saddle point, the system will inevitably converge to one of two states: either cooperators dominate in the most abundant environmental state, or defectors dominate in a depleted environmental state. Moreover, a higher frequency of cooperation or a more abundant environmental state is conducive to the system tending toward a more favorable stable state.

In the next subsection, we will present the results of the stability analysis for cases where the system has seven equilibrium points.

When $1-k<\frac{T-R}{P-S}<\frac{1}{1-k}$, the system (2) has seven equilibria, resulting in more complicated dynamical behaviors. Compared to the case with five equilibria, there are two more equilibria, $(x^{*},0)$ and $(x^{*},1)$. They represent the coexistence of cooperators and defectors in the depleted and replete environment, respectively. Correspondingly, the specific conclusions about the system dynamics are given by the following theorems.

The proof of Theorem 3 is provided in A.4. This theorem indicates that the system (2) exhibits a coexistence of persistent oscillations and the domination of cooperation or defection when the inequality $\frac{\theta k-\theta+1}{1-\theta-k}(P-S)<T-R<{\rm max}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$ is satisfied. More specifically, the interior equilibrium point $\left(\frac{1}{1+\theta},\frac{1}{2}\right)$ is a center, similar to the result in Theorem 1, and is surrounded by periodic solutions. Additionally, the system has a stable equilibrium at either $(0,0)$, representing that defectors dominate in a depleted environment, or $(1,1)$, representing that cooperators dominate in a replete environment.

Similar to Proposition 1, a heteroclinic cycle can also be obtained under the conditions of Theorem 3.

The proof of Proposition 2 is provided in A.4. Indeed, all trajectories starting from $(x_{0},m_{0})$ within the interior of the heteroclinic cycle, but not at the center, form closed periodic orbits. In other words, inside the heterocyclic cycle, both the frequency of cooperation and the environmental state display persistent oscillations. By combining Theorem 3 and Proposition 2, it is obtained that when $T-R>0$, $P-S<0$, and $1-k<\frac{T-R}{P-S}<\frac{\theta k-\theta+1}{1-\theta-k}$, the system initialized with a low frequency of cooperators will inevitably reach a state where defectors dominate in the depleted environmental state, even if its initial environment state is favorable. Conversely, when $T-R<0$, $P-S>0$, and $\frac{\theta k-\theta+1}{1-\theta-k}<\frac{T-R}{P-S}<\frac{1}{1-k}$, the system initialized with a high frequency of cooperators will converge to a state where cooperators dominate in the replete environment, regardless of the initial environmental state. Furthermore, as indicated by Proposition 1 and 2, a heteroclinic cycle exists if and only if the system’s interior equilibrium point $\left(\frac{1}{1+\theta},\frac{1}{2}\right)$ is a center.

It is worth noting that due to the existence of a heteroclinic cycle, it is clear to observe the global dynamics of the system with oscillatory behavior and easier to analyze the impact of population structure on the system dynamics. For instance, when $T-R>0$ and $P-S>0$, the frequency of cooperators and the environmental state oscillate cyclically and indefinitely, regardless of the number of neighbors $k$. However, when $\frac{\theta k-\theta+1}{1-\theta-k}(P-S)<T-R<{\rm max}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$, the system exhibits two distinct behaviors: oscillations within the heteroclinic cycle and dominance of cooperation or defection outside of it. Therefore, the size of the regions where these two dynamics occur separately is influenced by the number of neighbors. Specifically, the number of neighbors impacts the value of $x^{*}$, thereby determining the size of the heteroclinic cycle.

We also find an alternative dynamic of the system (2) that includes seven equilibrium points as described in the following theorem.

The proof details of Theorem 4 are given in A.5. This theorem suggests that when the inequality ${\rm min}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}<T-R<\frac{\theta k-\theta+1}{1-\theta-k}(P-S)$ is satisfied, there exist two stable equilibria: a stable coexistence of cooperation and defection, and a dominant of defection or cooperation. There are two stable manifolds converging to the interior saddle point, each originating from a separate unstable node. When $T-R>0$, $P-S<0$, and $\frac{\theta k-\theta+1}{1-\theta-k}<\frac{T-R}{P-S}<\frac{1}{1-k}$, the coexisting stable equilibrium is $(x^{*},1)$, and the stable dominance equilibrium is $(0,0)$. If the initial condition $(x_{0},m_{0})$ lies below the two stable manifolds originating from $(0,1)$ and $(x^{*},0)$, the solution will converge to $(0,0)$. Otherwise, it will converge to $(x^{*},1)$. In other words, if the initial cooperative frequency is low or the environmental state is poor, the system will evolve to a state where defectors dominate in a depleted environment. Otherwise, cooperators and defectors will coexist, and the environmental state will remain replete. When $T-R<0$, $P-S>0$, and $1-k<\frac{T-R}{P-S}<\frac{\theta k-\theta+1}{1-\theta-k}$, the stable coexisting equilibrium of the system is $(x^{*},0)$, and the stable dominance equilibrium is $(1,1)$. In this case, if the initial cooperative frequency is high or the environmental state is rich, the system will eventually evolve into a state where cooperators dominate in the replete environmental state.

In summary, the coevolutionary dynamics of feedback-evolving games in structured populations are more complicated compared to those in well-mixed populations (i.e., in the limit case of $k\to\infty$). This fully reflects the significant impact of the number of neighbors on the coevolutionary dynamics. Moreover, according to the theoretical results obtained above, it is found that, $\theta$, a parameter representing the ratio of enhancement rates to degradation rates for cooperators and defectors, does not alter the diversity of dynamics in the system.

If not specified, set $\epsilon=0.1$, $\theta=1.8$ and $k=3$. Under these conditions, the interior equilibrium of the system (2) is $(\frac{5}{14},\frac{1}{2})$.

To verify the results of Theorem 1 and Proposition 1, we provide a numerical example in Fig. 1. For the example, we assume $R=3$, $S=0$, $T=5$, and $P=1$, which satisfy the conditions $T-R>0$ and $P-S>0$.

In Fig. 1(a), the vector fields on the phase plane are depicted. We observe that trajectories along the boundary of the region form a neutral heteroclinic cycle. Additionally, any trajectories within this cycle, except at the point $(\frac{5}{14},\frac{1}{2})$, represents a closed orbit. Furthermore, Fig. 1(b) shows a time series plot about the fraction of cooperators and the environmental state. We further verify that the solution, starting from the initial condition, completes one full counterclockwise revolution and returns to its starting point. This indicates that the frequency of cooperators and the environmental state exhibit persistent oscillations between 0 and 1, which further verifies Theorem 1.

To verify the result of Theorem 2, we provide an example in Fig. 2. Assume that $R=3$, $S=0$, $T=5$, and $P=1$, which satisfy the conditions $T-R<0$ and $P-S<0$.

In Fig. 2(a), the phase plane is depicted with two green curves. One of the curves connects the points $(0,1)$ and $(1,0)$, representing two stable manifolds of the saddle point $(\frac{5}{14},\frac{1}{2})$, with the vector field directed toward the saddle point. The other green curve connects the points $(0,0)$ and $(1,1)$, representing two unstable manifolds, with the vector field pointing away from the saddle point. Clearly, solutions above the stable manifolds eventually converge towards the stable equilibrium $(1,1)$, representing a state where cooperators dominate in the most abundant environmental state. Conversely, solutions below the unstable manifolds converge to the stable equilibrium $(0,0)$, indicating a state where defectors dominate in the most scarce environmental state. The time series plots in Fig. 2(b) and (c) confirm the evolutionary outcomes of the system, where either cooperators dominate in a favorable environmental state or defectors prevail in an unfavorable environmental state. Therefore the result of Theorem 2 is verified.

Here, we will verify the results of the coevolutionary dynamics when the system has seven equilibria. As usual, set $\epsilon=0.1$, $\theta=1.8$, and $k=3$. Then, the interior equilibrium of the system (2) is $(\frac{5}{14},\frac{1}{2})$.

To verify Theorem 3, we consider the parameter values as $R=3$, $S=1$, $T=4.9$, and $P=0$, which satisfy $T-R>0$, $P-S<0$, and $1-k<\frac{T-R}{P-S}<\frac{\theta k-\theta+1}{1-\theta-k}$. As shown in Fig. 3(a)-(c), the system exhibits a single stable equilibrium at $(0,0)$, corresponding to a state where defectors dominate in a depleted environment. Next, we set $R=3.6$, $S=0$, $T=3$, and $P=1$, which satisfy $T-R<0$, $P-S>0$, and $\frac{\theta k-\theta+1}{1-\theta-k}<\frac{T-R}{P-S}<\frac{1}{1-k}$. As depicted in Fig. 3(d)-(f), the system has a single stable equilibrium $(1,1)$, representing a state where cooperators dominate in a replete environmental state. In both cases, if the initial conditions lie within a heteroclinic cycle, the frequency of cooperators and the environmental state exhibit persistent oscillations. The heteroclinic cycle is characterized by the trajectories $(1,0)\to(1,1)\to(\frac{1}{29},1)\to(\frac{1}{29},0)\to(1,0)$ and $(0,0)\to(\frac{7}{8},0)\to(\frac{7}{8},1)\to(0,1)\to(0,0)$, respectively. Therefore, the system demonstrates a coexistence of oscillation and dominance when $\frac{\theta k-\theta+1}{1-\theta-k}(P-S)<T-R<{\rm max}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}$. This verifies the result of Theorem 3.

To validate Theorem 4, we consider the parameter values as $R=3$, $S=1$, $T=3.6$, and $P=0$, which satisfy $T-R>0$, $P-S<0$, and $\frac{\theta k-\theta+1}{1-\theta-k}<\frac{T-R}{P-S}<\frac{1}{1-k}$. As shown in the top row of Fig. 4, the system has two stable equilibria, $(0,0)$ and $(\frac{7}{8},1)$. This indicates that, regardless of the initial conditions, the environmental state will evolve toward one of two extreme states: depletion or repletion. Furthermore, when the system reaches a replete environmental state, the majority of individuals tend to cooperate. Next, we set $R=4.9$, $S=0$, $T=3$, and $P=1$, which satisfy $T-R<0$, $P-S>0$, and $1-k<\frac{T-R}{P-S}<\frac{\theta k-\theta+1}{1-\theta-k}$. As shown in the bottom row of Fig. 4, the system also has two stable equilibria, $(\frac{1}{29},0)$ and $(1,1)$. In this case, a similar result is observed: the system evolves towards either a depleted or replete environmental state, with full cooperation occurring in the replete environmental state. Thus, the system exhibits a coexistence of cooperation and defection when ${\rm min}\left\{\frac{P-S}{1-k},(1-k)(P-S)\right\}<T-R<\frac{\theta k-\theta+1}{1-\theta-k}(P-S)$. This confirms Theorem 4.

In order to help readers easily overview the evolutionary outcomes of our dynamical system for different parameter regions, we provide several schematic plots to elucidate the diversity of coevolutionary dynamics in structured populations. As shown in Fig. 5, different colors correspond to the various dynamics previously described.

Specifically, gray regions indicate oscillations, green regions indicate bistability, red regions represent the coexistence of oscillation and domination, and blue regions mean the coexistence of cooperation and defection. As the degree $k$ approaches infinity, the population tends to become well-mixed, exhibiting only two dynamical behaviors when $\frac{T-R}{P-S}<0$: the coexistence of oscillation and domination, and the coexistence of cooperation and defection. This reveals that the dynamics in a structured population are more diverse and complicated than those in a well-mixed population, and this disparity on dynamics becomes more pronounced as the parameter $k$ decreases. Furthermore, comparing the first and second columns of Fig. 5, we find that the value of $\theta$ does not impact the diversity of the dynamics, but it can influence their distribution.

In addition, we can learn that the population structure plays a crucial role in shaping coevolutionary dynamics. In particular, the number of neighbors is a key parameter in the spatially structured population. Subsequently, we will investigate how the neighbor number influences the coevolutionary dynamics, particularly its effect on the evolution of the environmental state.

Note that the significant impact of neighborhood size $k$ on the system dynamics with fixed payoff parameters is complicated due to the diverse dynamics involved. Here, our primary focus is on the dynamics featuring two stable equilibria: bistability and the coexistence of cooperation and defection. Specifically, we analyze the basin of attraction of the equilibrium characterized by the replete environmental state and its dependence on the number of neighbors. This analysis is crucial for understanding the maintenance of environmental state in a population. The outcomes shown in Fig. 6 and Fig. 7 are obtained through Monte Carlo simulations. To reach the expected accuracy for each value of the parameter $k$, we generate $10^{6}$ random points in the whole regions of $0\leq x\leq 1$ and $0\leq m\leq 1$ as initial states. The state trajectory starting from each initial state point is then computed by using numerical calculations and the system along each trajectory will converge to a specific equilibrium, either $(1,1)$ or $(x^{*},1)$, where $x^{*}$ depends on the given value of parameter $k$. We finally compute the fraction of the trajectories which eventuate in the $(1,1)$ state as the proportion of attraction region of the stable state $(1,1)$, and the fraction of the trajectories which eventuate in the $(x^{*},1)$ state is evaluated as the ratio of attraction region of the stable state $(x^{*},1)$.

When the system exhibits bistability (corresponding to the case of $T-R<0$ and $P-S<0$), the impact of $k$ on dynamics is closely related to $\theta$. There is an ideal stable state, $(1,1)$, where cooperators dominate in the replete environmental state. As depicted in Fig. 6, for higher values of $\theta$, a larger number of neighbors is more likely to be beneficial for environmental development. Conversely, a smaller number of neighbors tends to be more advantageous when the value of $\theta$ is lower.

When the system exhibits coexistence of cooperation and defection, we focus only on situations where the number and type of equilibria remain unchanged as $k$ varies. In this case, a favorable stable equilibrium, $(x^{*},1)$ or $(1,1)$, is characterized by a replete environmental state. Besides, the size of the attraction basin for this stable equilibrium is closely related to $x^{*}$, which is independent of $\epsilon$. For simplicity, we assume $\epsilon=0.1$ without loss of generality. As shown in Fig. 7, we find that the role of $\theta$ in the correlation between the size of the attraction basin and the parameter $k$ can be negligible. Instead, the value of $\frac{T-R}{P-S}$, representing the ratio of temptation minus reward to punishment minus sucker’s payoff, plays a pivotal role. Specifically, when the value of $\frac{T-R}{P-S}$ is relatively high, the attraction basin of the stable equilibrium $(x^{*},1)$ or $(1,1)$ expands as $k$ increases (see Fig. 7(a) and (b)). Conversely, when the value of $\frac{T-R}{P-S}$ is relatively low, the attraction basin shrinks as $k$ increases (see Fig. 7(c) and (d)). In other words, when the ratio of temptation minus reward to punishment minus sucker’s payoff is closer to zero (but still negative), it is beneficial to have more neighbors. However, when the ratio of temptation minus reward to punishment minus sucker’s payoff is more negative, it is advantageous to have fewer neighbors.