{'fraction_non_alphanumeric': 0.05248539583557323, 'fraction_numerical': 0.016145217360140487, 'mean_word_length': 4.100255894717602, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Limit Cycles Analysis and Control of Evolutionary Game Dynamics with Environmental Feedback

Recently, an evolutionary game dynamics model taking into account the environmental feedback has been proposed to describe the co-evolution of strategic actions of a population of individuals and the state of the surrounding environment; correspondingly a range of interesting dynamic behaviors have been reported. In this paper, we provide new theoretical insight into such behaviors and discuss control options. Instead of the standard replicator dynamics, we use a more realistic and comprehensive model of replicator-mutator dynamics, to describe the strategic evolution of the population. After integrating the environment feedback, we study the effect of mutations on the resulting closed-loop system dynamics. We prove the conditions for two types of bifurcations, Hopf bifurcation and Heteroclinic bifurcation, both of which result in stable limit cycles. These limit cycles have not been identified in existing works, and we further prove that such limit cycles are in fact persistent in a large parameter space and are almost globally stable. In the end, an intuitive control policy based on incentives is applied, and the effectiveness of this control policy is examined by analysis and simulations.


Introduction
Evolutionary game theory studies the evolution of strategic decision-making processes of individuals in one or multiple populations. A range of mathematical models have been extensively investigated, among which the replicator dynamics model is the most wellknown [10]. Powerful theories and tools originating from evolutionary game theory facilitate the study of complex system behaviors of biological, ecological, social, and engineering fields [19], [12], [14], [22], [30]. In the classic game setting, the payoffs in each two-player game are usually predetermined in the form of constant payoff matrices. However, in many applications, especially those involving common resources in the environment that are consumed by groups of individuals, it is recognized that the payoffs can change over time or be affected directly by external environmental factors. Thus game-playing individuals' decisions can influence the surrounding environment, and the environment, especially its richness in resources, also acts back on the payoff distributions, which forms the socalled bi-directional game-environment feedback [35]. Such a feedback mechanism is believed to be able to explain rich real-world complex population dynamics, including human decision-making, social culture evolution, plant nutrient acquisition, and natural resource harvesting. In particular, evolutionary dynamics with game-environment feedback have attracted great attention in recent years, due to its high relevance in biological and sociological systems [24], [14], [32].
This line of research is originated in [35], where the replicator dynamics of a two-player two-strategy game are coupled with the logistic environment resource dynamics. Utilizing an elaborate environment-dependent payoff matrix, the joint game-environment dynamics exhibit interesting system behaviors. In most cases, the strategic states of the population and the environment converge to an equilibrium point on the boundary of the phase space with the zero value of environment state, which reflects the tragedy of the commons. The system dynamics may also show cyclic oscillations which correspond to closed periodic orbits. Particularly when the system dynamics converge to a heteroclinic cycle on the boundary, the environment state cycles between low and high values, which is referred to as the oscillating tragedy of commons. Since then, [35] has been followed by several extensions and variations [7], [32], [4], [11], [18], [15]. The game-environment feedback mechanism has been studied more broadly in different applications, while its theoretical settings have been adapted in different research directions, e.g. implementing different types of environment resource models and considering the structures of the populations. New meaningful dynamical behaviors have been identified, and these results are in turn helpful to understand real-world applications.
The majority of the existing studies on the gameenvironment framework chooses replicator dynamics to model the strategic dynamics of the game-environment system. While the replicator dynamics have been proved to be a powerful model in analyzing a variety of classical games, they do not take into account the mutation or exploration in strategies which exists ubiquitously in biological and sociological systems. Genetic mutations typically occur with small but non-negligible probabilities, and random exploration of available strategies is common in behavioral experiments [33]. Mutations and explorations can be captured by allowing individuals to spontaneously switch from one strategy to another in small probability, which results in the so-called replicator-mutator dynamics. It has also played a prominent role in evolutionary game theory and has been applied frequently in different fields [19], [12], [23], [22].
We note that although rich system behaviors, such as convergence to equilibria or heteroclinic cycles on the boundary and existence of neutral periodic orbits, have been identified in related works, the results on the limit cycle dynamics have been rarely reported in the gameenvironment feedback setting. Only in [32], it is shown that the time-scale difference between game and environment dynamics can result in limit cycles. When mutations are taken into account, it is of great interest to study if the coupled game and environment dynamics can exhibit limit cycles.
We also emphasize that research on how to design reasonable and effective control policies to achieve expected system behaviors in evolutionary games has attracted increasing attention [17], [37], [25], [27]. Specifically, using game-environment feedback, [21] considers optimal control problems where control takes the form of the incentives and opinions respectively. The adopted control policies can maximize accumulated resource over time. However, since the controlled system still exhibits heteroclinic oscillations, repeated collapses of the resource are inevitable. [34] considers adjusting the law of feedback from population states to the environment based on a slightly different co-evolutionary model of public goods games. An involved nonlinear control law is proved to be able to steer the system to evolve towards the designed behaviors, but the biological implication of such control policies still remains open to debate.
Motivated by all the unsolved issues just mentioned, in this work we study the simultaneous co-evolution of game and environment using an integrated model consisting of replicator-mutator dynamics and logistic resource dynamics. Although similar linear environmentdependent payoff matrices have been studied before, we focus on investigating the significant impact of mutations on the overall system dynamics using bifurcation theory as the main tool. In particular, we show that two different types of bifurcations, which lead to limit cycles, can take place. In sharp comparison to the marginal oscillating behaviors reported before, such as neutral periodic orbits and heteroclinic cycles, we clarify that the identified limit cycle conditions correspond to a large space of the mutation parameter and such limit cycles are almost globally stable. In this sense, the limit cycle behavior is robust and persistent, which agrees with biological observations. Furthermore, when mutations are not too rare, the corresponding limit cycle will not be close to the boundary, and consequently the (oscillating) tragedy of commons can be averted.
Enabled by the new insight into co-evolutionary dynamics of the game and environment, we consider further a simple but effective incentive-based control policy, where the individuals receive some external incentives in the game interactions if they choose the cooperation strategy. It is shown that the controlled system dynamics can converge to certain equilibria with a high level of the environmental resource. Thus, the effectiveness of the proposed control policy is guaranteed. Compared with the proposed control policies in [ [21]] and [ [34]], our proposed approach is not only more effective, but also easier to implement than influencing the opinions of the individuals or realizing more complicated feedback between strategists and the environment.
Our preliminary conference paper [5] considers Hopf bifurcation and limit cycles of the system dynamics in a special case where the enhancement effect and the degradation effect is balanced, but in this paper we consider general bifurcation conditions and provide a comprehensive account of the results. Moreover, the latter part about incentive-based control of this paper is not in the conference version.
The remainder of this paper is organized as follows. Section 2 provides the framework of the co-evolutionary dynamics, and introduces the integrated model with gameenvironment feedback. The analysis of system dynamics is presented in Section 3. The incentive-based control is applied and its performance is provided in Section 4. Finally, concluding remarks and discussion points are given in Section 5.

Problem Formulation
Co-evolutionary game and environment dynamics occur because evolutionary game dynamics are sometimes coupled with the dynamics of the surrounding environment (see Fig. 1 for an illustrative drawing). The coupled game and environment dynamics form a closed-loop system, which consists of bi-directional actions between the environment and the population's game play: an individual's fitness depends on not only the population state, but also the state of the environment, and the state of the environment is influenced by the distribution of the choices of strategies in the population. In the following subsections, we will introduce each component of this system in turn, and then integrate them to give a complete description of our model.

Environment-dependent payoffs
We consider games played in a changing environment, which is characterized by the richness r ∈ [0, 1] of a resource of interest; such a resource has direct influence on the payoffs that players receive. Correspondingly, the payoffs become dynamic and depend on r. We consider a two-player game with two strategies, Cooperation and Defection (C, D). Following [35], we assume that the environment-dependent payoff matrix is a convex combination of the payoff matrices from the classic Prisoner's Dilemma, i.e., where R 1 > T 1 , S 1 > P 1 , R 2 < T 2 and S 2 < P 2 . The entries of A(r) represent the payoffs to the players at the given environment state r. For example, a C player receives the payoff of R 1 + (R 2 − R 1 )r and S 1 + (S 2 − S 1 )r when the opponent uses strategies C and D, respectively. Because of the inequalities below (1), the two matrices R 2 S 2 T 2 P 2 and R 1 S 1 T 1 P 1 correspond to the prisoner's dilemma and its inverse game, respectively. Note that the payoff matrix A(r) is the linear interpolation between the above two matrices. Thus, when r approaches 0 and 1, the mutual cooperation and defection tend to be the Nash Equilibrium of this dynamic game respectively. The biological interpretation is that: the individuals tend to cooperate facing the depleted environmental resource, and become more selfish and choose to defect while the resource becomes abundant.

Replicator-Mutator equations
The two-player game governed by the payoff matrix A(r) in (1)  We focus on the dynamics of x, which can be described by the following replicator-mutator equatioṅ where (A(r)x) 1 is the first entry of A(r)x and represents the fitness of choosing strategy C, and x A(r)x is the average fitness at the population state x. The term −µx represents the mutation from strategy C towards D, while µ(1 − x) captures the mutation in the opposite direction, from D towards C.

Closed-loop system with game-environment feedback
To model the evolution of the environmental change, we use the standard logistic model where the parameter θ > 0 represents the ratio between the enhancement effect due to cooperation and degradation effect due to defection. The influence of the population on the environment state is captured by the function [θx − (1 − x)], and the dynamics are restricted in the unit interval by r(1 − r). When θ = 1, the enhancement effect and the degradation effect are balanced. And if θ > 1 or 0 < θ < 1, the enhancement effect is respectively stronger or weaker than the opposite.
Substituting (1) into (2), and then combining with (3), we obtain the closed-loop planar system describing the evolutionary game dynamics under the gameenvironment feedback: where to simplify notations, we have used a = R 1 − T 1 , b = S 1 − P 1 , c = T 2 − R 2 , and d = P 2 − S 2 , all of which are positive because of the inequalities below (1). Note that system (4) reduces to the model considered in [35] when µ = 0. Therefore, system (4) generalizes the previous model.
In view of the payoff matrix (1) and system equations (4), these four parameters for re-parameterization have intuitive interpretations with respect to strategy changes when the resource is either abundant or depleted: the parameter a quantifies the incentive to stick to strategy C given that all individuals are following strategy C and the system is currently in a resource-poor state; b quantifies the incentive to switch to strategy C when the resource is depleted and all individuals enforce strategy D; in contrast, c quantifies the incentive to switch to strategy D given that all individuals enforce strategy C and the system is in a state of sufficient resource; d quantifies the incentive to stick to strategy D when the resource is abundant and all individuals are following strategy D.

System Dynamics and Bifurcations
The state space of system (4) is the unit square I = [0, 1] 2 . The interior region of this square is denoted by int(I) = (0, 1) 2 . Four sides of the unit square form the boundary ∂I. For convenience of exposition, we denote the four sides by: Then, one has ∂I := B t ∪ B b ∪ B l ∪ B r . It is easy to check that I is invariant under the dynamics (4) for all µ ∈ [0, 1] by Nagumo's theorem [1,Theorem 4.7]. System (4) can have equilibria in the interior int(I); we call them interior equilibria. System (4) also has equilibria on the boundary ∂I (boundary equilibria), which will be discussed in detail in Section 3.4. It is easy to check that the interior equilibria, if they exist, must lie on the line {(x, r) ∈ int(I) : x = 1/(θ + 1)}. Substituting such x into the x-dynamics in (4), we have which leads to Thus if the following condition holds system (4) has a unique interior equilibrium θ(a + c + θb + θd) .
(7) Otherwise, no interior equilibrium exists. The interior equilibrium supports the coexistence of the two strategies under a non-depleted resource state. In view of (6), one can see that for a given θ, the existence of the interior equilibrium is determined by not only the payoff parameters but also the mutation rate.
Note that when θ = 1, (x * , r * ) is independent of the parameter µ, and (6) holds trivially for all µ ∈ [0, 1]. This specific case has been studied in our previous conference paper [5]. In this paper we do not require θ = 1 and study the generic θ that takes any positive value satisfying (6).
Towards the goal of establishing conditions for the existence of limit cycles, we first present some preliminary results on the non-existence of closed orbits or limit cycles.

Non-existence of closed orbits
First, we define the C 1 function ϕ(x, r) : int(I) → R by with the parameters α, β, γ, and δ to be determined later. Note that ϕ(x, r) is strictly positive in int(I) for any values of these parameters. Denote the right hand side of (4) by f (x, r) = [f 1 (x, r), f 2 (x, r)] , and then the divergence of the modified vector field ϕf (x, r) on int(I) is given by Let all the coefficients of the powers of x, r in (9) (except for the term containing µ) be zero, and we obtain the following set of equations We first show that the equation set (10) always has solutions.

PROOF. See Appendix A. 2
In fact, if one of −c + d − a + b and a − b is not zero, (10) will have exactly one solution with ζ = a + c + θb + θd > 0. Otherwise, (10) will have infinitely many solutions including (11). So we fix the parameters α, β, γ, and δ at the values given in (11) such that the divergence (9) can be greatly simplified.
Then we show the non-existence of closed orbits for system (4).

Lemma 2
For system (4), the following statements hold: (1) when ad − bc > 0, there exists some µ 0 > 0 such that there are no closed orbits in I for µ ∈ [µ 0 , 1]; (2) when ad − bc < 0, there are no closed orbits in I for µ ∈ [0, 1]; (3) when ad − bc = 0, there are no closed orbits in I for µ ∈ (0, 1]. PROOF. The proof is direct by applying Bendixson-Dulac criterion [36, Theorem 4.1.2]. Since the boundary ∂I contains equilibria (as shown in Section 3.4), it is not possible to have periodic orbits on it. And the sides are either invariant or repelling, so they cannot intersect with closed orbits. Now let us consider the interior int(I) which is simply connected.
For the second case, it is easy to prove since the right hand side of (14) will always be negative for any µ ≥ 0. For the third case, the right hand side of (14) is always negative when µ > 0. 2 In Lemma 2, it is noted that the conditions for the nonexistence of periodic orbits depend on the relationship between ad and bc. The reason is as follows: when ad − bc < 0, the divergence div ϕf (x, r) will remain positive for any µ ∈ [0, 1] in view of (14), which excludes the possibility of periodic orbits in the phase space; when ad − bc > 0 or ad − bc = 0, there exists some µ such that div ϕf (x, r) equals 0. As a consequence, it is possible for the system to have periodic orbits in these two cases.
Another point worth noting is that the only difference between the second and third cases of Lemma 2 is that system (4) may have periodic solutions for µ = 0 when ad − bc = 0. This possibility has been shown in [35], [4]. Despite the possibility of closed orbits under certain conditions, for µ = 0 it can be proven that limit cycles (isolated closed orbits) are impossible to exist in (4) for all the parameter space as summarized below in Lemma 3.
Lemma 3 System (4) has no limit cycles in I when µ = 0.
The full proof is presented in [5]. See Appendix B for a sketch of the proof.

Hopf bifurcation at the interior equilibrium
In this section we take µ as the bifurcation parameter to study the effect of different mutation rates on the system dynamics. We investigate the possible bifurcations of the dynamics (4) as µ varies in [0, 1]. We analyze the stability of the interior equilibrium and prove that it involves a Hopf bifurcation which leads to periodic orbits.
We first invoke the Hopf bifurcation theorem [6] which will be used to prove the existence of stable limit cycles.
Then we make an assumption about some parameters of system (4).
The following theorem demonstrates that the interior equilibrium (x * , r * ) of system (4) undergoes a Hopf bifurcation.
Theorem 6 Under Assumption 5, the interior equilibrium (x * , r * ) of system (4) with the bifurcation parameter µ undergoes a supercritical Hopf bifurcation at µ = , which leads to the existence of stable limit cycles for µ < µ 1 in the vicinity of µ 1 .

PROOF.
To prove the existence of Hopf bifurcation, we need to show that the two conditions of Theorem 4 are satisfied in system (4). The Jacobian of the vector field of (4) at (x * , r * ) is given by Note that the off-diagonal entries because of (6). One can calculate the eigenvalues of J * λ(µ) = J 11 ± J 2 11 + 4J 12 J 21 2 = (λ(µ)) ± (λ(µ))i. (17) It is easy to check that when µ = µ 1 := θ 2 (ad−bc) ∆ the eigenvalues are purely imaginary, i.e., (λ(µ 1 )) = 0, which implies that the first condition listed in Theorem 4 is satisfied. In addition, the second condition is also satisfied because Then one can conclude that the dynamics (4) undergo a Hopf bifurcation at µ = µ 1 . Since µ varies in the interval (0, 1], to guarantee the occurrence of Hopf bifurcation, the critical parameter value µ 1 should also be constrained in this interval, i.e., 0 < µ 1 ≤ 1, which leads to which is exactly Assumption 5. The Hopf bifurcation theorem implies that a family of periodic orbits bifurcate from (x * , r * ) for some µ in the vicinity of µ 1 . Whether the bifurcated periodic orbits are stable or not is determined by the so-called first Lyapunov coefficient 1 at (x * , r * ) when µ = µ 1 . If 1 < 0, then these periodic orbits are stable limit cycles; if 1 > 0, the periodic orbits are repelling. Following the computation procedure in Appendix C, we obtain 1 (µ 1 ) as below which is negative because of (6) and Assumption 5. Thus, the bifurcated periodic orbits are stable limit cycles.
Furthermore, in view of (17), for µ < µ 1 , (x * , r * ) is unstable since the eigenvalues have positive real parts. On the other hand, the interior equilibrium is asymptotically stable for µ > µ 1 with the negative eigenvalues. Hence, the stable limit cycles exist for µ < µ 1 in the vicinity of µ 1 . Since the bifurcation is associated with stable limit cycles, it is a supercritical Hopf bifurcation. 2 The relationship between ad and bc has played an important role in Theorem 6 in the sense that Hopf bifurcation can only happen when ad − bc > 0. In view of the Jacobian (16), when ad − bc ≤ 0, J 11 will remain nonpositive no matter what value µ takes. Thus, as µ varies in [0, 1], the stability of the interior equilibrium does not change, which excludes the existence of a Hopf bifurcation. On the other hand, when ad − bc > 0, to ensure there is a Hopf bifurcation for µ ∈ [0, 1], Assumption 5 has to be satisfied as proved.
From the Hopf bifurcation theory, we know that the limit cycles generated from Hopf bifurcation have the amplitude of O( √ µ 1 − µ). Since Hopf bifurcation is a local bifurcation, the obtained results only hold for µ in the vicinity of µ 1 . Next, we will show that system (4) exhibits another type of bifurcation.

Heteroclinic bifurcation
In the case µ = 0, it has been shown in Lemma 3 that the system (4) cannot exhibit limit cycles for any parameters. We show in this section that limit cycles may occur for µ > 0 in the vicinity of 0.
When µ = 0, system (4) has four equilibria on the boundary ∂I. Denote the four equilibria respectively by E 1 = (0, 0), E 2 = (1, 0), E 3 = (1, 1), and E 4 = (0, 1). Then, one can calculate the Jacobians at these points, which are It is easy to see that the eigenvalues of these equilibria are indeed the diagonal entries of the associated Jacobian matrices, and thus all of these equilibria are saddle points. Next we note that the trace of each Jacobian, which is denoted by T i , is A heteroclinic cycle on the boundary denoted by Λ, which connects the four equilibria and has the orienta- can be identified. Its stability is determined by the payoff parameters.
We refer to the reference [35] for the complete proof for the case of the stable heteroclinic cycle. And the opposite case can be proved using the analogous method according to [9,Corollary 2].
When µ > 0, the four corners are no longer equilibria, and the heteroclinic cycle Λ is "broken". Now we are going to show that the system (4) undergoes a heteroclinic bifurcation at µ = 0 when certain condition holds, and a limit cycle with the same stability of the corresponding heteroclinic cycle may be generated from Λ for µ > 0 sufficiently close to 0.

Definition 8
The internal σ-neighborhood of Λ is the set of all points in I that are at a distance less than σ > 0 from Λ.
PROOF. The proof is straightforward by applying the theorem of the heteroclinic bifurcation [ In addition, all the traces T i of the Jacobians of the equilibria are negative, while their product exist sufficiently small σ, > 0 such that for 0 < µ < , system (4) has at most one limit cycle in the internal σneighbourhood of Λ, and this limit cycle (if it exists) is stable due to the stability of Λ. 2 When ad < bc, according to the second statement of Lemma 7, the heteroclinic cycle is unstable. It is impossible for limit cycles to be generated from this unstable heteroclinic cycle in view of Lemma 2, which claims that there are no closed orbits in I for any µ ≥ 0 when ad − bc < 0.
Compared with the result of Hopf bifurcation, the conclusion of Heteroclinic bifurcation is relatively weaker since Theorem 9 only states that there may exist a unique limit cycle in some neighborhood of Λ. To study the exact conditions for the existence of the unique limit cycle is beyond this work.
Another difference between the limit cycles generated from Hopf bifurcation and Heteroclinic bifurcation is that the limit cycle arising from Hopf bifurcation is small in size in the phase space, as demonstrated in the simulation results in Fig. 2, while the other can be large. In these simulations the payoff matrices are given by The parameter θ is chosen to be 1, and thus the interior equilibrium is fixed and independent of µ. In view of Theorem 6, the Hopf bifurcation point is µ 1 = 0.1543.

Global dynamics
As stated before, the limit cycle analysis for Hopf and Heteroclinic bifurcations are limited to the vicinity of the bifurcation points since the used methods rely mainly on linearization. To study if the limit cycle persists for µ ∈ (0, µ 1 ), it is necessary to analyze the dynamics of (4) for all such µ.

Boundary equilibria and their stability
Recall that the condition (6), which involves µ, guarantees that there is a unique interior equilibrium. For some given payoff parameters and fixed values of θ, there may be some µ such that there is no interior equilibrium. In this case, the system dynamics are relatively trivial since the trajectories will converge to some limit sets on the boundary according to the Poincaré-Bendixson theorem. In the following part, we concentrate on the situation where an interior equilibrium always exists for all µ ∈ [0, 1]. Then, we have the following assumption.
Assumption 10 We assume hereafter that −(θa + Assumption 10 ensures that there is a unique interior equilibrium point for all µ ∈ [0, 1] in system (4). It is noted that when θ = 1, Assumption 10 is satisfied trivially.
It is easy to check that the boundary ∂I : B t ∪B b ∪B l ∪B r itself and each of its sides are invariant under (4) in the case µ = 0. When 0 < µ ≤ 1, the sides B t and B b remain invariant, while B l and B r do not. As µ varies, all the possible boundary equilibria lie on the sides B t and B b , but the positions of these equilibria change. Let us try to determine the locations of the equilibria on B t and B b .

PROOF. See Appendix D. 2
One can examine the stability of these boundary equilibria by evaluating the corresponding Jacobian matrices.

PROOF.
We only show the case on the side B b . For the case of B t , one can obtain the same result analogously. Denote the equilibrium on the side B b by (x * b , 0). We obtain its Jacobian matrix as below where the symbol stands for − As the matrix is upper triangular, the entry (θ+1)x * b −1 is one of its eigenvalues. According to Lemma 11, one has x * b ∈ (max{1/2, 1/(θ + 1)}, 1). Then the eigenvalue (θ + 1)x * b − 1 will be positive for any θ > 0, which implies that the equilibrium is unstable. 2 By checking the sign ofẋ on the left side B l and the right side B r , one observes that it is always positive and negative respectively. It can be easily verified that the vectors on these two sides point inwards to int(I). The sides B b and B t are positively invariant under system (4), and thus the dynamics on each side of B b and B t are easy to analyze. For example, on the side B b , one hasṙ = 0, x > 0 at x = 0 andṙ = 0,ẋ < 0 at x = 1. As there is only one equilibrium (x * b , 0) on the side B b , one can easily verify thatẋ > 0 for . The similar result can be obtained for the dynamics on B t .
In addition, we can show that the boundary ∂I is repelling for µ ∈ (0, 1].

PROOF. See Appendix E. 2
With Lemma 13 in hand, we are ready to present some global results of the system dynamics in I.
PROOF. Recall the two critical parameter values µ 0 and µ 1 , which appear respectively in Lemma 2 and Theorem 6. One can check that µ 0 ≥ µ 1 when Assumption 5 holds. Again according to Lemma 2, it is impossible for system (4) to have closed orbits in I when µ > µ 0 . In addition, the interior equilibrium is locally asymptotically stable in view of (16). Due to the fact that the boundary ∂I is repelling as shown in Lemma 13, any trajectories starting from int(I) \ {(x * , r * )} cannot converge to ∂I. Then according to the Poincaré-Bendixson theorem, because of the non-existence of closed orbits or other equilibria in int(I), the trajectories with any initial points in int(I) \ {(x * , r * )} will converge to (x * , r * ). Similarly, trajectories starting from the sides B l and B r also converge to (x * , r * ). 2 So the interior equilibrium is almost globally stable in system (4), because almost all trajectories in I (except for those starting from B t ∪ B b ) converge to the stable equilibrium (x * , r * ). We apply similar terminology to the stability of the limit cycle when it is stable and almost all trajectories converge to it.
3.4.2 Uniqueness of limit cycles in the case θ = 1 Assumption 15 We assume that θ = 1 henceforth.
Under Assumption 15, the enhancement effect and the degradation effect are balanced. In this case, the system can be transformed to a generalized Liénard system. However, the existence and uniqueness of limit cycles for µ ∈ (0, µ 1 ] in the case θ = 1 is still open. Lemma 16 Under Assumption 5, when a + c = b + d, the system dynamics in (4) admit at most one limit cycle in int(I) for µ ∈ (0, µ 1 ).
The detailed expressions of the functions involved in the above conditions can be found in Appendix F. It is not difficult to check these conditions, so the process is omitted here.
Then according to [28, Theorem 1, Corrollary 1], system (24) has at most one limit cycle inĨ. Therefore, the original system (22) admits at most one limit cycle in int(I). 2 It is noted that when a + c = b + d, the coefficient of the cross term xr is 0, and thus this cross term disappears in the system equations of (22). So the nonlinearity of this system is reduced. In this case, when system (22) is transformed into the Liénard form, the properties (1) − (4) listed in the proof of Lemma 16 allow us to use the Liénard system theory to prove that there is at most one limit cycle in int(I) for µ ∈ (0, µ 1 ). Unfortunately, those properties are not satisfied when a + c = b + d, and thus the existence and the number of limit cycles in this case remain unclear.
Now we are in the position to present the last result of this section about the uniqueness of limit cycle and stability of the interior equilibrium.
PROOF. Under Assumption 5, when a + c = b + d, according to Lemma 16, the system (4) has at most one limit cycle for 0 < µ < µ 1 . As shown in Theorem 6, the eigenvalues of the Jacobian at (x * , r * ) have positive real parts, implying that it is repulsive. Then all trajectories starting in the small neighborhood of this equilibrium (excluding the equilibrium) will diverge, and the existence of infinite many periodic orbits in the neighborhood is excluded. From Lemma 13, one also knows that the boundary ∂I is repelling. According to the Poincaré-Bendixson theorem, for any initial states in int(I) \ {(x * , r * )}, its ω-limit set must be a closed orbit. If there are no limit cycles in int(I) \ {(x * , r * )}, then every trajectory starting from int(I) \ {(x * , r * )} will form a periodic orbit, and an infinite number of periodic orbits will fill int(I) \ {(x * , r * )}. This contradicts the repulsiveness of the interior equilibrium and the boundary [28].
Suppose in int(I) \ {(x * , r * )}, there is some connected region which consists of an infinite number of neutral periodic orbits and is enclosed by two periodic orbits denoted by Γ 1 and Γ 2 . Then Γ 1 and Γ 2 must be limit cycles, which contradicts Lemma 16 regarding the maximal number of limit cycles. Therefore, one can conclude that there is exactly one stable limit cycle in I, and it is attractive for I \ (B t ∪ B b ∪ (x * , r * )).
Note that when a+c = b+d and θ = 1, we have µ 0 = µ 1 . Thus, the second claim is easy to verify according to Theorems 6 and 14. 2 Fig. 3 shows the plot of the unique limit cycle of system (4) when µ changes from 0 to 0.4. The payoff matrices used in the simulations are shown below, and they satisfy the condition a + c = b + d.
The Hopf bifurcation point is computed to be µ 1 = 0.1633. It can be observed that the limit cycle "collides" with the heteroclinic cycle lying on the boundary and the interior equilibrium as µ approaches 0 and 0.1633 respectively. In addition, it disappears after µ passes 0.1633.

Incentive-based Control
In real biological and social systems, the mutation or exploration rate is usually small. For example, an extremely high mutation rate for a biological entity may reach 1/400 per site as reported in [3]. As known from Fig. 3. The plot of a unique limit cycle of dynamics (4) with the given payoff matrices in (25). The x-axis is the mutation rate µ, the blue and red curves are stable and unstable equilibria respectively, and the magenta curves show the shape of stable limit cycles.
Section 3, although the Hopf bifurcation point µ 1 depends on the payoff parameters, limit cycles exist in system (4) for µ < µ 1 . It means that in many situations stable limit cycles persist in the co-evolutionary dynamics of the game and environment. On the one hand, the shapes and positions of limit cycles in nonlinear systems are generally difficult to identify analytically; on the other hand, such limit cycle oscillations, although different than the (oscillating) tragedy of the commons, are not desirable from a control point of view. In addition, considering the public resource management, the issues of how to best govern common-pool resources have long been a common concern in the field of economics [20]. Thus, in this section we intend to study how to design suitable control policies such that the following two objectives are attained: (a). alleviate or even eliminate the oscillation; (b). increase the environmental resource stock.
One needs to choose what the control input should be. Although multiple possible approaches have been proposed in the evolutionary games field, perhaps the most plausible control input is to offer suitable incentives to be added to the payoff of the desired strategies [26]. Now we consider that in each game interaction, the cooperator receives certain constant incentive u ≥ 0 from an external regulating authority.
Recall that the entries on the first row of A(r) are the payoffs to a C player. The external incentive will be added to the first row of (1). Thus, the control input u is incorporated into the previously defined payoff matrix (1) as follows (26) Substituting (26) into equations (4) results in the following system with the control input u To emphasize the effect of the control input u on the system behaviors, we have assumed that the enhancement effect and the degradation effect are balanced by fixing the parameter θ to be 1 (under Assumption 15). In the following part we will study whether such a constant incentive can achieve the desired control objectives.

Stable equilibrium
System (27) has a unique interior equilibrium which is given by (x * c , r * c ) = (1/2, (a+b+2u)/(a+c+b+d)) when 0 ≤ u < (c + d)/2. It is noted that the interior equilibrium is shifted up vertically, i.e., the x-coordinate of this equilibrium is the same as the equilibrium of the system without control, while the value of the r-coordinate increases. Then we evaluate the Jacobian at (x * c , r * c ), which is given by For this matrix, it is easy to identify that the real parts of its eigenvalues are of the same sign, which is the same as the first entry.
When u ≥ (c + d)/2, no equilibria exist in the interior of I, and all the equilibria are located on the sides B t and B b . For system (27), if u ≥ 0, it is easy to identify that the x-coordinate of the equilibrium on B b is in the range x ∈ (1/2, 1) by using a similar approach in Lemma 11. Then the same argument of Lemma 12 is applicable, so one can ensure that this equilibrium is unstable for int(I). Similarly, on the side B t , when 0 < u < (c+d)/2, the equilibria are located in the range x ∈ (0, 1/2), and thus they are unstable. However, when u > (c + d)/2, the situation will be different, and we have the following statement for the equilibria.

PROOF. See Appendix G. 2
It is easy to check the stability of these equilibria by evaluating the Jacobian.

Lemma 19
When u > (c + d)/2, the boundary equilibrium (x * t , 1) is locally asymptotically stable under system (27), while the other equilibria (if they exist) are unstable.
PROOF. One can check the local stability of (x * t , 1) by examining the Jacobian which is given by . Note that the eigenvalues of J * t are actually the two diagonal entries. According to Lemma 18, obviously the eigenvalue For the other eigenvalue Thus, the eigenvalue λ 2 is negative. The fact of two negative eigenvalues ensures that the equilibrium (x * t , 1) is locally asymptotically stable. The instability of other possible equilibria is easy to check since there is always a positive eigenvalue in view of (29). Then, we define a Dulac function which is similar to that in Section 3.1 where the parameters α, β, γ, and δ are the same as in (11) with θ = 1. Denote the vector field of (27) by f c (x, r). After the multiplication of ϕ c (x, r), we compute the divergence of the modified vector field ϕ c f c (x, r) in int(I), which yields with g(x) being the same as in (12), which takes the maximum value atx as shown in (13).
Next we discuss in two cases and show that there always exists some u such that system (27) has a stable equilibrium for the whole interior. Theorem 20 For system (27), when a+c > b+d, there exists some u 1 := ad−bc−(4−g(x))µ(a+b+c+d) , so the coefficient of the term with u in (31), i.e., 1.5 + α, is negative. Thus, if the right hand side of (31) is always negative, which implies that there are no closed orbits in I.
By checking the Jacobian (28), one can see that the interior equilibrium (x * c , r * c ) is stable when As g(x) ≥ 0, we have u 1 ≥ u 2 . Then if u 1 < (c+d)/2, one can choose some u ∈ (u 1 , (c+d)/2) such that the interior equilibrium (x * c , r * c ) exists and is asymptotically stable. Due to the non-existence of closed orbits in I and other stable equilibria, and also because of the repulsiveness of ∂I, one can conclude that almost all trajectories will converge to (x * c , r * c ) asymptotically.
On the other hand, if u 1 ≥ (c + d)/2, one can choose some u > (c + d)/2 which leads to the occurrence of the equilibrium (x * t , 1). In view of (29), the equilibrium (x * t , 1) is locally asymptotically stable. The other equilibria on B t and B b are all unstable, and they are located on the sets whose sufficient close neighborhoods in int(I) are repelling. Hence, it is impossible for trajectories starting from int(I) to converge to these points. As a consequence, (x * t , 1) is the only ω-limit set for all initial states in int(I) from the Poincaré-Bendixson theorem. Similarly, this result extends to the sides B l and B r by following the preceding analysis. 2 Theorem 21 For system (27), when a + c ≤ b + d, the boundary equilibrium (x * t , 1) is almost globally stable under u > (c + d)/2.
PROOF. When a + c ≤ b + d, the first entry of (28) will always be positive for any u if ad−bc 2(a+b+c+d) − 2µ > 0, in which case the equilibrium will always be unstable (saddle). Thus, it is impossible to stabilize the interior equilibrium with any positive u if it is unstable in the uncontrolled system. However, one can choose some u > (c + d)/2 such that there is an equilibrium (x * t , 1) on B t which can be proved to be almost globally stable by using the same argument in the second case of Theorem 20. 2 In Theorems 20 and 21, the relationship between a + c and b + d determines whether the interior equilibrium can be stabilized or not under the positive control input u. In view of (28), one can observe that the stability of the interior equilibrium can be changed from unstable to stable when a + c > b + d under a suitable u, while the stability is not qualitatively affected by u when a + c ≤ b + d.
When there exists an equilibrium that is globally stable for I, obviously the limit cycle oscillation is eliminated. No matter whether this stable equilibrium is in the interior or on the top side. The second designed goal is also achieved, since the r-coordinates of the equilibria are higher. One can find some illustrations in Fig. 4.

Reduced amplitude of the limit cycle
For the case of a + c = b + d, although it is impossible to stabilize the interior equilibrium with a suitable u, one can observe that the amplitude of the oscillation on r-coordinate decreases under some small u.
As a + c = b + d, it is noted that the control input does not affect the stability of equilibrium (x * c , r * c ) according to the Jacobian (28). However, one can transform system (27) into a form of Liénard system as (24) in the same manner. Since the listed conditions are satisfied straightforwardly, according to Lemma (16), for µ < (0, µ 1 ), system (27) has a unique limit cycle when u ∈ [0, (c + d)/2). When there is some intermediate control input u ∈ (0, (c + d)/2), the corresponding interior equilibrium is shifted towards the top side B t . To some extent, the vector field between the interior equilibrium and the top side is "compressed". The position of the limit cycle is shifted up, and the amplitude of the cyclic oscillation on r-coordinate is reduced at the same time. In this sense, the oscillation is alleviated. Some numerical results can be found in Fig. 5.

Conclusions and Future Work
We have investigated the co-evolutionary dynamics of the game and environment when strategies' mutations are taken into account. We used the replicator-mutator model to depict the strategic evolution. By using bifurcation theory we showed that the system at the interior equilibrium undergoes a Hopf bifurcation which generates stable limit cycles; on the other hand, we showed there is a heteroclinic bifurcation which may also generate a stable limit cycle. Our analysis highlights the important role that mutations play in the integrated system dynamics. For the co-evolutionary system, the stable limit cycle corresponds to a sustained oscillation of population's decisions and richness of the environmental resource. Compared with the neutral periodic oscillation or heteroclinic cycle oscillation that have been exhibited in previous studies using the replicator model, such limit cycle oscillations can better explain the phenomena that appear in real biological or social systems. Moreover, we proved that under certain parameter conditions, the system always admits a unique stable limit cycle which is attractive for almost the entire domain. Such a robust limit cycle oscillation may provide implications and theoretical support for relevant studies in biology and sociology. For example, in the microbial context, some cooperating bacteria can produce costly public good, while the defecting bacteria only profit. Thus, the richness of the public good is enhanced by the cooperating bacteria and degraded by the defectors. In the process of bacterial reproduction, mutations are common. The amount of mutant offsprings can have important impact on the dynamic outcomes of the biological systems. For such microbial systems, [2] shows that the dynamics can converge to an interior equilibrium, while [35] shows that persistent heteroclinic oscillations can exist. However, the results obtained in this work theoretically reveal that a moderate limit cycle oscillation can exist in addition to the convergence to a steady state or a heteroclinic cycle. Thus, this work complements the previous studies.
We also considered the incentive-based control applied into the co-evolutionary dynamics, and showed the effectiveness of such control on maintaining the environmental resource. Namely, compared with the uncontrolled system, the incentive will increase the level of environmental resource in most cases. In particular, one can apply some intermediate incentive to let the system have an interior equilibrium that is stable for the whole interior of domain. If the interior equilibrium cannot be stabilized, some large incentive could be put into use to make the system dynamics converge to a boundary equilibrium with full environmental resource.
The obtained results in the control part of this work can provide valuable insights into many practical issues in nature, such as the common-pool resource management, environmental and sustainable development, and ecological diversity conservation. For example, the global climate change is one of the biggest challenges of our times. The strategic decisions of individuals, corporations, and governments have long-term environmental consequence that will, in turn, alter the strategic landscape those parties face [32]. To mitigate the global warming, all concerned parties are called upon to reduce greenhouse gas emissions. In such a socio-ecological system, the reduction of greenhouse gas emissions can be taken as the environmental factor of interest. Policy-makers can design and implement economic incentives aiming at adapting individual decisions to collectively agreed goals [31]. Hence, the incentive-based policy proposed in this work shows the power of incentive-based control policy and leads to some directions in mitigating the global warming issue.
In this paper, we have used a more general and realistic model, replicator-mutator dynamics, to describe the strategic evolution, but the environment dynamics is still described by the simplified logistic model. Thus, it is of great interest to study the game dynamics under different types of feedback of the practical environmental resource, such as renewing or decaying resource. The replicator or replicator-mutator equations describe the evolutionary game dynamics in an infinite well-mixed population. Due to this shortcoming, they are not suitable to study game dynamics in finite or structured populations. Therefore, another direction of future studies is to develop more appropriate individual-based models to study the game dynamics under the environmental feedback. Last but not least, in addition to the incentivebased control, there could exist other kinds of control policies that can be more effective in various specific situations.
A Proof of Lemma 1 PROOF. The equations (10) can be taken as a set of linear equations with respect to unknown variables α, β, γ, and δ. It suffices to show that (10) always has solutions. This can be done by checking the ranks of the coefficient matrix C and augmented matrix [C|D], which are given by One can easily check that rank C ≡ rank [C|D] ≤ 4, which implies that (10) has at least one set of solutions. Actually when −c + d − a + b = 0 and a − b = 0, (10) will have infinite many solutions; otherwise, it will have exactly one solution. PROOF. For the boundary of I, the same argument presented in the proof of Lemma 2 is also applicable to exclude the existence of periodic orbits that lie on or intersect with the boundary. Thus we only need to focus on the case of the interior int(I).
By definition, a limit cycle is an isolated periodic orbit, which is the α-limit set or ω-limit set for some initial points that do not lie on the orbit [ [8]]. Suppose there is a limit cycle Γ 0 in int(I), and an arbitrary trajectory starts from the initial point (x 0 , r 0 ) and spirals toward Γ 0 as t → +∞ or t → −∞. Then according to [8,Corollary 10.1], (x 0 , r 0 ) must have a neighborhood Ω such that Γ 0 is the ω-limit set for all points in it. Due to the fact that V (x, r) is the first integral for system (4), V (x, r) must be constant in Ω because of the continuity. However, the partial derivatives of V (x, r) only equal 0 at (x * , r * ), which leads to a contradiction. Therefore, system (4) cannot have limit cycles in I when µ = 0. 2

C Computation of first Lyapunov coefficient
The following is the process of calculating 1 as presented in [13]. Consider a planar autonomous systeṁ y = f (y, µ) where y ∈ R 2 and µ ∈ R. Let A 0 be the Jacobian matrix evaluated at the equilibrium y 0 and the Hopf bifurcation point µ h . Suppose A 0 has two purely imaginary complex conjugate eigenvalues, given by ±iω 0 . The algorithm for computing 1 can be summarized as follows: (1) compute the complex vectors q, p ∈ C 2 such that where p, q =p T q is the inner product; (2) compose the expression with a complex variable z H(z,z) = p, f (y 0 + zq +zq, µ h ) , and find the coefficients g 20 , g 11 , and g 21 in the formal Taylor series (3) the first Lyapunov coefficient 1 (y 0 , µ h ) is given by (ig 20 g 11 + ω 0 g 21 ).
Consider system (4). When µ = µ 1 , we have Then following the above process, we obtain PROOF. One only needs to check the case of B b , since the other case is analogous. When r(0) = 0 initially, we haveṙ = 0, and thus one can focus on the following onedimensional dynamics about x: Thus, the equilibrium's x-coordinates are the solutions of the following equation The left hand side (LHS) of (D.2) is non-negative in the interval x ∈ [0, 1] and is exactly equal to 0 at x = 0, 1, while the right hand side (RHS) of (D.2), µ(2x − 1), is negative for x ∈ [0, 1/2) and positive for x ∈ (1/2, 1] when µ ∈ (0, 1]. Thus, the LHS and RHS of (D.2) will be equal at some x ∈ (1/2, 1). Now, we will show that the LHS and RHS of (D.2) will equal for exactly one value of x in the range at (1/2, 1). Denote the LHS of (D.2) by φ(x). Its first derivative is φ (x) = −3(a − b)x 2 + 2(a − 2b)x + b, and the second derivative is φ (x) = 2(a − 2b) − 6(a − b)x. We now discuss in the following cases: 1). When a − b = 0, one has φ (x) = −2b < 0 as b > 0, which implies φ(x) is concave. In this case, φ(x) equals the RHS of (D.2) exactly once when When θ ≥ 1, the latter statement of the lemma holds trivially since 1/(θ + 1) ≤ 1/2. For the case θ < 1, one can substitute x = 1/(θ + 1) into the RHS of (D.1), which yields µ(θ−1)(θ+1) 2 +θ(a+θb) . This value is positive for any µ ∈ (0, 1] because of (6). Since the RHS of (D.1) is negative at x = 1, then due to continuity, the RHS of (D.1) equals 0 when x ∈ (1/(θ + 1), 1). As stated before, the vectors on the left side B l and the right side B r point inwards to int(I). This implies that these two sides are repelling. Therefore, the subsequent discussion is restricted to the top and bottom sides B t and B b , which are positively invariant.
We define a partial distance function F : [0, 1] → R as below: which indicates the distance of the system states (x, r) to the top and bottom sides B t and B b . It is positive in int(I) and only equals zero on these two sides. In addition, as r increases from 0 to 1, the function value F (r) monotonically increases from 0 to the maximal value at r = 1/2 and then monotonically decreases to 0 at r = 1.
The derivative of F with respect to time t iṡ F (r) =ṙ(1 − 2r) = F (r)Φ(x, r), where Φ(x, r) := (1 − 2r)[(θ + 1)x − 1]. We want to show that there exists aδ > 0 such that lim inf t→∞ F (r(t)) > δ, which implies that B t and B b are repelling. We only discuss the case of B b here, as the analysis for B t is similar and thus omitted. Denote the unique interior equilibrium by q * = (x * , r * ), as given before. Choose 0 < r < min{r * , 1/2}, and let m = F (r ) = r (1 − r ). We define a lower region by I l (m) := {(x, r) ∈ I : 0 < F (r) ≤ m, 0 < r ≤ r }. Note that I l (m) ∩ B b = ∅, and the condition 0 ≤ F (r) ≤ m is redundant in this definition, but it simply indicates that the set is related to the partial distance function evaluated at r = r . Denote the solution of the dynamics (4) by ξ(t) := (x(t), r(t)), given some initial condition ξ(0) = (x(0), r(0)) ∈ I.
First we show that if ξ(0) ∈ I l (m), then there exists a time instant T > 0, such that ξ(T ) / ∈ I l (m). Define a line L c in the lower region I l (m) by L c := {(x, r) ∈ I l (m) : x = 1/(θ + 1)}. Note that on this line,ṙ = 0 andẋ ≥ > 0 for some positive constant , and thus the vectors point horizontally towards the positive xdirection. The line L c divides the lower region I l (m) into the left region J L := {(x, r) ∈ I l (m) : 0 ≤ x < 1/(θ+1)} and the right region J R := {(x, r) ∈ I l (m) : 1/(θ + 1) < x ≤ 1}. Note thatḞ = F (r)Φ(x, r) > 0 in J R whilė F < 0 in J L . Therefore, if ξ(0) ∈ J R , the trajectory of r(t) will increase, and there exists a T > 0 such that ξ(T ) / ∈ I l (m) (due to the directions of vectors on the line L c , the right side B r andḞ (r) ≥ > 0 for some positive constant in a sufficiently large compact set Ω satisfying {(x, r) ∈ J R : r = r } ⊂ Ω ⊂ J R and ξ(0) ∈ Ω). If ξ(0) ∈ J L , then the trajectory ξ(t) will not always stay in J L due to the Poincaré-Bendixson theorem, since there are no equilibria in J L . AsḞ (r) < 0 in J L , the trajectory of r(t) will decrease, and x(t) will cross the center line L c , and enter the right region J R . Due to the uniqueness of solutions of (4), similar to the previous analysis of the right region, there exits a T > 0 such that ξ(T ) / ∈ I l (m).
Let n = me ΦminT . Therefore we have shown that the trajectory will not reach I l (n) for t ∈ (t 0 , t 0 + T ). Since T is the supremum of the duration when the solution could stay in the lower region I l (m), then there exists T ∈ (t 0 , t 0 +T ] such that ξ(T ) / ∈ I l (m). Therefore, the trajectory ξ(t) will leave the lower region I l (m) without reaching I l (n) ⊂ I l (m). Repeating this argument, we conclude that the trajectory ξ(t) will never reach I l (n). Thus we have shown that there exists 0 <δ < n such that lim inf t→∞ F (r(t)) >δ, and thus the lower side B b is repelling. (1/4 −x 2 ) 2 (r + r * )(1 −r − r * ) ; (4) F (x,r) = 2µx − 8µ 1x (1/4 −x 2 )) (1/4 −x 2 )(r + r * )(1 −r − r * ) .

G Proof of Lemma 18
PROOF. In view of the system equations (27), the xcoordinates of the equilibria on B t are obtained by solving the following polynomial equation The RHS of (G.1), µ(2x − 1), is a linear function of x which equals 0 at x = 1/2. It is negative for x ∈ (0, 1/2) and positive for x ∈ (1/2, 1) when µ ∈ (0, 1]. The LHS of (G.1) equals 0 at the two end points x = 0, 1.
When u > (c + d)/2 and c = d, the LHS will be always positive for x ∈ (0, 1), such that the root of (G.1) is in the range (1/2, 1). The value of LHS will change the sign from positive to negative (resp. negative to positive) as x increases from 0 to 1 when c > d (resp. c < d), and it equals 0 also at x = (d − u)/(d − c). One has (d − u)/(d − c) < 1/2 and (d − u)/(d − c) > 1/2 respectively for c < d and c > d. For both cases, there will be some x ∈ (1/2, 1) such that the RHS and LHS are equal as shown geometrically in Fig. (G.1). Therefore, one can conclude that there is always one equilibrium with xcoordinate being 1/2 < x * t < 1.
Next, only in the case c < d, the curves of the LHS and RHS can have one or two intersections in the range 0 < x < 1/2. Thus, the other possible equilibria on B t can only be located in x ∈ (0, 1/2).