Limit Cycles Analysis and Control of Evolutionary Game Dynamics with Environmental Feedback

We consider games played in a changing environment, which is characterized by the richness $r\in[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 $\begin{bmatrix}R_{2}&S_{2}\\ T_{2}&P_{2}\end{bmatrix}$ and $\begin{bmatrix}R_{1}&S_{1}\\ T_{1}&P_{1}\end{bmatrix}$ 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.

The two-player game governed by the payoff matrix $A(r)$ in (1) are played in a well-mixed and infinite population. We denote the proportion of individuals choosing C by the variable $x\in[0,1]$. Since there are only two strategies, the population state can be represented by the vector $\mathbf{x}=[x~{}~{}1-x]^{\top}$. Assume the individuals with different strategies can mutate into each other with an identical probability $\mu\in[0,1]$: a fraction of cooperators, $\mu x$, spontaneously change to choose strategy D; a fraction of defectors, $\mu(1-x)$, change to choose strategy C at the same time. We focus on the dynamics of $x$, which can be described by the following replicator-mutator equation

where $(A(r)\mathbf{x})_{1}$ is the first entry of $A(r)\mathbf{x}$ and represents the fitness of choosing strategy C, and $\mathbf{x}^{\top}A(r)\mathbf{x}$ is the average fitness at the population state $\mathbf{x}$. The term $-\mu x$ represents the mutation from strategy C towards D, while $\mu(1-x)$ captures the mutation in the opposite direction, from D towards C.

To model the evolution of the environmental change, we use the standard logistic model

where the parameter $\theta>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 $[\theta x-(1-x)]$, and the dynamics are restricted in the unit interval by $r(1-r)$. When $\theta=1$, the enhancement effect and the degradation effect are balanced. And if $\theta>1$ or $0<\theta<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 game-environment 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 $\mu=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.

First, we define the $\mathbf{C}^{1}$ function $\varphi(x,r):\text{int}(\mathcal{I})\rightarrow\mathbb{R}$ by

with the parameters $\alpha$, $\beta$, $\gamma$, and $\delta$ to be determined later. Note that $\varphi(x,r)$ is strictly positive in $\text{int}(\mathcal{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)]^{\top}$, and then the divergence of the modified vector field $\varphi f(x,r)$ on $\text{int}(\mathcal{I})$ is given by

Let all the coefficients of the powers of $x$, $r$ in (9) (except for the term containing $\mu$) be zero, and we obtain the following set of equations

We first show that the equation set (10) always has solutions.

See Appendix A. ∎

In fact, if one of $-c+d-a+b$ and $a-b$ is not zero, (10) will have exactly one solution

with $\zeta=a+c+\theta b+\theta d>0$. Otherwise, (10) will have infinitely many solutions including (11). So we fix the parameters $\alpha$, $\beta$, $\gamma$, and $\delta$ 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).

The proof is direct by applying Bendixson-Dulac criterion [36, Theorem 4.1.2]. Since the boundary $\partial\mathcal{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 $\text{int}(\mathcal{I})$ which is simply connected.

In view of the fact that $\alpha$, $\beta$, $\gamma$, and $\delta$ are given as in (11), one can calculate the divergence of $\varphi f(x,r)$

From (11) one knows $\alpha\in(-2,-1)$, and thus the term $g(x):=\frac{\mu(1-2x)(3x+\alpha)}{x(1-x)}$ takes its maximum value at $\bar{x}\in(0,1)$ given by

Rearranging (12) yields

It is easy to check that $\frac{(1-2\alpha)\bar{x}+\alpha-4\bar{x}^{2}}{\bar{x}(1-\bar{x})}$ is negative for any value of $\alpha\in(-2,-1)$. Since $\varphi(x,r)$ is positive in $\text{int}(\mathcal{I})$ by definition, if $ad-bc>0$, from (14) one has $\text{div}~{}\varphi f(x,r)\leq 0$ when

And because the divergence is not identically zero in $\text{int}(\mathcal{I})$, from the Bendixson-Dulac criterion, one knows that no closed orbits can exist in $\text{int}(\mathcal{I})$ for $\mu\geq\mu_{0}$.

For the second case, it is easy to prove since the right hand side of (14) will always be negative for any $\mu\geq 0$. For the third case, the right hand side of (14) is always negative when $\mu>0$. ∎ In Lemma 2, it is noted that the conditions for the non-existence of periodic orbits depend on the relationship between $ad$ and $bc$. The reason is as follows: when $ad-bc<0$, the divergence $\text{div}~{}\varphi f(x,r)$ will remain positive for any $\mu\in[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 $\mu$ such that $\text{div}~{}\varphi 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 $\mu=0$ when $ad-bc=0$. This possibility has been shown in [35], [4]. Despite the possibility of closed orbits under certain conditions, for $\mu=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.

The full proof is presented in [5]. See Appendix B for a sketch of the proof.

In this section we take $\mu$ 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 $\mu$ 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.

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

where

and

Note that the off-diagonal entries

and

because of (6). One can calculate the eigenvalues of $J^{*}$

It is easy to check that when $\mu=\mu_{1}:=\frac{\theta^{2}(ad-bc)}{\Delta}$ the eigenvalues are purely imaginary, i.e.,

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 $\mu=\mu_{1}$. Since $\mu$ varies in the interval $(0,1]$, to guarantee the occurrence of Hopf bifurcation, the critical parameter value $\mu_{1}$ should also be constrained in this interval, i.e., $0<\mu_{1}\leq 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 $\mu$ in the vicinity of $\mu_{1}$. Whether the bifurcated periodic orbits are stable or not is determined by the so-called first Lyapunov coefficient $\ell_{1}$ at $(x^{*},r^{*})$ when $\mu=\mu_{1}$. If $\ell_{1}<0$, then these periodic orbits are stable limit cycles; if $\ell_{1}>0$, the periodic orbits are repelling. Following the computation procedure in Appendix C, we obtain $\ell_{1}(\mu_{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 $\mu<\mu_{1}$, $(x^{*},r^{*})$ is unstable since the eigenvalues have positive real parts. On the other hand, the interior equilibrium is asymptotically stable for $\mu>\mu_{1}$ with the negative eigenvalues. Hence, the stable limit cycles exist for $\mu<\mu_{1}$ in the vicinity of $\mu_{1}$. Since the bifurcation is associated with stable limit cycles, it is a supercritical Hopf bifurcation. ∎ 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\leq 0$, $J_{11}$ will remain non-positive no matter what value $\mu$ takes. Thus, as $\mu$ 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 $\mu\in[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(\sqrt{\mu_{1}-\mu})$. Since Hopf bifurcation is a local bifurcation, the obtained results only hold for $\mu$ in the vicinity of $\mu_{1}$. Next, we will show that system (4) exhibits another type of bifurcation.

In the case $\mu=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 $\mu>0$ in the vicinity of $0$.

When $\mu=0$, system (4) has four equilibria on the boundary $\partial\mathcal{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 $\Lambda$, which connects the four equilibria and has the orientation $E_{1}\rightarrow E_{2}\rightarrow E_{3}\rightarrow E_{4}\rightarrow E_{1}$, 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 $\mu>0$, the four corners are no longer equilibria, and the heteroclinic cycle $\Lambda$ is “broken”. Now we are going to show that the system (4) undergoes a heteroclinic bifurcation at $\mu=0$ when certain condition holds, and a limit cycle with the same stability of the corresponding heteroclinic cycle may be generated from $\Lambda$ for $\mu>0$ sufficiently close to $0$.

The proof is straightforward by applying the theorem of the heteroclinic bifurcation [16, Theorem 2.3.2, 2.3.3]. When $a>\theta$, $b<1$, $c<\theta$, and $d>1$ hold, which implies $ad>bc$. Then according to Lemma 7, there is a heteroclinic cycle when $\mu=0$ and it is stable.

In addition, all the traces $T_{i}$ of the Jacobians of the equilibria are negative, while their product $\prod_{i=1}^{4}T_{i}$ is positive. Then according to [16, Theorem 2.3.3], there exist sufficiently small $\sigma,\epsilon>0$ such that for $0<\mu<\epsilon$, system (4) has at most one limit cycle in the internal $\sigma$-neighbourhood of $\Lambda$, and this limit cycle (if it exists) is stable due to the stability of $\Lambda$. ∎ 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 $\mathcal{I}$ for any $\mu\geq 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 $\Lambda$. 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 $\theta$ is chosen to be $1$, and thus the interior equilibrium is fixed and independent of $\mu$. In view of Theorem 6, the Hopf bifurcation point is $\mu_{1}=0.1543$.

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 $\mu\in(0,\mu_{1})$, it is necessary to analyze the dynamics of (4) for all such $\mu$.