Bifurcation Analysis of a Discrete-Time Two-Species Model

Abstract

We study the local dynamics and bifurcation analysis of a discrete-time modified Nicholson–Bailey model in the closed first quadrant ${\mathbb{R}}_{+}^2$. It is proved that model has two boundary equilibria: $O\left(0,0\right),A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$, and a unique positive equilibrium $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ under certain parametric conditions. We study the local dynamics along their topological types by imposing method of Linearization. It is proved that fold bifurcation occurs about the boundary equilibria: $O\left(0,0\right),A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$. It is also proved that model undergoes a Neimark–Sacker bifurcation in a small neighborhood of the unique positive equilibrium $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ and meanwhile stable invariant closed curve appears. From the viewpoint of biology, the stable closed curve corresponds to the period or quasi-periodic oscillations between host and parasitoid populations. Some simulations are presented to verify theoretical results. Finally, bifurcation diagrams and corresponding maximum Lyapunov exponents are presented for the under consideration model.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The usual framework for the discrete-time host–parasite models is: $$\begin{array}{}\text{(1)}& \left.\begin{array}{c}{X}_{n+1}=b{X}_{n}f\left(X_n,Y_n\right)\hfill \\ Y_{n+1}=c{X}_n\left(1-f\left(X_n,Y_n\right)\right)\hfill \end{array}\right\},\end{array}$$

where $X_n$ and $Y_n$ represent the population size of the host and parasite in successive generations $n$ and $n+1$ respectively. The parameter $b$ is the host finite rate of increase in the absence of parasites, $c$ is the biomass conversion constant and $f$ is the function defining the fractional survival of hosts from parasitism. The simplest version of this model is that of Nicholson, and Nicholson and Bailey who explored in depth a model in which the proportion of hosts escaping parasitism is given by the zero term of the Poisson distribution [1–3]: $$\begin{array}{}\text{(2)}& f\left(X_n,Y_n\right)=e^{-a{Y}_n},\end{array}$$

where $a{Y}_n$ are the mean encounters per host. Thus, $1-e^{-a{Y}_n}$ is the probability of a host will be attacked. Using (2) in (1), one gets $$\begin{array}{}\text{(3)}& \left.\begin{array}{c}{X}_{n+1}=b{X}_n{e}^{-a{Y}_n}\hfill \\ Y_{n+1}=c{X}_n\left(1-e^{-a{Y}_n}\right)\hfill \end{array}\right\}.\end{array}$$

In 2014, Qureshi et al. [4] have investigated the asymptotic behavior of the following Nicholson–Bailey model: $$\begin{array}{}\text{(4)}& \left.\begin{array}{c}{X}_{n+1}=R{x}_n{e}^{-a\sqrt{y_n}}\hfill \\ Y_{n+1}=x_n\left(1-e^{-a\sqrt{y_n}}\right)\hfill \end{array}\right\},\end{array}$$

where $R$ , $a$ and initial conditions $x_0,y_0$ are positive real numbers. Further in 2015, Khan and Qureshi [5] have investigated the dynamics of the following modified Nicholson–Bailey model: $$\begin{array}{}\text{(5)}& \left.\begin{array}{c}{X}_{n+1}=\frac{b{X}_n{e}^{-a{Y}_n}}{1+d{X}_n}\hfill \\ Y_{n+1}=c{X}_n\left(1-e^{-a{Y}_n}\right)\hfill \end{array}\right\},\end{array}$$

where $a,b,c,d$ and initial conditions $x_0,y_0$ are positive real numbers. Our aim in this paper is to explore the local dynamics along with topological classification and bifurcation analysis of the model (5). First, we make the following rescaling transformations: $$\begin{array}{}\text{(6)}& X_n=\frac{x_n}{b},Y_n=\frac{y_n}{a},\end{array}$$

then system (5) becomes $$\begin{array}{}\text{(7)}& \left.\begin{array}{c}{x}_{n+1}=\frac{b{x}_n{e}^{-y_n}}{1+\left(d,/,b\right)x_n}\hfill \\ y_{n+1}=\frac{ac}{b}{x}_n\left(1-e^{-y_n}\right)\hfill \end{array}\right\}.\end{array}$$

For simplicity, we assume that $b=ac$ , and then model (7) becomes: $$\begin{array}{}\text{(8)}& \left.\begin{array}{c}{x}_{n+1}=\frac{{\zeta }_1{x}_n{e}^{-y_n}}{1+{\zeta }_2{x}_n}\hfill \\ y_{n+1}=x_n\left(1-e^{-y_n}\right)\hfill \end{array}\right\},\end{array}$$

where ${\zeta }_1=b>0$ and ${\zeta }_2=d/b>0$ .

The rest of the paper is organized as follows: Section 2 deals with the study of existence of equilibria of the model (8). In Section 3, we study the local dynamics and existence of bifurcations about equilibria: $O\left(0,0\right),A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right),B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of the model. Section 4 deals with the study of Neimark–Sacker bifurcation about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of the model (8). Numerical simulations along with discussion are presented in the last Section.

2. Existence of Equilibria of the Discrete-Time Model (8)

In this Section, we study the existence of equilibria of the model (8) in ${\mathbb{R}}_{+}^2$ . The results about the existence of equilibria are summarized as follows:

Lemma 1.

Discrete-time model (8) has at least two boundary equilibria and the unique positive equilibrium point in ${\mathbb{R}}_{+}^2$ . More precisely,

(i)

For all parametric values ${\zeta }_1$ and ${\zeta }_2$ , model (8) has a unique equilibrium point: $O\left(0,0\right)$ ;

(ii)

If ${\zeta }_1>1$ then model has boundary equilibrium point: $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ ;

(iii)

Suppose that $$\begin{array}{}\text{(9)}& z=G\left(\widehat{y}\right)=e^{\widehat{y}}+{\zeta }_2\frac{\widehat{y}{e}^{2\widehat{y}}}{e^{\widehat{y}}-1},\end{array}$$ and $z={\zeta }_1$ and when ${\zeta }_1>{\zeta }_2+1$ , the curve $z=G\left(\widehat{y}\right)$ intersect the line $z={\zeta }_1$ at $r$ , say. If ${\zeta }_1>{\zeta }_2+1$ then there exist a unique $r$ such that $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ has a unique positive equilibrium point of (8).

Proof.

For finding number of equilibria of the model (8), we have to solve the following system of equations: $$\begin{array}{}\text{(10)}& \left.\begin{array}{c}\widehat{x}=\frac{{\zeta }_1\widehat{x}{e}^{-\widehat{y}}}{1+{\zeta }_2\widehat{x}}\hfill \\ \widehat{y}=\widehat{x}\left(1-e^{-\widehat{y}}\right)\hfill \end{array}\right\}.\end{array}$$

(i)

Let $\widehat{x}=0$ , then $1^{st}$ equation of system (10) satisfied identically and from $2^{nd}$ equation we obtain $\widehat{y}=0$ . So system (10) has always equilibrium $O\left(0,0\right)$ for all parameter values ${\zeta }_1,{\zeta }_2>0$ .

(ii)

Let $\widehat{y}=0$ , then $2^{nd}$ equation of (10) satisfied identically and from $1^{st}$ equation we obtain $\widehat{x}=\left({\zeta }_1-1\right)/{\zeta }_2$ . Hence system has boundary equilibrium $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ if ${\zeta }_1>1$ .

(iii)

Now we locate the unique positive equilibrium of (10) in ${\mathbb{R}}_{+}^2$ . For this, let $\widehat{x}\ne 0$ , then (10) becomes

$$\begin{array}{}\text{(11)}& \left.\begin{array}{c}1=\frac{{\zeta }_1{e}^{-\widehat{y}}}{1+{\zeta }_2\widehat{x}}\hfill \\ \widehat{y}=\widehat{x}\left(1-e^{-\widehat{y}}\right)\hfill \end{array}\right\}.\end{array}$$

Now eliminating $\widehat{x}$ from (11), one gets $$\begin{array}{}\text{(12)}& {\zeta }_1=e^{\widehat{y}}+{\zeta }_2\frac{\widehat{y}{e}^{2\widehat{y}}}{e^{\widehat{y}}-1}.\end{array}$$

Denote, $$\begin{array}{}\text{(13)}& G\left(\widehat{y}\right)=e^{\widehat{y}}+{\zeta }_2\frac{\widehat{y}{e}^{2\widehat{\mathrm{y}}}}{e^{\widehat{y}}-1}.\end{array}$$

Then the $\widehat{y}$ -coordinates of positive equilibria of (8) are the corresponding $\widehat{y}$ -coordinates of the point of intersection of $z=G\left(\widehat{y}\right)$ and $z={\zeta }_1$ with $\widehat{y}>0$ . By calculating derivative of $F\left(\widehat{y}\right)$ , one get $$\begin{array}{}\text{(14)}& \begin{array}{cc}{G}^{\prime }\left(\widehat{y}\right)\hfill & =e^{\widehat{y}}+{\zeta }_2{e}^{2\widehat{y}}\frac{\left(e^{\widehat{y}}-1\right)\left(2\widehat{y}+1\right)-\widehat{y}{e}^{\widehat{y}}}{{\left(e^{2\widehat{y}}-1\right)}^2}>0,\forall \widehat{y}>0.\end{array}\end{array}$$

Moreover $$\begin{array}{}\text{(15)}& \underset{\widehat{y}\to {\mathrm{0}}^{+}}{\lim }G\left(\widehat{y}\right)=\underset{\widehat{y}\to {\mathrm{0}}^{+}}{\lim }\left(e^{\widehat{y}}+{\zeta }_{\mathrm{2}}\frac{\widehat{y}{e}^{\mathrm{2}\widehat{y}}}{e^{\widehat{y}}-\mathrm{1}}\right),={\zeta }_{\mathrm{2}}+\mathrm{1}.\end{array}$$

So, if ${\zeta }_1\le {\zeta }_2+1$ , then there exists no intersection point of $z=G\left(\widehat{y}\right)$ and $z={\zeta }_1$ . This implies that model (8) has no positive equilibria if ${\zeta }_1\le {\zeta }_2+1$ . And if ${\zeta }_1>{\zeta }_2+1$ , then there exists a unique point of intersection $\left({\zeta }_1,r\right)$ of $z=G\left(\widehat{y}\right)$ and $z={\zeta }_1$ with $r>0$ (see Figure 1). Therefore, if ${\zeta }_1>{\zeta }_2+1$ then (8) has positive equilibrium point and the positive equilibrium point of (8) is unique. We denote it by $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ where $r$ is the positive solution of (12).

[figures omitted; refer to PDF]

3. Local Dynamics and Existence of Bifurcations about Equilibria: $O\left(0,0\right)$ , $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ , $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of the Model (8)

In this Section, we will study the local dynamics of (8) about $O\left(0,0\right),A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ , and $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ . The Jacobian matrix $J_{\left(\widehat{x},\widehat{y}\right)}$ of (8) about equilibrium $\left(\widehat{x},\widehat{y}\right)$ becomes $$\begin{array}{}\text{(16)}& J_{\left(\widehat{x},\widehat{y}\right)}=\left(\begin{array}{cc}\frac{{\zeta }_1{e}^{-\widehat{y}}}{{\left(1+{\zeta }_2\widehat{x}\right)}^2}& -\frac{{\zeta }_1\widehat{x}{e}^{-\widehat{y}}}{1+{\zeta }_2\widehat{x}}\\ 1-e^{-\widehat{y}}& \widehat{x}{e}^{-\widehat{y}}\end{array}\right).\end{array}$$

And its characteristic equation is $$\begin{array}{}\text{(17)}& {\kappa }^2-p\left(\widehat{x},\widehat{y}\right)\kappa +q\left(\widehat{x},\widehat{y}\right)=0,\end{array}$$

where $$\begin{array}{}\text{(18)}& \begin{array}{c}p\left(\widehat{x},\widehat{y}\right)=\frac{{\zeta }_1{e}^{-\widehat{y}}}{{\left(1+{\zeta }_2\widehat{x}\right)}^2}+\widehat{x}{e}^{-\widehat{y}},\\ q\left(\widehat{x},\widehat{y}\right)=\frac{{\zeta }_1\widehat{x}{e}^{-2\widehat{y}}}{{\left(1+{\zeta }_2\widehat{x}\right)}^2}+\frac{{\zeta }_1\widehat{x}{e}^{-\widehat{y}}\left(1-e^{-\widehat{y}}\right)}{1+{\zeta }_2\widehat{x}}.\end{array}\end{array}$$

Lemma 2.

For equilibrium $O$ , the following holds:

(i)

$O$ is a sink if ${\zeta }_1<1$ ;

(ii)

$O$ is never source;

(iii)

$O$ is a saddle if ${\zeta }_1>1$ ;

(iv)

$O$ is nonhyperbolic if ${\zeta }_1=1$ .

From Lemma 2, we can see that one of the eigenvalues about the equilibrium $O\left(0,0\right)$ is 1. So fold bifurcation may occurs when parameter vary in the small neighborhood of ${\zeta }_1=1$ .

Lemma 3.

For $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ , the following holds:

(i)

$A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ is a sink if ${\zeta }_2>{\zeta }_1-1$ ;

(ii)

$A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ is never source;

(iii)

$A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ is a saddle if ${\zeta }_2<{\zeta }_1-1$ ;

(iv)

$A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ is nonhyperbolic if ${\zeta }_2={\zeta }_1-1$ .

We can easily see that if condition (iv) of Lemma 3 hold then one of the eigenvalues about equilibrium $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ is 1. So fold bifurcation may occur when parameters vary in a small neighborhood of ${\zeta }_2={\zeta }_1-1$ . And we denote the parameters satisfying ${\zeta }_2={\zeta }_1-1$ as $$\begin{array}{}\text{(19)}& F_A=\left\{\left({\zeta }_1,{\zeta }_2\right):{\zeta }_2={\zeta }_1-1,{\zeta }_1,{\zeta }_2>0\right\}.\end{array}$$

Hereafter, we will investigate the local dynamics of (8) about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ by using Lemma 2.2 of [6]. The Jacobian matrix $J_{B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)}$ of linearized system of (8) about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is $$\begin{array}{}\text{(20)}& {\kappa }^2-p\kappa +q=0,\end{array}$$

where $$\begin{array}{}\text{(21)}& \begin{array}{c}p=\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}+\frac{r}{e^r-1},\\ q=\frac{r}{e^r-1+{\zeta }_{2}r{e}^r}+r.\\ \end{array}\end{array}$$

Moreover eigenvalues of $J_{B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)}$ about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is given by $$\begin{array}{}\text{(22)}& {\kappa }_{1,2}=\frac{-p\pm \sqrt{\mathrm{\Delta }}}{2},\end{array}$$

where $$\begin{array}{}\text{(23)}& \begin{array}{c}\mathrm{\Delta }=p^2-4q\\ ={\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}+\frac{r}{e^r-1}\right)}^2-4\left(\frac{l}{e^r-1+{\zeta }_{2}r{e}^r}+r\right)\\ ={\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r.\\ \end{array}\end{array}$$

Hereafter, we will give the topological classification of (8) about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ according to the sign of $\mathrm{\Delta }={\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r$ .

Lemma 4.

For $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ , the following holds:

(i) $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is Locally Asymptotically Node if

$$\begin{array}{}\text{(24)}& {\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r\ge 0,0<{\zeta }_2<\frac{{\left(e^r-1\right)}^2}{e^r\left(r-\left(r+1\right)\left(e^r-1\right)\right)};\end{array}$$

(ii) $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is Unstable Node if

$$\begin{array}{}\text{(25)}& {\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r\ge 0,{\zeta }_2>\frac{{\left(e^r-1\right)}^2}{e^r\left(r-\left(r+1\right)\left(e^r-1\right)\right)};\end{array}$$

(iii) $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is nonhyperbolic if

$$\begin{array}{}\text{(26)}& {\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r\ge 0,{\zeta }_2=\frac{{\left(e^r-1\right)}^2}{e^r\left(r-\left(r+1\right)\left(e^r-1\right)\right)}.\end{array}$$

Lemma 5.

For $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ , following statements holds:

(i) $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is Locally Asymptotically Focus if

$$\begin{array}{}\text{(27)}& {\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r<0,0<{\zeta }_2<\frac{e^r\left(1-r\right)-1}{r{e}^r\left(r-1\right)};\end{array}$$

(ii) $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is Unstable Focus if

$$\begin{array}{}\text{(28)}& {\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r<0,{\zeta }_2>\frac{e^r\left(1-r\right)-1}{r{e}^r\left(r-1\right)};\end{array}$$

(iii) $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ is nonhyperbolic if

$$\begin{array}{}\text{(29)}& {\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r<0,{\zeta }_2=\frac{e^r\left(1-r\right)-1}{r{e}^r\left(r-1\right)}.\end{array}$$

If condition (iii) of Lemma 5 holds then we obtain that eigenvalues of $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ are a pair of conjugate complex numbers with modulus one. So Neimark–Sacker bifurcation exists by the variation of parameter in a small neighborhood of ${\zeta }_2=\left(e^r\left(1-r\right)-1\right)/\left(r{e}^r\left(r-1\right)\right)$ . For simplicity, we denote the parameters satisfying ${\zeta }_2=\left(e^r\left(1-r\right)-1\right)/\left(r{e}^r\left(r-1\right)\right)$ as $$\begin{array}{}\text{(30)}& N_B=\left\{\left({\varsigma }_1,{\varsigma }_2\right):{\left(\frac{e^r-1}{e^r-1+{\varsigma }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r<0,{\varsigma }_2=\frac{e^r\left(1-r\right)-1}{r{e}^r\left(r-1\right)},0<r<r^{\ast }\right\}.\end{array}$$

4. Bifurcation Analysis about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of the Model (8)

This Section deals with the study of Neimark–Sacker bifurcation of the model (8) about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ . Consider parameter ${\zeta }_2$ in a small neighborhood of ${\zeta }_2^{\ast }$ , i.e., ${\zeta }_2={\zeta }_2^{\ast }+ϵ$ , where $ϵ\ll 1$ , then (8) becomes: $$\begin{array}{}\text{(31)}& \left.\begin{array}{c}{x}_{n+1}=\frac{{\zeta }_1{x}_n{e}^{-y_n}}{1+\left({\zeta }_2^{\ast }+ϵ\right)x_n}\hfill \\ y_{n+1}=x_n\left(1-e^{-y_n}\right)\hfill \end{array}\right\}.\end{array}$$

The characteristic equation of $J_{B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)}$ about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of (31) is $$\begin{array}{}\text{(32)}& {\kappa }^2-q_1\left(ϵ\right)\kappa +q_2\left(ϵ\right)=0,\end{array}$$

where $$\begin{array}{}\text{(33)}& \begin{array}{c}{q}_1\left(ϵ\right)=\frac{e^r-1}{e^r-1+\left({\zeta }_2^{\ast }+ϵ\right)r{e}^r}+\frac{r}{e^r-1},\\ q_2\left(ϵ\right)=\frac{r}{e^r-1+\left({\zeta }_2^{\ast }+ϵ\right)r{e}^r}+r.\\ \end{array}\end{array}$$

The roots of characteristic equation of $J_{B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)}$ about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ are $$\begin{array}{}\text{(34)}& \begin{array}{c}{\kappa }_{1,2}=\frac{q_1\left(ϵ\right)\pm \iota \sqrt{4q_2\left(ϵ\right)-q_1^2\left(ϵ\right)}}{2}\\ =\frac{e^r-1}{2\left(e^r-1+\left({\zeta }_2^{\ast }+ϵ\right)r{e}^r\right)}+\frac{r}{2\left(e^r-1\right)}\pm \frac{\iota }{2}\sqrt{4r-{\left(\frac{e^r-1}{e^r-1+\left({\zeta }_2^{\ast }+ϵ\right)r{e}^r}-\frac{r}{e^r-1}\right)}^2}.\\ \end{array}\end{array}$$ $$\begin{array}{}\text{(35)}& \left|{\kappa }_{1,2}\right|=\sqrt{q_2\left(ϵ\right)},{\left.\frac{d\left|{\kappa }_{1,2}\right|}{dϵ}\right|}_{ϵ=0}=-\frac{1}{2}{e}^r{\left(r-1\right)}^2<0.\end{array}$$

Additionally, we required that when $ϵ=0$ , ${\kappa }_{1,2}^m\ne 1$ , $m=1,2,3,4$ , which corresponds to $q_1\left(0\right)\ne -2,0,1,2$ . Since $q_1{\left(0\right)}^2-4q_2\left(0\right)<0$ and $q_2\left(0\right)=1$ . Thus $q_1{\left(0\right)}^2<4$ and hence $q_1\left(0\right)\ne \pm 2$ . So we only require that $q_1\left(0\right)\ne 0,1$ . By computation, we get $$\begin{array}{}\text{(36)}& \frac{r}{e^r-1}\ne \frac{\left(e^r-1\right)\left(r-1\right)}{r},\frac{e^r\left(r-1\right)+1}{r}.\end{array}$$

Let $u_n=x_n-x^{\ast },v_n=y_n-y^{\ast }$ then equilibrium $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of system (8) transform into $O\left(0,0\right)$ . By calculating, we obtain $$\begin{array}{}\text{(37)}& \left.\begin{array}{c}{u}_{n+1}=\frac{{\zeta }_1\left(u_n+x^{\ast }\right)e^{-\left(v_n+y^{\ast }\right)}}{1+\left({\zeta }_2^{\ast }+ϵ\right)\left(u_n+x^{\ast }\right)}-x^{\ast }\hfill \\ v_{n+1}=\left(u_n+x^{\ast }\right)\left(1-e^{-\left(v_n+y\ast \right)}\right)-y\ast \hfill \end{array}\right\},\end{array}$$

where $x^{\ast }=\left(r{e}^r\right)/\left(e^r-1\right),y^{\ast }=r$ . Hereafter, when $ϵ=0$ , normal form of system (37) is studied. Expanding (37) up to third order about $\left(u_n,v_n\right)=\left(0,0\right)$ by Taylor series, we get $$\begin{array}{}\text{(38)}& \left.\begin{array}{c}{u}_{n+1}={\mathrm{\Lambda }}_{11}{u}_n+{\mathrm{\Lambda }}_{12}{v}_n+{\mathrm{\Lambda }}_{13}{u}_n^2+{\mathrm{\Lambda }}_{14}{u}_n{v}_n+{\mathrm{\Lambda }}_{15}{v}_n^2+{\mathrm{\Lambda }}_{16}{u}_n^3+{\mathrm{\Lambda }}_{17}{u}_n^2{v}_n+{\mathrm{\Lambda }}_{18}{u}_n{v}_n^2+\hfill \\ {\mathrm{\Lambda }}_{19}{v}_n^3+o\left({\left(\left|u_n\right|+\left|v_n\right|\right)}^3\right)\hfill \\ u_{n+1}={\mathrm{\Lambda }}_{21}{u}_n+{\mathrm{\Lambda }}_{22}{v}_n+{\mathrm{\Lambda }}_{23}{u}_n{v}_n+{\mathrm{\Lambda }}_{24}{v}_n^2+{\mathrm{\Lambda }}_{25}{u}_n{v}_n^2+{\mathrm{\Lambda }}_{26}{v}_n^3+o\left({\left(\left|u_n\right|+\left|v_n\right|\right)}^3\right)\hfill \end{array}\right\},\end{array}$$

where $$\begin{array}{}\text{(39)}& \begin{array}{c}{\mathrm{\Lambda }}_{11}=\frac{1}{1+{\zeta }_2^{\ast }{x}^{\ast }},{\mathrm{\Lambda }}_{12}=-x^{\ast },{\mathrm{\Lambda }}_{13}=-\frac{{\zeta }_2^{\ast }}{{\left(1+{\zeta }_2^{\ast }{x}^{\ast }\right)}^2},\\ {\mathrm{\Lambda }}_{14}=-\frac{1}{1+{\zeta }_2^{\ast }{x}^{\ast }},{\mathrm{\Lambda }}_{15}=\frac{1}{2}{x}^{\ast },{\mathrm{\Lambda }}_{16}=\frac{{\zeta }_2^{\ast }}{{\left(1+{\zeta }_2^{\ast }{x}^{\ast }\right)}^3},\\ {\mathrm{\Lambda }}_{17}=\frac{{\zeta }_2^{\ast }}{{\left(1+{\zeta }_2^{\ast }{x}^{\ast }\right)}^2},{\mathrm{\Lambda }}_{18}=\frac{1}{2\left(1+{\zeta }_2^{\ast }{x}^{\ast }\right)},{\mathrm{\Lambda }}_{19}=-\frac{1}{6}{x}^{\ast },\\ {\mathrm{\Lambda }}_{21}=\frac{e^{y^{\ast }}-1}{e^{y^{\ast }}},{\mathrm{\Lambda }}_{22}=\frac{x^{\ast }}{e^{y^{\ast }}},{\mathrm{\Lambda }}_{23}=\frac{1}{e^{y^{\ast }}},{\mathrm{\Lambda }}_{24}=-\frac{1}{2}\frac{x^{\ast }}{e^{y^{\ast }}},\\ {\mathrm{\Lambda }}_{25}=-\frac{1}{2}\frac{1}{e^{y^{\ast }}},{\mathrm{\Lambda }}_{26}=\frac{1}{6}\frac{x^{\ast }}{e^{y^{\ast }}}.\\ \end{array}\end{array}$$

Now, let $$\begin{array}{}\text{(40)}& \eta =\frac{e^r-1}{2\left(e^r-1+{\zeta }_2^{\ast }r{e}^r\right)}+\frac{r}{2\left(e^r-1\right)},\zeta =\frac{1}{2}\sqrt{4r-{\left(\frac{e^r-1}{e^r-1+{\zeta }_2^{\ast }r{e}^r}-\frac{r}{e^r-1}\right)}^2},\end{array}$$

and invertible matrix $T$ defined by $$\begin{array}{}\text{(41)}& T=\left(\begin{array}{cc}{\mathrm{\Lambda }}_{12}& 0\\ \eta -{\mathrm{\Lambda }}_{11}& -\zeta \end{array}\right).\end{array}$$

Using following translation $$\begin{array}{}\text{(42)}& \left(\begin{array}{c}{u}_n\\ v_n\end{array}\right)=\left(\begin{array}{cc}{\mathrm{\Lambda }}_{12}\hfill & 0\\ \eta -{\mathrm{\Lambda }}_{11}& -\zeta \end{array}\right)\left(\begin{array}{c}{\hat{X}}_n\hfill \\ {\hat{Y}}_n\hfill \end{array}\right),\end{array}$$

(38) gives: $$\begin{array}{}\text{(43)}& \left(\begin{array}{c}\widehat{X_{n+1}}\hfill \\ \widehat{Y_{n+1}}\hfill \end{array}\right)=\left(\begin{array}{cc}\eta \hfill & -\zeta \\ \zeta \hfill & \eta \end{array}\right)\left(\begin{array}{c}\widehat{X_n}\hfill \\ \widehat{Y_n}\hfill \end{array}\right)\left(\begin{array}{c}\hat{\mathrm{\Phi }}\left(\widehat{X_n},\widehat{Y_n}\right)\hfill \\ \hat{\mathrm{\Psi }}\left(\widehat{X_n},\widehat{Y_n}\right)\hfill \end{array}\right),\end{array}$$

where $$\begin{array}{}\text{(44)}& \begin{array}{c}\hat{\mathrm{\Phi }}\left(\widehat{X_n},\widehat{Y_n}\right)={\mathrm{\Omega }}_{11}{\widehat{X_n}}^2+{\mathrm{\Omega }}_{12}\widehat{X_n}\widehat{Y_n}+{\mathrm{\Omega }}_{13}{\widehat{Y_n}}^2+{\mathrm{\Omega }}_{14}{\widehat{X_n}}^3+{\mathrm{\Omega }}_{15}{\widehat{X_n}}^2\widehat{Y_n}+{\mathrm{\Omega }}_{16}\widehat{X_n}{\widehat{Y_n}}^2+\\ {\mathrm{\Omega }}_{17}{\widehat{Y_n}}^3+o\left({\left(\left|\widehat{X_n}\right|+\left|\widehat{Y_n}\right|\right)}^3\right),\\ \hat{\mathrm{\Psi }}\left(\widehat{X_n},\widehat{Y_n}\right)={\mathrm{\Omega }}_{21}{\widehat{X_n}}^2+{\mathrm{\Omega }}_{22}\widehat{X_n}\widehat{Y_n}+{\mathrm{\Omega }}_{23}{\widehat{Y_n}}^2+{\mathrm{\Omega }}_{24}{\widehat{X_n}}^3+{\mathrm{\Omega }}_{25}{\widehat{X_n}}^2\widehat{Y_n}+{\mathrm{\Omega }}_{26}\widehat{X_n}{\widehat{Y_n}}^2+\\ {\mathrm{\Omega }}_{27}{\widehat{Y_n}}^3+o\left({\left(\left|\widehat{X_n}\right|+\left|\widehat{Y_n}\right|\right)}^3\right),\\ \end{array}\end{array}$$

$$\begin{array}{}\text{(45)}& \begin{array}{c}{\mathrm{\Omega }}_{11}={\mathrm{\Lambda }}_{12}{\mathrm{\Lambda }}_{13}+\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left[{\mathrm{\Lambda }}_{14}+\frac{{\mathrm{\Lambda }}_{15}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}\right],{\mathrm{\Omega }}_{12}=-\zeta \left[{\mathrm{\Lambda }}_{14}+\frac{2{\mathrm{\Lambda }}_{15}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}\right],\\ {\mathrm{\Omega }}_{13}=\frac{{\zeta }^2{\mathrm{\Lambda }}_{15}}{{\mathrm{\Lambda }}_{12}},{\mathrm{\Omega }}_{14}={\mathrm{\Lambda }}_{12}\left({\mathrm{\Lambda }}_{12}{\mathrm{\Lambda }}_{16}+{\mathrm{\Lambda }}_{17}\left(\eta -{\mathrm{\Lambda }}_{11}\right)\right)+{\left(\eta -{\mathrm{\Lambda }}_{11}\right)}^2\left({\mathrm{\Lambda }}_{18}+\frac{{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}\right),\\ {\mathrm{\Omega }}_{15}=-\zeta \left[{\mathrm{\Lambda }}_{12}{\mathrm{\Lambda }}_{17}+\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left(2{\mathrm{\Lambda }}_{18}+\frac{3{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}\right)\right],\\ {\mathrm{\Omega }}_{16}={\zeta }^2\left[{\mathrm{\Lambda }}_{18}+\frac{3{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}\right],{\mathrm{\Omega }}_{17}=-\frac{{\zeta }^3{\mathrm{\Lambda }}_{19}}{{\mathrm{\Lambda }}_{12}},\\ {\mathrm{\Omega }}_{21}=\frac{\eta -{\mathrm{\Lambda }}_{11}}{\zeta }\left[{\mathrm{\Lambda }}_{12}\left({\mathrm{\Lambda }}_{13}+{\mathrm{\Lambda }}_{14}\left(\eta -{\mathrm{\Lambda }}_{11}\right)-{\mathrm{\Lambda }}_{23}\right)+\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left(\frac{{\mathrm{\Lambda }}_{15}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{24}\right)\right],\\ {\mathrm{\Omega }}_{22}=2\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left(\frac{{\mathrm{\Lambda }}_{15}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{24}\right)-{\mathrm{\Lambda }}_{12}{\mathrm{\Lambda }}_{14}\left(\eta -{\mathrm{\Lambda }}_{11}-{\mathrm{\Lambda }}_{23}\right),\\ {\mathrm{\Omega }}_{23}=\zeta \left[\frac{{\mathrm{\Lambda }}_{15}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{24}\right],\\ {\mathrm{\Omega }}_{24}=\frac{\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{\zeta }\left[{\mathrm{\Lambda }}_{12}\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left({\mathrm{\Lambda }}_{12}{\mathrm{\Lambda }}_{16}+{\mathrm{\Lambda }}_{17}\left(\eta -{\mathrm{\Lambda }}_{11}\right)+\left(\frac{{\mathrm{\Lambda }}_{18}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{25}\right)\right)+\right.\\ \left.{\left(\eta -{\mathrm{\Lambda }}_{11}\right)}^2\left(\frac{{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{26}\right)\right],\\ {\mathrm{\Omega }}_{25}=-\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left[{\mathrm{\Lambda }}_{12}\left(\frac{{\mathrm{\Lambda }}_{17}}{\zeta }+2\left(\frac{{\mathrm{\Lambda }}_{18}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{25}\right)\right)+\right.\\ \left.3\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left(\frac{{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{26}\right)\right],\\ {\mathrm{\Omega }}_{26}=\zeta \left[{\mathrm{\Lambda }}_{12}\left(\frac{{\mathrm{\Lambda }}_{18}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{25}\right)+3\left(\eta -{\mathrm{\Lambda }}_{11}\right)\left(\frac{{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{26}\right)\right],\\ {\mathrm{\Omega }}_{27}=-{\zeta }^2\left[\frac{{\mathrm{\Lambda }}_{19}\left(\eta -{\mathrm{\Lambda }}_{11}\right)}{{\mathrm{\Lambda }}_{12}}-{\mathrm{\Lambda }}_{26}\right].\\ \end{array}\end{array}$$

In addition, $$\begin{array}{}\text{(46)}& \begin{array}{c}{\left.{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{11},{\left.{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}={\mathrm{\Omega }}_{12},{\left.{\widehat{\mathrm{\Phi }}}_{\widehat{Y_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{13},{\left.{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}\widehat{X_n}}\right|}_{\left(0,0\right)}=6{\mathrm{\Omega }}_{14},\\ {\left.{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{15},{\left.{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{Y_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{16},{\left.{\widehat{\mathrm{\Phi }}}_{\widehat{Y_n}\widehat{Y_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=6{\mathrm{\Omega }}_{17},\\ \end{array}\end{array}$$ $$\begin{array}{}\text{(47)}& \begin{array}{c}{\left.{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{21},{\left.{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}={\mathrm{\Omega }}_{22},{\left.{\widehat{\mathrm{\Psi }}}_{\widehat{Y_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{23},{\left.{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}\widehat{X_n}}\right|}_{\left(0,0\right)}=6{\mathrm{\Omega }}_{24},\\ {\left.{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{25},{\left.{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{Y_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=2{\mathrm{\Omega }}_{26},{\left.{\widehat{\mathrm{\Psi }}}_{\widehat{Y_n}\widehat{Y_n}\widehat{Y_n}}\right|}_{\left(0,0\right)}=6{\mathrm{\Omega }}_{27}.\\ \end{array}\end{array}$$

In order for (43) to undergo Neimark–Sacker bifurcation, it is required that following discriminatory quantity, i.e., $\mathrm{\Psi }\ne 0$ (see [6–13]). $$\begin{array}{}\text{(48)}& \mathrm{\Psi }=-Re\left[\frac{\left(1-2\overline{\kappa }\right){\overline{\kappa }}^2}{1-\kappa }{\tau }_{11}{\tau }_{20}\right]-\frac{1}{2}{∥{\tau }_{11}∥}^2-{∥{\tau }_{02}∥}^2+Re\left(\overline{\kappa }{\tau }_{21}\right),\end{array}$$

where $$\begin{array}{}\text{(49)}& \begin{array}{c}{\tau }_{02}={\left.\frac{1}{8}\left[{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}}-{\widehat{\mathrm{\Phi }}}_{\widehat{Y_n}\widehat{Y_n}}+2{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{Y_n}}+\iota \left({\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}}-{\widehat{\mathrm{\Psi }}}_{\widehat{Y_n}\widehat{Y_n}}+2{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{Y_n}}\right)\right]\right|}_{\left(0,0\right)},\\ {\tau }_{11}={\left.\frac{1}{4}\left[{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}}+{\widehat{\mathrm{\Phi }}}_{\widehat{Y_n}\widehat{Y_n}}+\iota \left({\mathrm{\Psi }}_{\widehat{X_n}\widehat{X_n}}+{\widehat{\mathrm{\Psi }}}_{\widehat{Y_n}\widehat{Y_n}}\right)\right]\right|}_{\left(0,0\right)},\\ {\tau }_{20}={\left.\frac{1}{8}\left[{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}}-{\widehat{\mathrm{\Phi }}}_{\widehat{Y_n}\widehat{Y_n}}+2{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{Y_n}}+\iota \left({\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}}-{\widehat{\mathrm{\Psi }}}_{\widehat{Y_n}\widehat{Y_n}}-2{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{Y_n}}\right)\right]\right|}_{\left(0,0\right)},\\ {\tau }_{21}=\frac{1}{16}\left[{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}\widehat{X_n}}+{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{Y_n}\widehat{Y_n}}+{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}\widehat{Y_n}}+{\widehat{\mathrm{\Psi }}}_{\widehat{Y_n}\widehat{Y_n}\widehat{Y_n}}+\right.\\ {\left.\left.\iota \left({\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{X_n}\widehat{X_n}}+{\widehat{\mathrm{\Psi }}}_{\widehat{X_n}\widehat{Y_n}\widehat{Y_n}}-{\widehat{\mathrm{\Phi }}}_{\widehat{X_n}\widehat{X_n}\widehat{Y_n}}-{\widehat{\mathrm{\Phi }}}_{\widehat{Y_n}\widehat{Y_n}\widehat{Y_n}}\right)\right]\right|}_{\left(0,0\right)}.\\ \end{array}\end{array}$$

After calculating, we get $$\begin{array}{}\text{(50)}& \begin{array}{c}{\tau }_{02}=\frac{1}{4}\left[{\mathrm{\Omega }}_{11}-{\mathrm{\Omega }}_{13}+{\mathrm{\Omega }}_{22}+\iota \left({\mathrm{\Omega }}_{21}-{\mathrm{\Omega }}_{23}+{\mathrm{\Omega }}_{12}\right)\right],\\ {\tau }_{11}=\frac{1}{2}\left[{\mathrm{\Omega }}_{11}+{\mathrm{\Omega }}_{13}+\iota \left({\mathrm{\Omega }}_{21}+{\mathrm{\Omega }}_{23}\right)\right],\\ {\tau }_{20}=\frac{1}{4}\left[{\mathrm{\Omega }}_{11}-{\mathrm{\Omega }}_{13}+{\mathrm{\Omega }}_{22}+\iota \left({\mathrm{\Omega }}_{21}-{\mathrm{\Omega }}_{23}-{\mathrm{\Omega }}_{12}\right)\right],\\ {\tau }_{21}=\frac{1}{8}\left[3{\mathrm{\Omega }}_{14}+{\mathrm{\Omega }}_{16}+{\mathrm{\Omega }}_{25}+3{\mathrm{\Omega }}_{27}+\iota \left(3{\mathrm{\Omega }}_{24}+{\mathrm{\Omega }}_{26}-{\mathrm{\Omega }}_{15}-3{\mathrm{\Omega }}_{17}\right)\right].\\ \end{array}\end{array}$$

Based on this analysis and Neimark–Sacker bifurcation Theorem discussed in [12, 13], we arrive at the following Theorem:

Theorem 1.

If $\mathrm{\Psi }\ne 0$ then model (8) undergoes a Neimark–Sacker bifurcation about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ as the parameters $\left({\zeta }_1,{\zeta }_2\right)$ go through $N_B$ . Additionally, attracting (respectively repelling) invariant closed curve bifurcate from $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ if $\mathrm{\Psi }<0$ (respectively $\mathrm{\Psi }>0$ ).

According to Neimark–Sacker bifurcation discussed in [12, 13], the bifurcation is called supercritical Neimark–Sacker bifurcation if the discriminatory quantity $\mathrm{\Psi }<0$ . In the following Section, numerical simulations guarantee that supercritical Neimark–Sacker bifurcation occurs for the model (8). Biologically, attracting closed curve indicates that both parasitoid and host populations will coexist under the periodic or quasi-periodic oscillations with long time.

5. Numerical Simulations and Discussion

This work deals with the study of local dynamics and bifurcation analysis of a discrete-time two-species model in ${\mathbb{R}}_{+}^2$ . We proved that the model has two boundary equilibria: $O,A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ and a unique positive equilibrium point $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ under certain parametric conditions. We studied the local dynamics along with topological classification about equilibria: $O\left(0,0\right)$ , $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ , $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right),$ and conclusion is presented in Table 1. We proved that about $A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$ there may exist a fold bifurcation when parameters of (8) are located in the set: $F_A=\left\{\left({\zeta }_1,{\zeta }_2\right):{\zeta }_2={\zeta }_1-1,{\zeta }_1,{\zeta }_2>0\right\}$ . We also proved that if ${\zeta }_2=\left(e^r\left(1-r\right)-1\right)/\left(r{e}^r\left(r-1\right)\right)$ then eigenvalues $J_{B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)}$ about $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ are pair of complex conjugate with modulus one and thus in particular supercritical Neimark–Sacker bifurcation occurs under the bifurcation curve: $$\begin{array}{}\text{(51)}& N_B=\left\{\left({\zeta }_1,{\zeta }_2\right):{\left(\frac{e^r-1}{e^r-1+{\zeta }_{2}r{e}^r}-\frac{r}{e^r-1}\right)}^2-4r<0,{\zeta }_2=\frac{e^r\left(1-r\right)-1}{r{e}^r\left(r-1\right)},0<r<r^{\ast }\right\}.\end{array}$$

Table 1

Number of equilibria along their qualitative behavior of (8).

E.P	Corresponding behavior
$O\left(0,0\right)$	Sink if ${\zeta }_1<1$ ; never source; saddle if ${\zeta }_1>1$ ; nonhyperbolic if ${\zeta }_1=1$ .
$A\left(\left({\zeta }_1-1\right)/{\zeta }_2,0\right)$	Sink if ${\zeta }_2>{\zeta }_1-1$ ; never source; saddle if ${\zeta }_2<{\zeta }_1-1$ ; nonhyperbolic if ${\zeta }_2={\zeta }_1-1$ .
$B\left(r{e}^r/\left(e^r-1\right),r\right)$	Locally asymptotically node if
	${\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r\ge 0\text{and}0<{\zeta }_2<\left({(e^r-1)}^2\right)/\left(e^r\left(r-\left(r+1\right)\left(e^r-1\right)\right)\right);$
	Unstable node if
	${\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r\ge 0\text{and}{\zeta }_2>{(e^r-1)}^2/e^r\left(r-\left(r+1\right)\left(e^r-1\right)\right);$
	Nonhyperbolic (for real eigenvalues) if
	${\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r\ge 0\text{and}{\zeta }_2=\left({(e^r-1)}^2\right)/\left(e^r\left(r-\left(r+1\right)\left(e^r-1\right)\right)\right);$
	Locally asymptotically focus if
	${\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r<0\text{and}0<{\zeta }_2<\left(e^r\left(1-r\right)-1\right)/\left(r{e}^r\left(r-1\right)\right);$
	Unstable focus if
	${\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r<0\text{and}{\zeta }_2>\left(e^r\left(1-r\right)-1\right)/\left(r{e}^r\left(r-1\right)\right);$
	Nonhyperbolic (for complex eigenvalues) if
	${\left(\left(e^r-1\right)/\left(e^r-1+{\zeta }_{2}r{e}^r\right)-\left(r\right)/\left(e^r-1\right)\right)}^2-4r<0\text{and}{\zeta }_2=\left(e^r\left(1-r\right)-1\right)/\left(r{e}^r\left(r-1\right)\right).$

Biologically, existence of stable closed curves implies that there exist the periodic or quasiperiodic oscillations between host and parasitoid populations. Finally, numerical simulations are provided to verify theoretical discussion. These numerical simulations presented in Figures 2–5 agree with our theoretical discussion. Figure 2 shows that $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of the model (8) is Locally Asymptotically Focus when ${\zeta }_2<{\zeta }_2^{\ast }$ , where ${\zeta }_2^{\ast }\approx 0.5000046770$ as presented in Figures 2(a)–2(i) by choosing ${\zeta }_1=2.9899$ . But when ${\zeta }_2$ goes through the bifurcation value ${\zeta }_2^{\ast }$ , equilibrium $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of (8) is Unstable Focus. Meanwhile, an attracting closed invariant curve bifurcates from $B\left(\left(r{e}^r\right)/\left(e^r-1\right),r\right)$ of the model (8) as presented in Figures 3(a)–3(i). Moreover, bifurcation diagrams along with Maximum Lyapunov Exponent in this case, are plotted and drawn in Figure 4. Finally 3D bifurcation diagrams are also plotted and drawn in Figure 5.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Disclosure

The author declares that he got no funding on any part of this research.

Conflicts of Interest

The author declares that they have no conflicts of interest.

Acknowledgments

This research is partially supported by the Higher Education Commission of Pakistan.

References

[1] A. J. Nicholson, "The balance of animal populations," The Journal of Animal Ecology, vol. 2 no. 1, pp. 131-178, DOI: 10.2307/954, 1933.

[2] A. Nicholson, V. Bailey, "The balance of animal population," Proceedings of the Zoological Society of London, vol. 105 no. 3, pp. 551-598, DOI: 10.1111/j.1096-3642.1935.tb01680.x, 1935.

[3] Y. Kang, D. Armbruster, Y. Kuang, "Dynamics of a plant-herbivore model with applications to gypsy moth outbreaks," Journal of Biological Dynamics, 2007.

[4] M. N. Qureshi, A. Q. Khan, Q. Din, "Asymptotic behavior of a Nicholson–Bailey model," Advances in Difference Equations, vol. 1,DOI: 10.1186/1687-1847-2014-62, 2014.

[5] A. Q. Khan, M. N. Qureshi, "Dynamics of a modified Nicholson–Bailey host-parasitoid model," Advances in Difference Equations, vol. 2015,DOI: 10.1186/s13662-015-0357-2, 2015.

[6] X. Liu, D. Xiao, "Complex dynamic behaviors of a discrete-time predator–prey system," Chaos, Solitons & Fractals, vol. 32 no. 1, pp. 80-94, DOI: 10.1016/j.chaos.2005.10.081, 2007.

[7] A. Q. Khan, M. N. Qureshi, "Global dynamics and bifurcations analysis of a two-dimensional discrete-time Lotka-Volterra model," Complexity,DOI: 10.1155/2018/7101505, 2018.

[8] Z. Hu, Z. Teng, L. Zhang, "Stability and bifurcation analysis of a discrete predator–prey model with nonmonotonic functional response," Nonlinear Analysis: Real World Applications, vol. 12 no. 4, pp. 2356-2377, DOI: 10.1016/j.nonrwa.2011.02.009, 2011.

[9] Z. Jing, J. Yang, "Bifurcation and chaos in discrete-time predator–prey system," Chaos, Solitons & Fractals, vol. 27 no. 1, pp. 259-277, DOI: 10.1016/j.chaos.2005.03.040, 2006.

[10] C.-H. Zhang, X.-P. Yan, G.-H. Cui, "Hopf bifurcations in a predator–prey system with a discrete delay and a distributed delay," Nonlinear Analysis: Real World Applications, vol. 11 no. 5, pp. 4141-4153, DOI: 10.1016/j.nonrwa.2010.05.001, 2010.

[11] M. Sen, M. Banerjee, A. Morozov, "Bifurcation analysis of a ratio-dependent prey-predator model with the Allee effect," Ecological Complexity, vol. 11, pp. 12-27, DOI: 10.1016/j.ecocom.2012.01.002, 2012.

[12] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields,DOI: 10.1007/978-1-4612-1140-2, 1983.

[13] Y. A. Kuznetson, Elements of Applied Bifurcation Theory, 2004.