Bimatrix Replicator Dynamics with Periodic Impulses

Abstract

This paper investigates the bimatrix replicator dynamics with periodic impulses. In biological system, impulsive perturbations may due to the occurrence of an unfavorable physical environment, or due to the seasonal life history effects in the physiological and reproductive mechanisms of the population. We show that impulsive perturbations can lead to complicated dynamical behaviors. On the one hand, the system can have multiple $\tau$-periodic solutions, where the lower bound of the number of solutions is increasing linearly in the impulsive period $\tau$. On the other hand, for shorter impulsive period, we provide a differential approximation for the impulsive dynamical system. By analyzing the resulting differential dynamics, we show that the interior equilibrium (which corresponds to a $\tau$-periodic solution) must be globally stable if it exists. Furthermore, when the impulsive effect is weak, all interior trajectories of the impulsive dynamics evolve to a small neighborhood of the interior equilibrium of the bimatrix replicator dynamics. In summary, longer impulsive period causes an increase in the complexity of the evolutionary process and shorter period promotes evolutionary stability of the interior equilibrium where multiple strategies can coexist.

Appendix

1.1 Note A1. Population Dynamics with Periodic Impulses

To derive the time evolution of population sizes with periodic impulses, we assume that all A-individuals have the same density-dependent background fitness $W_A(n)$, where n is the size of population A, and all B-individuals have the same background fitness $W_B(m)$, where m is the size of population B [18, 31]. Specifically, we consider $W_A(n)=R_A(1-\frac{n}{K_A})$ and $W_B(m)=R_B(1-\frac{m}{K_B})$, i.e., the background fitness has the form of logistic growth [18]. Let $n_1$ and $n_2$ be the numbers of $A_1$ and $A_2$ individuals in population A, respectively, and $m_1$ and $m_2$ be the numbers of $B_1$ and $B_2$ individuals in population B, respectively. Then the time evolution of $n_1, n_2, m_1, m_2$ with periodic impulses can be written as

$$\begin{aligned} \begin{aligned}&\frac{\hbox {d}n_1}{\hbox {d}t}=n_1\left( a_{11}\frac{m_1}{m}+a_{12}\frac{m_2}{m}+R_A\left( 1-\frac{n}{K_A}\right) \right) , \quad \text {for } t\ne {k\tau },\\&\frac{\hbox {d}n_2}{\hbox {d}t}=n_2\left( a_{21}\frac{m_1}{m}+a_{22}\frac{m_2}{m}+R_A\left( 1-\frac{n}{K_A}\right) \right) , \quad \text {for } t\ne {k\tau },\\&\frac{\hbox {d}m_1}{\hbox {d}t}=m_1\left( b_{11}\frac{n_1}{n}+b_{21}\frac{n_2}{n}+R_B\left( 1-\frac{m}{K_B}\right) \right) , \quad \text {for } t\ne {k\tau },\\&\frac{\hbox {d}m_2}{\hbox {d}t}=m_2\left( b_{12}\frac{n_1}{n}+b_{22}\frac{n_2}{n}+R_B\left( 1-\frac{m}{K_B}\right) \right) , \quad \text {for } t\ne {k\tau },\\&\Delta n_1=-\gamma _1 n_1, \quad \text {for } t={k\tau },\\&\Delta n_2=-\gamma _2 n_2, \quad \text {for } t={k\tau },\\&\Delta m_1=-\psi _1 m_1, \quad \text {for } t={k\tau },\\&\Delta m_2=-\psi _1 m_2, \quad \text {for } t={k\tau }. \end{aligned} \end{aligned}$$

(22)

We note that the frequency dynamics Eq. (3) is well defined only for $m>0$ and $n>0$. Thus, it is essential to provide sufficient conditions that $m=0$ or $n=0$ are unstable under Eq. (22). For $n \rightarrow 0$, the change of $n_i$ in an impulsive interval can be approximated as $\frac{\hbox {d}n_i}{n_i}=(a_{i1}\frac{m_1}{m}+a_{i2}\frac{m_2}{m}+R_A)\hbox {d}t$ with integral

$$\begin{aligned} \ln \frac{n_i(\tau )}{n_i(0)}=\int _{0}^{\tau }(a_{i1}\frac{m_1}{m}+a_{i2}\frac{m_2}{m}+R_A)\hbox {d}t \ge \tau (\check{a_i}+R_A), \end{aligned}$$

(23)

where $\check{a_i}=\min \{a_{i1},a_{i2}\}$. The boundary $n=0$ is unstable if $n_1(\tau )+n_2(\tau )-n_1(0)-n_2(0)>\gamma _1 n_1(\tau )+\gamma _2 n_2(\tau )$, i.e., the increase in the population size in an impulsive interval is greater than the decrease at the moment of jump [2]. This inequality holds for $\gamma _i<1-\exp (-\tau (\check{a_i}+R_A))$ with $i=1,2$. Similarly, the boundary $m=0$ is unstable if $\psi _j<1-\exp (-\tau (\check{b_j}+R_B))$ with $j=1,2$, where $\check{b_j}=\min \{b_{1j},b_{2j}\}$.

Lemma A1

Let us denote the cyclic orbit of Eq. (2) that passes $(p^*, q^*)$ by $H(p, q)=M^*$, and denote the period of this orbit by $T(p^*, q^*)$. (1) $T(p^*, q^*)$ takes its minimum $T_{min}=2\pi \sqrt{-\frac{(a_1+a_2)(b_1+b_2)}{a_1a_2b_1b_2}}$ when $(p^*, q^*)=(p_0, q_0)$. (2) Furthermore, for $q^*=q_0$ and $p^*>p_0$, T is monotonically increasing in $p^*$ and $\lim \limits _{p^* \rightarrow 1}T(p^*, q^*)=\infty$.

Proof

Part (b) of Lemma 1 is directly from [39] Theorem 1 and example 2. Thus, $T(p^*, q^*)$ takes its minimum when $(p^*, q^*)=(p_0, q_0)$. We next calculate $T(p^*, q^*)$ at this equilibrium.

When the dynamical system Eq. (2) is near the interior equilibrium, it could be approximated as

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}p}{\hbox {d}t}=a(q_0-q),\\ \frac{\hbox {d}q}{\hbox {d}t}=b(p_0-p), \end{aligned} \end{aligned}$$

(24)

where $a=p_o(1-p_0)(a_1+a_2)$ and $b=q_o(1-q_0)(b_1+b_2)$. By setting $y=q_0-q$ and $x=p_0-p$, we obtain

$$\begin{aligned} \left( \begin{array}{c}\frac{\hbox {d}x}{\hbox {d}t}\\ \frac{\hbox {d}y}{\hbox {d}t}\end{array}\right) =\left( \begin{array}{cc}0 &{} -a \\ -b &{} 0\end{array}\right) \left( \begin{array}{c}x\\ y\end{array}\right) , \end{aligned}$$

(25)

where the eigenvalues of the matrix $\left( \begin{array}{cc}0 &{} -a \\ -b &{} 0\end{array}\right)$ are $\lambda _1,\lambda _2=\pm {i\sqrt{-ab}}$. Clearly, the solution of Eq. (25) can be written as

$$\begin{aligned} \left( \begin{array}{c}x(t) \\ y(t)\end{array}\right) =x_1\exp (i\sqrt{-ab}t)+x_2\exp (-i\sqrt{-ab}t), \end{aligned}$$

(26)

where $x_1$ and $x_2$ are eigenvectors corresponding to eigenvalues $\lambda _1$ and $\lambda _2$, and the period of this solution is $2\pi /\sqrt{-ab}$. Thus, the period of cyclic orbits near the interior equilibrium $(p_0,q_0)$ is

$$\begin{aligned} T_{\hbox {min}}=2\pi \sqrt{-\frac{(a_1+a_2)(b_1+b_2)}{a_1a_2b_1b_2}}. \end{aligned}$$

(27)

$\square$

Lemma A2

Define ${{\widetilde{H}}}(p,q)=H(p+\Delta p,q+\Delta q)$. ${{\widetilde{H}}}(p,q)=M$ (with $M<0$) is a closed curve. Furthermore, there exists a unique $M_0$ such that $H(p,q)=M$ and ${{\widetilde{H}}}(p,q)=M$ have at least two intersection points for $M<M_0$.

Proof

For convenience, we denote $p+\Delta p$ and $q+\Delta q$ by $p^+$ and $q^+$, respectively. For the relationship between p, q and $p^+, q^+$,

$$\begin{aligned} \begin{aligned} p^+=\frac{(1-\gamma _1)p}{1-\gamma _1p-\gamma _2(1-p)},\quad q^+=\frac{(1-\psi _1)q}{1-\psi _1q-\psi _2(1-q)}, \end{aligned} \end{aligned}$$

(28)

and

$$\begin{aligned} \begin{aligned} p=\frac{(1-\gamma _2)p^+}{1-\gamma _1+(\gamma _1-\gamma _2)p^+},\quad q=\frac{(1-\psi _2)q^+}{1-\psi _1+(\psi _1-\psi _2)q^+}. \end{aligned} \end{aligned}$$

(29)

Then, we have

$$\begin{aligned} \frac{\hbox {d}p^+}{\hbox {d}t}= & {} \frac{\hbox {d}p^+}{\hbox {d}p}\frac{\hbox {d}p}{\hbox {d}t} =\frac{(1-\gamma _2)(1-\gamma _1)p^+(1-p^+)}{1-\gamma _1+(\gamma _1-\gamma _2)p^+}\bigg [a_2-\frac{(1-\psi _2)q^+(a_1+a_2)}{1-\psi _1+(\psi _1-\psi _2)q^+}\bigg ],\nonumber \\ \frac{\hbox {d}q^{+}}{\hbox {d}t}= & {} \frac{\hbox {d}q^{+}}{\hbox {d}q}\frac{\hbox {d}q}{\hbox {d}t} =\frac{(1-\psi _2)(1-\psi _1)q^+(1-q^+)}{1-\psi _1+(\psi _1-\psi _2)q^+}\bigg [b_2-\frac{(1-\gamma _2)p^+(b_1+b_2)}{1-\gamma _1+(\gamma _1-\gamma _2)p^+}\bigg ].\nonumber \\ \end{aligned}$$

(30)

Specifically, Eq. (30) is equivalent to a Hamilton system

$$\begin{aligned} \begin{aligned}&\frac{\hbox {d}p^+}{\hbox {d}t}=\frac{a_2-\frac{q^+(1-\psi _2)(a_1+a_2)}{1-\psi _1+(\psi _1-\psi _2)q^+}}{(1-\psi _2)(1-\psi _1)\frac{q^+(1-q^+)}{1-\psi _1+(\psi _1-\psi _2)q^+}},\\&\frac{\hbox {d}q^+}{\hbox {d}t}=\frac{b_2-\frac{p^+(1-\gamma _2)(b_1+b_2)}{1-\gamma _1+(\gamma _1-\gamma _2)p^+}}{(1-\gamma _2)(1-\gamma _1)\frac{p^+(1-p^+)}{1-\gamma _1+(\gamma _1-\gamma _2)p^+}}. \end{aligned} \end{aligned}$$

(31)

where this system has an unique interior equilibrium $(p_0^+, q_0^+)=((\frac{b_{2}}{b_{1}+b_{2}})^+,(\frac{a_{2}}{a_{1}+a_{2}})^+)$ and a constant of motion $H(p^+,q^+)={{\widetilde{H}}}(p,q)$. Thus, ${{\widetilde{H}}}(p,q)=M$ (with $M<0$) is a closed curve that surrounds $(p_0^+, q_0^+)$.

We note that both $H(p,q)=M$ and ${{\widetilde{H}}}(p,q)=M$ are closed curves. Furthermore, they have no intersection point for $M_{\max }=H(p_0,q_0)$ because in this case $H(p,q)=M$ is close to $(p_0,q_0)$ and ${{\widetilde{H}}}(p,q)=M$ is close to $(p_0^+,q_0^+)$. Thus, there exists $M_0<M_{\max }$ such that $H(p,q)=M_0$ and ${{\widetilde{H}}}(p,q)=M_0$ are tangent. In addition, ${{\widetilde{H}}}(p,q)=M$ is never included in ${{\widetilde{H}}}(p,q)=M$ (and vice versa) since $p^+<p^-$ and $q^+<q^-$. Therefore, for $M<M_0$, $H(p,q)=M$ and ${{\widetilde{H}}}(p,q)=M$ have at least two intersection points. $\square$

Lemma A3

$\tau _{10}$ and $\tau _{20}$ can be seen as functions of T, and they change continuously in T.

Proof

Lemma A2 shows that $H(p,q)=M$ and ${{\widetilde{H}}}(p,q)=M$ have (at least) two intersection points for $M<M_0$. We denote the two intersection points by $(p_1^-, q_1^-)$ and $(p_2^-, q_2^-)$, and denote their integral times by $\tau _{10}$ and $\tau _{20}$ (with $0<\tau _{10}<\tau _{20}<T$), respectively.

Without loss of generality, we focus on the properties of $\tau _{10}$, and the properties of $\tau _{20}$ can be analyzed analogously. Since T is a continuously decreasing functions of M [39], we only need to prove that $\tau _{10}$ is a continuous function of M.

From the expressions of $H(p,q)=M$ and ${{\widetilde{H}}}(p,q)=M$, both of them are continuous functions of p and q. Specifically, $\frac{\partial H(p,q)}{\partial p}$, $\frac{\partial H(p,q)}{\partial q}$, $\frac{\partial {{\widetilde{H}}}(p,q)}{\partial p}$, and $\frac{\partial {{\widetilde{H}}}(p,q)}{\partial q}$ do not equal to 0 almost always. Thus, their intersection point $(p_1^-, q_1^-)$ must change continuously on M. Otherwise, a sudden change on $(p_1^-, q_1^-)$ will cause a jump on M. Besides, from Eq. (28), $(p_1^+, q_1^+)$ also changes continuously on M.

We note that $\tau _{10}$ can be expressed as the following integral on curve $H(p,q)=M$

$$\begin{aligned} \tau _{10}=\int _{(p_1^+, q_1^+)}^{((p_1^-, q_1^-)}\frac{\hbox {d}s}{\sqrt{P^2(p,q)+Q^2(p,q)}}, \end{aligned}$$

(32)

where s is the element of arc. Since both $(p_1^-, q_1^-)$ and $(p_1^+, q_1^+)$ change continuously on M, $\tau _{10}$ is a continuous function of M. This completes the proof. $\square$

Lemma A4

Equation (17) has no periodic solutions in region $(0,1)\times (0,1)$ for arbitrary values of $a_1, a_2, b_1, b_2$.

Proof

Let us define Dulac function $B=\frac{1}{p(1-p)q(1-q)}$. From the Bendixson–Dulac theorem, no periodic solutions can exist in region $(0,1)\times (0,1)$ since

$$\begin{aligned} \frac{\partial (BP)}{\partial p}+\frac{\partial (BQ)}{\partial q}=-\frac{U(p)^2}{q(1-q)\tau }-\frac{V(q)^2}{p(1-p)\tau }<0 \end{aligned}$$

(33)

for all possible $0<p<1$,$0<q<1$. $\square$

Lemma A5

For BS game with $a_2-\frac{\gamma _1-\gamma _2}{\tau (1-\gamma _1)}>0$ and $b_1+\frac{\psi _1-\psi _2}{\tau (1-\psi _1)}<0$, all trajectories for Eq. (17) with initial values in $(0,1)\times (0,1)$ will converge to $({\widetilde{p}},{\widetilde{q}})$.

Proof

For any initial value $(p^*, q^*) \in (0,1)\times (0,1)$, there exists $\varepsilon >0$ such that $(p^*, q^*)$ and $({\widetilde{p}},{\widetilde{q}})$ are in $D=[\varepsilon ,1-\varepsilon ] \times [\varepsilon ,1-\varepsilon ]$. Define

$$\begin{aligned} \begin{aligned} V(p, q)&=(b_1+b_2)\int \frac{{\widetilde{p}}-p}{p(1-p)}\hbox {d}p+(a_1+a_2) \int \frac{{\widetilde{q}}-q}{q(1-q)}\hbox {d}q \\&=(b_1+b_2)[{\widetilde{p}}\ln \frac{p}{1-p}+\ln (p-1)]-(a_1+a_2) [{\widetilde{q}}\ln \frac{q}{1-q}+\ln (q-1)]. \end{aligned} \end{aligned}$$

(34)

Clearly, V(p, q) is a continuous and bounded function in D. We next show that V(p, q) is a Lyapunov function for Eq. (17).

$$\begin{aligned} \frac{\hbox {d}V(p, q)}{\hbox {d}t}= & {} (b_1+b_2)\frac{{\widetilde{p}}-p}{p(1-p)}\frac{\hbox {d}p}{\hbox {d}t}-(a_1+a_2) \frac{{\widetilde{q}}-q}{q(1-q)}\frac{\hbox {d}q}{\hbox {d}t} \nonumber \\= & {} (b_1+b_2)({\widetilde{p}}-p)[a_2-q(a_1+a_2)+\frac{U(p)}{\tau }]\nonumber \\&-(a_1+a_2) ({\widetilde{q}}-q)[b_2-p(b_1+b_2)+\frac{V(q)}{\tau }]. \end{aligned}$$

(35)

Notice that ${\widetilde{p}}$ and ${\widetilde{q}}$ satisfy

$$\begin{aligned} a_2-{\widetilde{q}}(a_1+a_2)+\frac{U({\widetilde{p}})}{\tau }=b_2-{\widetilde{p}} (b_1+b_2)+\frac{V({\widetilde{q}})}{\tau }=0, \end{aligned}$$

(36)

Equation (35) can be rewritten as

$$\begin{aligned} \frac{\hbox {d}V(p, q)}{\hbox {d}t}=(b_1+b_2)({\widetilde{p}}-p)\frac{U(p)-U({\widetilde{p}})}{\tau }-(a_1+a_2)({\widetilde{q}}-q)\frac{V(q)-V({\widetilde{q}})}{\tau }.\nonumber \\ \end{aligned}$$

(37)

Since $\frac{\hbox {d}U(p)}{\hbox {d}p}=-(U(p))^2<0$ and $\frac{\hbox {d}V(q)}{\hbox {d}q}=-(V(q))^2<0$, we must have $({\widetilde{p}}-p)(U(p)-U({\widetilde{p}}))\ge 0$ and $({\widetilde{q}}-q)(V(q)-V({\widetilde{q}}))\ge 0$, where the equalities hold only for $(p,q)=({\widetilde{p}},{\widetilde{q}})$. This implies $\frac{\hbox {d}V(p, q)}{\hbox {d}t}\le 0$ in D, where the equality holds only for $(p,q)=({\widetilde{p}},{\widetilde{q}})$. Applying the Lasalle’s invariance principle [22], V(p, q) is a Lyapunov function for Eq. (17) and the trajectory with initial value $(p^*, q^*) \in D$ will converge to $({\widetilde{p}},{\widetilde{q}})$ as $t\rightarrow \infty$. $\square$