Discrete-time predator-prey model with flip bifurcation and chaos control

1.   Introduction

It is well-know that in mathematical biology, discrete-time models described by difference equations are more reasonable as compared to corresponding continuous models. The reasons are that in the case of non-overlapping generation discrete models are more realistic than continuous ones, and also these models provide more efficient computational models for numerical simulations as compared to continuous-time models ^[1,2,3]. Among these mathematical models, predator-prey systems have received reasonable attraction during the last few decades. For instance, Yan et al. ^[4] investigated the stability of fixed points, flip and Neimark-Sacker bifurcations, and chaotic behavior of the following $2$-dimensional discrete-time predator-prey model with allee effect in the prey:

Zhao and Yan ^[5] investigated the existence and local stability of fixed points, flip and Neimark-Sacker bifurcations of following discrete predator-prey model with modified Holling-Tanner functional response:

Fang and Li ^[6] investigated the existence and local stability of fixed points, bifurcations and complex dynamical behaviors of the following discrete-time predator-prey model with a strong Allee effect on the prey and a ratio-dependent functional response:

Kangalgi and Kartal ^[7] investigated the existence of equilibrium points, local stability and Neimark-Sacker bifurcation of the following host-parasitoid model with Hassell growth function:

Li and Shen ^[8] investigated the dynamics and bifurcations of the following continuous-time predator-prey model with double Allee effects and time delays:

Stápán ^[9] investigated the stability properties of the zero equilibrium solution, existence of stable and unstable oscillations of the following predator-prey model:

Cheng and Cao ^[10] investigated the existence and local stability of fixed points, fold bifurcation, flip and Neimark-Sacker bifurcations, and complex dynamics of the following discrete-time ratio-dependent predator-prey model with Allee effect:

Liu et al. ^[11] investigated the dynamics behaviors of the following discrete-time predator-prey bioeconomic model:

Liu et al. ^[12] investigated the stability criterion at equilibria, flip and Neimark-Sacker bifurcations, and chaos in the following host-parasitoid model with a lower bound for the host:

Rana ^[13] investigated the existence of fixed points and their stability, flip and Neimark-Sacker bifurcations, and chaos control in the following discrete-time ratio-dependent predator-prey model:

Santra et al. ^[14] investigated the existence of fixed points and their stability, Neimark-Sacker and flip bifurcations, existence of marottos chaos and chaos control in the following discrete-time predator-prey model with crowley-martin functional response incorporating proportional prey refuge:

Mareno and English ^[15] investigated the stability of fixed points, flip and Neimark-Sacker bifurcations of the following coupled logistic map system:

Khan ^[16] investigated the local dynamics and bifurcation analysis of the following discrete-time modified Nicholson-Bailey model:

In 2020, Znegui et al. ^[17] investigated the explicit expression of the Poincaré map for the passive dynamic walking of the compass-gait biped model. Motivated from aforementioned studies, we explore local dynamics, flip bifurcation, chaos control and existence of periodic point of the following predator-prey model with Allee effect on prey population in the interior of $\mathbb{R}^*{_+^2}$ ^[18]:

where $r$ and $a$ are positive constants, and $x_t$ and $y_t$ can be interpreted as the densities of prey and predator populations at time $t$, respectively. It is pointed here that in ^[18], Çelik and Duman have explored local dynamics about equilibria of the model (1.15) and also presented some numerical simulation for correctness of obtained results. Later, Khan et al. ^[19] have explored local dynamical properties about equilibria: $A(1, 0), \ O(0, 0)$ and $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, and Neimark-Sacker bifurcation about the unique positive equilibrium point: $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of the model (1.15). The purpose of present study is to investigate local dynamics, flip bifurcation, chaos control and existence of periodic point of the model (1.15). Our contribution in this article is as follows:

$(1)$ Topological classifications about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of the model.

$(2)$ Detailed analysis of the flip bifurcation about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ by bifurcation theory.

$(3)$ Verification of the obtained results by numerical simulation.

$(4)$ Study of the existence of periodic point of the model.

This article is organized as follows: In section 2, local dynamics and existence of bifurcation about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ are explored. In section 3, we present flip bifurcation about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$. Theoretical results are verified numerically in section 4, and the study of fractal dimension also include in this section. In section 5, chaos control is explored by feedback control method. The existence of period point for the model is explored in section 6 whereas concluding remarks are given in section 7.

2.   Local dynamics and existence of bifurcation about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of the model (1.15)

In the interior of $\mathbb{R}^*{_+^2}$, the result regarding occurrence of fixed point can be stated as a following lemma:

Lemma 2.1. $\forall \; a, r$, model (1.15) has fixed point: $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, which is positive and unique.

Proof. If fixed point of (1.15) is $P_{xy}(x, y)$ then

By straightforward computation, from (2.1) one gets: $x = y = \frac{r}{a+r}$. Therefore $\forall \; a, r$, model (1.15) has a fixed point: $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, which is unique and positive.

Hereafter in order to investigate local dynamics about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, first we present the linearized form of (1.15) about $P_{xy}(x, y)$. So, linearized form of (1.15) under the transformation: $(f, g)\rightarrow (x_{t+1}, y_{t+1})$ is

where

Now $J|_{P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)}$ about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is

And the corresponding characteristic equation of (2.4) is

where

Finally roots of (2.5) are

where

So, we have the following Lemma regarding local dynamics about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of model (1.15) if $\Delta\ge0$.

Lemma 2.2. If $\Delta = (2-r)^2-4\left(\frac{a+r-ar-r^2+ar^2}{a+r}\right)\ge0$ then for $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, the following classifications hold:

(i) $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is stable node if

(ii) $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is unstable node if

(iii) $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is non-hyperbolic if

Proof. $(i)$ If $\Delta = (2-r)^2-4\left(\frac{a+r-ar-r^2+ar^2}{a+r}\right)\ge0$ then from (2.7) real roots of the characteristic equation of $J|_{P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)}$ about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ are $\lambda_{1, 2} = \frac{r-2}{2}\pm\frac{1}{2}\sqrt{(2-r)^2-4\left(\frac{a+r-ar-r^2+ar^2}{a+r}\right)}$. Therefore, by Lemma 2.2 of ^[2], $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of the model (1.15) is stable if $|\lambda_{1, 2}| = \left|\frac{r-2}{2}\pm\frac{1}{2}\sqrt{(2-r)^2-4\left(\frac{a+r-ar-r^2+ar^2}{a+r}\right)}\right| < 1$, which gives $r < \frac{a-2-\sqrt{4+4a-3a^2}}{a-2}$. Therefore $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is a stable node if $r < \frac{a-2-\sqrt{4+4a-3a^2}}{a-2}$. Similarly by Lemma 2.2 of ^[2], it is easy to prove that $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is unstable (respectively non-hyperbolic) if (2.10) (respectively (2.11)) hold.

It is noted that if (2.11) holds then $\lambda_1|_{(2.11)} = -1$ but $\lambda_2|_{(2.11)}$ is neither $1$ nor $-1$ by computation. Hence flip bifurcation exists if $(a, r)$ goes through the following curve:

Remark 1: It is necessary to mention here that results regarding local dynamics along with topological classifications and Neimark-Sacker bifurcation are already published by Khan et al. ^[19] if $\Delta = (2-r)^2-4\left(\frac{a+r-ar-r^2+ar^2}{a+r}\right) < 0$.

3.   Flip bifurcation analysis

Here in this section, flip bifurcation about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is explored. It is noted that flip bifurcation take place if $(a, r)$ goes through the curve, which is depicted in (2.12). Hence considering $r$ in a small neighborhood of $r^*$, $i.e$, $r = r^*+\epsilon$ with $\epsilon < < 1$, then (1.15) becomes

with equilibrium point: $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$. Now one transform $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ into $P_{00}\left(0, 0\right)$ by the transformation:

where $x^* = \frac{r}{a+r}, \ y^* = \frac{r}{a+r}$. Utilizing (3.2) into (3.1) one gets:

where

Now (3.3) becomes

where

by

where

with

Now the center manifold $M^cP_{00}(0, 0)$ about $P_{00}(0, 0)$ for (3.5) is explored in a small neighborhood of $\epsilon$ where

After computation, one gets

Finally, map (3.5) restrict to $M^cP_{00}(0, 0)$ as follows

where

So, in order for (3.12) undergoes flip bifurcation, following holds ^[20,21]:

After calculating, one has

and

with $\alpha_{11} = \sqrt{r^2-2ar-3a^2}, \ \alpha_{12} = 2a+2r-ar-r^2+r\alpha_{11}$. So one has the following following theorem based on above analysis:

Theorem 3.1. About $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, model (1.15) undergoes flip bifurcation if $\epsilon$ varies in the neighborhood of $P_{00}(0, 0)$. Likewise, if $\Gamma_2 < 0 \; (\mathrm{resp.}\ \Gamma_2 > 0)$, then the period-$2$ points bifurcating from $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ are unstable $(\mathrm{resp.} \ \mathrm{stable})$.

4.   Numerical simulations

Here some simulation will be presented in order to justify correctness of obtained results in section 2 & 3. For instance, let $a = 0.95$ then from equation (2.11) one gets: $r = r^*$ where $r^* = 3.149197$. By lemma 2.2, $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is stable node if $r < r^*$, change the behavior if $r = r^*$, unstable node if $r > r^*$ and hence flip bifurcation take place if $r > r^*$. Let $r = 1.5 < r^*$ then from Figure 1a, unique positive fixed point: $P^+_{(0.612245, 0.612245)}(0.612245, 0.612245)$ of (1.15) is stable node. Similarly for rest of bifurcation values, $i.e.$, $r = 1.524, 2.12, 2.5, 2.81, 2.99$ the corresponding fixed point is also stable node (see Figure 1b-1f) and hence simulations agree with the obtained results in section 2. But if $r = 3.2 > 3.149197$ then $\lambda_1 = -1$ but $\lambda_2 = -0.149197\not = 1 \ {\mathrm or} \ -1$. Moreover if $a = 0.95$ and $r = 3.2$ then from (3.15) and (3.16) one gets: $\Gamma_1 = -0.943672$ and $\Gamma_2 = 2933.6713945 > 0$, which implies that stable period-$2$ points bifurcate from $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$. The $2D$ bifurcation diagrams with corresponding maximum Lyapunov exponents (M.L.E) are plotted in Figure 2a-2c. Finally, the trajectories associated with Figure 2a-2b are also plotted in Figure 3a-3f that indicates (1.15) exhibits complex dynamics having orbits of period-$2, 8, 11, 17, 20$ and $22$.

4.1.   Fractal dimension

This designated strange attractors of the system, which is represented as follows ^[22,23]:

with lyapunov exponents are $\lambda_1, \lambda_2, \cdots, \lambda_n$ and largest integer $\jmath$ s.t. $\sum_{i = 1}^{\jmath}\lambda_\jmath\ge0$ and $\sum_{i = 1}^{\jmath+1}\lambda_\jmath < 0$. For the model under consideration (1.15), the fractal dimension becomes of the following form:

Hereafter two Lyapunov exponents are numerically computed for parametric values of $a$ and $r$. If $a = 1.3$ then $\lambda_1 = 2.564421$ (resp. $\lambda_1 = 3.273967$) and $\lambda_2 = -0.584036$ (resp. $\lambda_2 = -0.977670$) for $r = 3.1$ (resp. $\alpha_1 = 3.5$). So fractal dimension for the model (1.15) are

For above chosen parametric values, strange attractors are also plotted in Figure 4a-4b that demonstrate (1.15) has a complex dynamical behavior if $r$ increases.

5.   Chaos control

By the method of state feedback control, here chaos control is explored motivated from existing literature ^[24,25]. By adding $u_t$ as a control force, model (1.15) becomes

where $l_i {(i = 1, 2)}$ represent feedback gains and $x^* = \frac{r}{a+r} = y^*$. The $JC|_{P^+_{xy}(x^*, y^*)}$ of (5.1) is

where

Not its auxiliary equation is

where

If $\lambda_{1, 2}$ are roots of (5.4) then

and

Now it is noted here that lines of marginal stability determine from solution of $\lambda_1 = \pm 1$ as well as $\lambda_1\lambda_2 = 1$, which gives $|\lambda_{1, 2}| < 1$. If $\lambda_1\lambda_2 = 1$ then from (5.7), one gets:

If $\lambda_1 = 1$ then from (5.6) and (5.7) one gets:

Finally if $\lambda_1 = -1$ then from (5.6) and (5.7) one gets:

Hence $L_1, \ L_2$ and $L_3$ in $(l_1, l_2)$- plane delimit the triangular region, which give $|\lambda_{1, 2}| < 1$ (see Figure 5a). Finally, Figures 5b-5c also presented which shows that about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ the chaotic trajectories is stabilized.

6.   Existence of periodic point

Existence of periodic point for (1.15) is explored in this section.

Theorem 6.1. $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of the model (1.15) is periodic point with prime period-$1$.

Proof. From (1.15),

where

After straightforward manipulation, one gets:

From (6.3) this complete the proof.

Remark 2: In similar way, it is easy to establish that $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ of (1.15) is a periodic point of period-$2, 3, \cdots, n$.

7.   Concluding remarks

In this work, we have explored the local dynamics, flip bifurcation, chaos control and existence of periodic points of the predator-prey model with Allee effect. We have proved that $\forall \ a, \ r$ model has a fixed point: $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, which is unique and positive. Further about $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$, local dynamics with topological classifications have investigated if $\Delta = (2-r)^2-4\left(\frac{a+r-ar-r^2+ar^2}{a+r}\right)\ge0$. In particular, we have proved that $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is stable and unstable node if respectively (2.9) and (2.10) hold. Moreover if (2.11) hold then $P^+_{xy}\left(\frac{r}{a+r}, \frac{r}{a+r}\right)$ is non-hyperbolic. Further, it is investigated that flip bifurcation occurs if $(a, r)$ passes through curve, that is represented in (2.12). Some numerical simulations not only exhibit our results with the theoretical analysis but also show the complex dynamical behaviors, such as the period-$2, 8, 11, 17, 20$ and $22$ orbits. Further, in order to show the presence of chaotic behavior in the model, M.L.E and fractal dimension are also computed numerically. By applying method of state feedback control, chaos control is also explored.

Acknowledgments

A. Q. Khan, I. Ahmed and A. Khaliq research is partially supported by the Higher Education Commission of Pakistan while work of H. S. Alayachi and M. S. M. Noorani is financial support by UKM Grant DIP-2017-011 and Ministry of Education Malaysia Grant Number FRGS/1/2017/STG06/UKM/01/1.

Conflict of interest

The author declares that he has no conflicts of interest regarding the publication of this paper.