Complex Dynamic Behaviors of a Modified Discrete Leslie–Gower Predator–Prey System with Fear Effect on Prey Species

Abstract

A discrete modified Leslie–Gower prey-predator model considering the effect of fear on prey species is proposed and studied in this paper. First, we discuss the existence of equilibria and the local stability of the model. Second, we use the iterative method and comparison principle to obtain the set of conditions which ensures the global attractivity of positive equilibrium point. The results show that prey and predator can coexist stably when the intrinsic growth rates of both prey and predator are maintained within a certain range. Then, we study the global attractivity of the boundary equilibrium point. Our results suggest that when the intrinsic rate of prey is small enough or the fear factor is large enough, the prey will tend to go extinct, while the predator can survive stably due to the availability of other food sources. Subsequently, we discuss flip bifurcation, transcritical bifurcation at the equilibria of the system, by using the center manifold theorem and bifurcation theory. We find that system changes from chaotic state to four-period orbit, two-period orbit, stable state, and finally prey species will be driven to extinction, while predator species survive in a stable state for enough large birth rate of prey species with the increasing of fear effect. Finally, we verify the feasibility of the main results by numerical simulations, and discuss the influence of the fear effect. The results show that the fear effect within a certain range can enhance the stability of the system.

1. Introduction

The dynamic quantitative relationship between predator and prey has always been a hot topic in ecology and mathematical ecology. The derivation of more realistic predator–prey dynamics models is a major trend of relevant theoretical research. Traditionally, all predator–prey models evolve from two basic models—the first is the famous Lotka–Volterra predator–prey model [1], which takes the form

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{dt}}& =& x(a-bx)-cxy,\hfill \\ \hfill {\displaystyle \frac{dy}{dt}}& =& -dy+k_{1}cxy.\hfill \end{array}$$

(1)

In this model, $k_1$ is the transform rate of the predator species. A total of $cxy$ prey populations are captured, and the energy gained by the predator will further generate $k_{1}cxy$ individuals in the predator population. The model can be adjusted to a more suitable predator–prey model by incorporating functional response, which takes the form

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{dt}}& =& x(a-bx)-c\Psi (x,y)y,\hfill \\ \hfill {\displaystyle \frac{dy}{dt}}& =& -dy+k_{1}c\Psi (x,y)y.\hfill \end{array}$$

(2)

There are numerous kinds of functional response, depending on the circumstance considered. These include the Holling II functional response [1], Beddington–DeAngelis functional response [2] and square root functional response [3], etc.

The second one is now referred to as the Leslie–Gower predator–prey model, which was first introduced by Leslie [4], where the carrying capacity of the predator’s environment is proportional to the number of prey:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{dt}}& =& x(r_1-a_{1}x-b_{1}y),\hfill \\ \hfill {\displaystyle \frac{dy}{dt}}& =& {\displaystyle y\left(r_2-\frac{a_{2}y}{x}\right).}\hfill \end{array}$$

(3)

The model seems very simple, however, only in 2001 did Korobeinikov [5] give a strict proof of the global stability property of the positive equilibrium. Since then, numerous studies had been conducted in this direction. For example, Chen, Chen, and Xie [6] incorporated prey refuge in the system (3) and showed that refuge has a complicated influence on the dynamic behaviors of the predator species. To investigate the influence of human disturbance, Chen and Chen [7] proposed a Leslie–Gower predator–prey model with feedback control; by constructing a suitable Lyapunov function, they obtained a set of sufficient conditions which ensured the existence of a unique globally asymptotically stable positive equilibrium of the system. Li, Han, and Chen [8] studied the stability property of the stage-structured predator–prey model. Recently, several scholars [9,10,11,12] studied the influence of Allee effect on the Leslie–Gower predator–prey system. They gave a detailed bifurcation analysis of the model proposed. Yu and Chen [13] and Yu [14] studied the mutual interferences of the predator species. Zhu and Kong [15] showed that nonlinear harvesting may lead to very complicated dynamic behaviors in the system. Zou and Guo [16] studied the influence of a spatially heterogeneous environment on a diffusive Leslie–Gower predator–prey model. Mondal, Pal, and Samanta [17] proposed a Leslie–Gower predator–prey eco-epidemiological model with disease in predator; they gave a thorough analysis of the dynamic behaviors of the system. Liang, Zeng, Pang, and Liang [18] showed that a Leslie–Gower predator–prey system with ratio-dependent and state impulsive feedback control could have a periodic solution, which is quite different from the dynamic behaviors of the system (3).

Aziz-Alaoui and Okiye [19] proposed a modified Leslie–Gower predator–prey model in 2003. The authors argued that in the case of severe scarcity, a predator species can switch over to other populations, but its growth will be limited by the fact that its favorite food (prey species) is not available in abundance. The modified Leslie–Gower predator–prey model with Holling I functional response can be described by the following equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}& =& x(r_1-a_{1}x-b_{1}y),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}& =& y\left(r-{\displaystyle \frac{fy}{m+x}}\right),\hfill \end{array}$$

(4)

where $x\left(t\right)$ and $y\left(t\right)$ represent the population density of the prey population and the predator population at time t, respectively; m represents the environmental protection rate for predators, which reflects the fact that the predators have other food resources; r represents the birth rate of the predator population; and f represents the mortality of the predator population due to intraspecific competition. Wu [20] further incorporated prey refuge to above system, and he showed that under some suitable assumption, the system may exist a unique positive equilibrium. The condition which ensure the existence of the positive equilibrium is enough to ensure its globally asymptotically stability.

Considering that the population of predators can not only directly reduce the population density of the prey species by catching prey as its food resource, the presence of predators creates fear in the prey species because the prey species is always alert to potential attacks, and the prey will engage in reverse feeding behavior in response to the perceived risk (e.g., foraging frequency decreases and decreased number of offspring foraging), which leads to a decrease in the prey species population density. The fear response of prey is called the fear effect. The first famous experiment was done by Zanette, White, Allen, and Clinchy [21]. They experimented on song sparrows over the course of an entire breeding season to see if perceived predation danger could, even in the absence of direct killing, influence reproduction. They discovered that having a fear of predators caused the number of offspring to reduce by 40%.

In 2016, based on system (1), Wang, Zanette and Zou [22] first proposed a continuous predator–prey model incorporating the fear effect:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}& =& r_{0}uf(k,v)-du-a{u}^2-puv,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}}& =& cpuv-mv,\hfill \end{array}$$

(5)

where $u\left(t\right)$ and $v\left(t\right)$ represent the population density of the prey population and the predator population at time t, respectively. k represents the level of fear that drives the anti-predatory behavior of prey, $r_0$ and d represent the birth rate and natural mortality rate of the prey population, respectively, a represents the mortality of the prey population due to intraspecific competition, and b represents maximum capture rate of predators. The authors showed that for system (5), fear effect had no influence on the existence or stability of the equilibria.

Zhu, Wu, Lai et al. [23] proposed the following predator–prey model:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}& =& {\displaystyle \frac{r_{0}u}{1+kv}-du-a{u}^2-puv,}\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}}& =& cpuv+mv-d_1{v}^2.\hfill \end{array}$$

(6)

In this model, without the prey species, the predator species takes the form

$${\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}=mv-d_1{v}^2.}$$

(7)

It is a logistic equation, which means that in system (6), the authors made the assumption that the predator species takes other species as food resources. Their study showed that the fear effect is one of the most important factors leading to the extinction of the prey species. This result is quite different from the results of Wang, Zanette and Zou [22].

It is natural to study the influence of the fear effect in the Leslie–Gower-type predator–prey system; indeed, many scholars [9,10,24,25,26] have conducted research in this direction. Firdiansyah [24] employed a Leslie–Gower predator–prey model with Beddington–DeAngelis functional response to examine the influence of fear effect; the model had the following form:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}& =& {\displaystyle \frac{r_{1}x}{1+Ky}-bx-p{x}^2-\frac{\alpha (1-m)xy}{a+b(1-m)x+cy},}\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}& =& {\displaystyle y\left(r_2-\frac{\beta y}{(1-m)x+\gamma }\right).}\hfill \end{array}$$

(8)

The author demonstrated that an increase in fear might reduce the population density of both species. However, in the event of a constant fear rate, prey refuge is beneficial to the survival of both species. For more works on continuous predator–prey models incorporating the fear effect, one could refer to [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38] and the references cited therein.

It is well-known that when species have non-overlapping generations or their population sizes are too small, discrete models described by difference equations are more appropriate than continuous-time ones. Generally speaking, discrete systems will have more complex dynamical behaviors. However, to this day, there are not many works on discrete predator–prey systems incorporating the fear effect.

Corresponding to system (5), Kundu, Pal, and Samanta [34] proposed the following discrete predator–prey system with consideration of fear effect in prey species:

$$\begin{array}{ccc}\hfill x(n+1)& =& {\displaystyle x\left(n\right)exp[\frac{r_0}{1+ky\left(n\right)}-d-ax\left(n\right)-py\left(n\right)],}\hfill \\ \hfill y(n+1)& =& y\left(n\right)exp\left[cpx\right(n)-m].\hfill \end{array}$$

(9)

The dynamic behaviors of this system are quite different from those of the system (5). Indeed, the authors showed that the stability of the coexistence equilibrium point of the corresponding discrete-time model changed (unstable to stable) with the increase of the cost of fear k.

Stimulated by the work of Wang, Zanette, and Zou [22]; Zhu, Wu, Lai et al. [23]; and Kundu, Pal, and Samanta [34], corresponding to system (6), Chen, He, and Chen [35] recently proposed and studied a discrete predator–prey model, in which the prey population has a fear effect and the predator population has other food sources:

$$\begin{array}{ccc}\hfill x(n+1)& =& {\displaystyle x\left(n\right)exp[\frac{r_0}{1+ky\left(n\right)}-d-ax\left(n\right)-by\left(n\right)],}\hfill \\ \hfill y(n+1)& =& y\left(n\right)exp[m-hy(n)+cbx(n\left)\right].\hfill \end{array}$$

(10)

The authors paid attention to the local and global stability property of the equilibria of the system. They showed that due to the fear of predation, the prey species will be driven to extinction while the predator species tends to be stable since it has other food resources.

Given this background, it is natural to propose a discrete Leslie–Gower discrete predator–prey model with fear effect in prey species, and to study the influence of fear effect. In this paper, based on the model (4), we propose the following discrete Leslie–Gower predator–prey model:

$$\left\{\begin{array}{cc}\hfill & x(n+1)=x\left(n\right)exp{\displaystyle \left(\frac{r_0}{1+ky\left(n\right)}-d-ax\left(n\right)-by\left(n\right)\right)},\hfill \\ \hfill & y(n+1)=y\left(n\right)exp{\displaystyle \left(h-\frac{y\left(n\right)}{m+x\left(n\right)}\right)},\hfill \end{array}\right.$$

(11)

where $x\left(n\right)$ and $y\left(n\right)$ are the population density of prey and predator at the nth-generation, respectively. k represents the level of fear that drives the anti-predatory behavior of prey, $r_0$ and d represent the birth rate and natural mortality rate of the prey population, respectively, a represents the mortality of the prey population due to intraspecific competition, b represents the maximum capture rate of predators, and h represents the intrinsic growth rate of predators. m represents the environmental protection rate for predators, which reflects the fact that the predators have other food resources. $r_0$, k, d, a, b, h, and m are all positive constants. Since there seems to be no research on discrete Leslie–Gower predator–prey systems with fear effect, this study would be an interesting attempt to complement existing research on Leslie–Gower predator–prey systems.

As far as system (11) is concerned, there are several issues needing to be solved: (1) Is it possible to give a thorough analysis of the local and global stability of the equilibrium? (2) Is it possible to investigate the bifurcation phenomenon of the system? (3) To what extent does fear effect k change the dynamic behaviors of the system?

We investigate these three issues in this paper. The rest of the paper is arranged as follows: We investigate the existence, local stability, and global stability of equilibria in Section 2, Section 3 and Section 4, respectively. We then investigate the bifurcation phenomenon in Section 5. Some numerical simulations are carried out to show the feasibility of the main results in Section 6. We end this paper with a brief discussion.

2. The Existence of Equilibria

The equilibria of system (11) are determined by the following equations:

$$\left\{\begin{array}{cc}\hfill & x=xexp{\displaystyle \left(\frac{r_0}{1+ky}-d-ax-by\right)},\hfill \\ \hfill & y=yexp{\displaystyle \left(h-\frac{y}{m+x}\right)}.\hfill \end{array}\right.$$

(12)

By solving Equation (12), we can obtain the following conclusions:

Theorem 1.

System (11) always has two boundary equilibria given by $E_0(0,0),E_1(0,hm)$; if $r_0>d$ holds, there is also a boundary equilibrium ${\displaystyle E_2\left(\frac{r_0-d}{a},0\right)}$.

Theorem 2.

Assuming that

$$r_0>(d+bhm)(1+hkm)$$

(13)

holds, system (11) admits a unique positive equilibrium $E^{*}\left(x^{*},y^{*}\right)$, where $y^{*}=h(m+x^{*})$, $x^{*}$ is the unique positive solution of the equation:

$$A_2{x}^2+A_{1}x+A_0=0,$$

(14)

with

$$\begin{array}{cc}\hfill & A_2=b{h}^{2}k+ahk>0,\hfill \\ \hfill & A_1=2b{h}^{2}mk+ahmk+dhk+bh+a>0,\hfill \\ \hfill & A_0=b{h}^{2}mk+dhmk+bhm-d-r_0=(d+bhm)(1+hkm)-r_0<0.\hfill \end{array}$$

(15)

Proof. 

Since $x\ne 0,y\ne 0$, Equation (12) is equivalent to:

$$\left\{\begin{array}{cc}\hfill & {\displaystyle \frac{r_{0}x}{1+ky}-d-ax-by=0,}\hfill \\ \hfill & {\displaystyle h-\frac{y}{m+x}=0.}\hfill \end{array}\right.$$

(16)

The second equation of (16) is equivalent to $y=h(m+x)$. Substituting it into the first equation and simplifying it, we obtain Equation (14). Obviously, both $A_1$ and $A_2$ are positive, so Equation (14) has a unique positive solution $x^{*}$ if and only if $A_0<0$, i.e., $r_0>(d+bhm)(1+hkm)$. Consequently, the system (11) has a unique positive equilibrium $E^{*}(x^{*},y^{*})$, where $y^{*}=h(m+x^{*}).$

The proof of Theorem 2 is finished. □

Remark 1.

One could easily see that (14) has no positive solution, consequently, if condition (13) does not hold, the system (11) does not admit any positive equilibrium. So we can draw the following conclusions:

The equilibria of the system (11) are:

(1) 

$E_0(0,0)$, $E_1(0,hm)$, if $r_0\le d$;

(2) 

$E_0(0,0)$, $E_1(0,hm)$, and ${\displaystyle E_2\left(\frac{r_0-d}{a},0\right)}$, if $d<r_0\le (d+bhm)(1+hkm)$;

(3) 

$E_0(0,0)$, $E_1(0,hm)$, ${\displaystyle E_2\left(\frac{r_0-d}{a},0\right)}$, and $E^{*}(x^{*},y^{*})$, if $r_0>(d+bhm)(1+hkm)$.

Remark 2.

Assume $x^{*}=x^{*}\left(k\right)$ is the unique positive root of Equation (14), i.e.,

$$A_2\left(k\right){\left(x^{*}\left(k\right)\right)}^2+A_1\left(k\right)x^{*}\left(k\right)+A_0\left(k\right)=0,$$

(17)

where $A_0\left(k\right)<0$. Computing the derivative of (17) with respect to k, we obtain

$$A_2^{\prime }\left(k\right){\left(x^{*}\left(k\right)\right)}^2+2A_2\left(k\right)x^{*}\left(k\right)x^{*\prime }\left(k\right)+A_1^{\prime }\left(k\right)x^{*}\left(k\right)+A_1\left(k\right)x^{*\prime }\left(k\right)+A_0^{\prime }\left(k\right)=0,$$

where "${}^{\prime }$" represents the derivative ${\displaystyle \frac{\mathrm{d}}{\mathrm{d}k}}$. Simple calculation and analysis can reveal that

$${\displaystyle \frac{\mathrm{d}{x}^{*}\left(k\right)}{\mathrm{d}k}=-\frac{A_2^{\prime }\left(k\right){\left(x^{*}\left(k\right)\right)}^2+A_1^{\prime }\left(k\right)x^{*}\left(k\right)+A_0^{\prime }\left(k\right)}{2A_2\left(k\right)x^{*}\left(k\right)+A_1\left(k\right)}<0,}$$

which means that the value of prey equilibrium $x^{*}$ is strictly deceasing function of k.

Due to

$$y^{*}\left(k\right)=h\left(m+x^{*}\left(k\right)\right),$$

we can conclude that the value of predator equilibrium $y^{*}$ is also is the decreasing function of k,however, $y^{*}\left(k\right)\ge hm.$

3. The Local Stability of Equilibria

The Jacobian matrix of system (11) at the equilibrium $E(x,y)$ is

$$J\left(E\right)=\left(\begin{array}{cc}{\displaystyle (1-ax)exp{\displaystyle \left(\frac{r_0}{1+ky}-d-ax-by\right)}}& {\displaystyle x\left(-\frac{r_{0}k}{{(1+ky)}^2}-b\right)exp{\displaystyle \left(\frac{r_0}{1+ky}-d-ax-by\right)}}\\ {\displaystyle \frac{y^2}{{(m+x)}^2}exp{\displaystyle \left(h-\frac{y}{m+x}\right)}}& {\displaystyle \left(1-\frac{y}{m+x}\right)exp{\displaystyle \left(h-\frac{y}{m+x}\right)}}\end{array}\right).$$

Next, we judge the local stability of the equilibria of system (11) according to the above Jacobian matrix and Lemmas 1 and 2 in [35].

3.1. The Local Stability of Boundary Equilibria $E_0,E_1,E_2$

Theorem 3.

$E_0(0,0)$ is

(1) 

A source if $r_0>d$;

(2) 

A saddle if $r_0<d$;

(3) 

Non-hyperbolic if $r_0=d$.

Proof. 

The Jacobian matrix of system (11) at $E_0(0,0)$ is

$$J\left(E_0\right)=\left(\begin{array}{cc}exp(r_0-d)& 0\\ 0& exp\left(h\right)\end{array}\right).$$

Obviously, the two eigenvalues of $J\left(E_0\right)$ are ${\lambda }_1=exp(r_0-d)>0$ and ${\lambda }_2=exp\left(h\right)>1.$ Hence, if $r_0>d$, i.e., ${\lambda }_1>1$, then $E_0(0,0)$ is a source according to Lemma 1 of [35]. Thus, (1) holds. Numbers (2) and (3) can be proved in the same way. Therefore, $E_0$ is always unstable.

The proof of Theorem 3 is finished. □

Theorem 4.

$E_1(0,hm)$ is

(1) 

A sink if $r_0<(d+bhm)(1+hkm)$ and $h<2$;

(2) 

A source if $r_0>(d+bhm)(1+hkm)$ and $h>2$;

(3) 

A saddle if one of the following conditions holds:

(i) 

$r_0<(d+bhm)(1+hkm)$ and $h>2$;

(ii) 

$r_0>(d+bhm)(1+hkm)$ and $h<2$;

(4) 

Non-hyperbolic if $r_0=(d+bhm)(1+hkm)$ or $h=2$.

Proof. 

The Jacobian matrix of system (11) at $E_1(0,hm)$ is

$$J\left(E_1\right)=\left(\begin{array}{cc}exp{\displaystyle \left({\displaystyle \frac{r_0}{1+hkm}}-d-bhm\right)}& 0\\ h^2& 1-h\end{array}\right);$$

the two eigenvalues of $J\left(E_1\right)$ are ${\lambda }_1=exp{\displaystyle \left({\displaystyle \frac{r_0}{1+hkm}}-d-bhm\right)}>0$ and ${\lambda }_2=1-h<1.$ Hence, if $r_0<(d+bhm)(1+hkm)$ and $h<2$, we have $0<{\lambda }_1<1$ and $|{\lambda }_2|<1$, then $E_1(0,hm)$ is a sink according to Lemma 1 of [35]. Thus, (1) holds, and (2)–(4) can be proved in the same way. Therefore, $E_1$ is stable if and only if $r_0<(d+bhm)(1+hkm)$ and $h<2$.

The proof of Theorem 4 is finished. □

Theorem 5.

If the predator-free equilibrium ${\displaystyle E_2\left({\displaystyle \frac{r_0-d}{a}},0\right)}$ of system (11) exists (i.e., $r_0>d$ holds), we can easily obtain that $E_2$ is

(1) 

A source if $r_0-d>2$;

(2) 

A saddle if $r_0-d<2$;

(3) 

Non-hyperbolic if $r_0-d=2$.

Proof. 

The Jacobian matrix of system (11) at ${\displaystyle E_2\left({\displaystyle \frac{r_0-d}{a}},0\right)}$ is

$$J\left(E_2\right)=\left(\begin{array}{cc}1-(r_0-d)& {\displaystyle -\frac{(r_0-d)(k{r}_0+b)}{a}}\\ 0& exp\left(h\right)\end{array}\right).$$

Obviously, the two eigenvalues of $J\left(E_2\right)$ are ${\lambda }_1=1-(r_0-d)<1$ and ${\lambda }_2=exp\left(h\right)>1$. Similar to the proof in Theorem 4, the above results can be easily obtained.

The proof of Theorem 5 is finished. □

The above is the correlation analysis of the local stability of the boundary equilibria of the system (11). Next, we discuss the local stability of the positive equilibrium of system (11).

3.2. The Local Stability of the Positive Equilibrium $E^{*}$

Theorem 6.

If the positive equilibrium $E^{*}(x^{*},y^{*})$ of system (11) exists, i.e., $r_0>(d+bhm)(1+hkm)$, we can easily obtain that $E^{*}$ is

(1) 

A sink if $P<Q⩽2$ or $0<2Q-4<P<Q$;

(2) 

A source if one of the following conditions holds:

(i) 

$Q⩽2$ and $P>Q$;

(ii) 

$Q>2$ and $P>max(Q,2Q-4)$;

(3) 

A saddle if $Q>2$ and $P<2Q-4$;

(4) 

Non-hyperbolic if $P=2Q-4$, where

$$\begin{array}{cc}\hfill & P=h\left({\displaystyle a+\frac{h{r}_{0}k}{{(1+k{y}^{*})}^2}+bh}\right)x^{*},\hfill \\ \hfill & Q=a{x}^{*}+h.\hfill \end{array}$$

Hence, $E^{*}$ is locally asymptotically stable if $P<Q⩽2$ or $0<2Q-4<P<Q$; otherwise, it is unstable.

Proof. 

The Jacobian matrix of system (11) at $E^{*}(x^{*},y^{*})$ is

$$J\left(E^{*}\right)=\left(\begin{array}{cc}1-a{x}^{*}& {\displaystyle -x^{*}\left(\frac{r_{0}k}{{(1+k{y}^{*})}^2}+b\right)}\\ h^2& 1-h\end{array}\right).$$

Therefore, the characteristic equation of $J\left(E^{*}\right)$ is $F\left(\lambda \right)={\lambda }^2+B\lambda +C=0$, where

$$\begin{array}{cc}\hfill & B=-2+a{x}^{*}+h,\hfill \\ \hfill & C=(1-a{x}^{*})(1-h)+h^2{x}^{*}\left({\displaystyle \frac{r_{0}k}{{(1+k{y}^{*})}^2}+b}\right).\hfill \end{array}$$

It can be obtained through calculation that:

$$\begin{array}{cc}\hfill & F\left(1\right)=1+B+C=h\left({\displaystyle a+\frac{h{r}_{0}k}{{(1+k{y}^{*})}^2}+bh}\right)x^{*}=P>0,\hfill \\ \hfill & F(-1)=1-B+C=4-2h-2a{x}^{*}+h\left({\displaystyle a+\frac{h{r}_{0}k}{{(1+k{y}^{*})}^2}+bh}\right)x^{*}=4-2Q+P,\hfill \\ \hfill & C-1=-h-a{x}^{*}+h\left({\displaystyle a+\frac{h{r}_{0}k}{{(1+k{y}^{*})}^2}+bh}\right)x^{*}=P-Q,\hfill \end{array}$$

where

$$\begin{array}{c}\hfill {\displaystyle x^{*}={\displaystyle \frac{-A_1+\sqrt{A_1^2-4A_2{A}_0}}{2A_2}},y^{*}=h(m+x^{*}),}\end{array}$$

(18)

and $A_2,A_1,A_0$ are given in (15).

If $P<Q⩽2$ or $0<2Q-4<P<Q$ holds, then $F(-1)>0$, $C<1$. Hence, according to Lemmas 1 and 2 in [35], we obtain that $E^{*}$ is a sink, which is stable. Numbers (2)–(4) can be proved in the same way.

The proof of Theorem 6 is finished. □

Remark 3.

From Remark 1, it is easy to know that if $k=0$ and $r_0>d+bhm$ hold, the system (11) admits a unique positive equilibrium $E^{*}(x^{*},y^{*})$, where

$${\displaystyle x^{*}=-{\displaystyle \frac{A_0}{A_1}}={\displaystyle \frac{r_0-d-bhm}{a+bh}},y^{*}=h(m+x^{*}),}$$

which can not be obtained by equalities (18).

4. The Global Attractivity of Equilibria

4.1. The Global Attractivity of the Positive Equilibrium $E^{*}$

According to Theorem 6, the positive Equilibrium $E^{*}$ is locally asymptotically stable when $P<Q⩽2$ or $0<2Q-4<P<Q$. Next, we further explore global attractivity by using Lemmas 3–5 from [35].

Theorem 7.

If $(1+k{U}_1^y)(d+b{U}_1^y)<r_0\le d+1$ and $(1+a)/b<h\le 1$, the positive equilibrium $E^{*}$ is globally attractive, i.e.,

$$\underset{n\to \infty }{lim}x\left(n\right)=x^{*},\underset{n\to \infty }{lim}y\left(n\right)=y^{*},$$

(19)

where $U_1^y=h(m+U_1^x)+\epsilon$, $U_1^x=(r_0-d)/a+\epsilon$, ε is an arbitrarily small positive number.

Proof. 

Let $\left(x\left(n\right),y\left(n\right)\right)$ be any positive solution of system (11) with $x\left(0\right)>0$ and $y\left(0\right)>0$. Let

$$\begin{array}{cc}\hfill & S_1=\underset{n\to \infty }{lim\; sup}x\left(n\right),I_1=\underset{n\to \infty }{lim\; inf}x\left(n\right),\hfill \\ \hfill & S_2=\underset{n\to \infty }{lim\; sup}y\left(n\right),I_2=\underset{n\to \infty }{lim\; inf}y\left(n\right).\hfill \end{array}$$

Therefore, to prove (19), we need only prove the following equations:

$$S_1=I_1=x^{*},S_2=I_2=y^{*}.$$

(20)

If there exist four sequences $\left\{U_n^x\right\}$, $\left\{U_n^y\right\}$, $\left\{V_n^x\right\}$, $\left\{V_n^y\right\}$ satisfying:

$$V_n^x\le I_1\le S_1\le U_n^x,V_n^y\le I_2\le S_2\le U_n^y,$$

(21)

and

$$\underset{n\to \infty }{lim}{U}_n^x=\underset{n\to \infty }{lim}{V}_n^x=x^{*},\underset{n\to \infty }{lim}{U}_n^y=\underset{n\to \infty }{lim}{V}_n^y=y^{*},$$

(22)

then Equation (20) is clearly established according to the squeezing of the limit. Below, we construct these four sequences by the iterative method:

(I) 

Through the first iteration, we can prove that there exist $U_1^x,U_1^y>0$, such that $S_1\le U_1^x$ and $S_2\le U_1^y$.

(i) 

Given by the first equation of system (11), we get

$$\begin{array}{cc}\hfill x(n+1)& =x\left(n\right)exp\left({\displaystyle \frac{r_0}{1+ky\left(n\right)}-d-ax\left(n\right)-by\left(n\right)}\right)\hfill \\ \hfill & \le x\left(n\right)exp\left(r_0-d-ax\left(n\right)\right),n=0,1,2,\cdots \hfill \end{array}$$

Consider the auxiliary equation

$$u(n+1)=u\left(n\right)exp\left(r_0-d-au\left(n\right)\right);$$

(23)

since $r_0-d\le 1$, we obtain $u\left(n\right)\le 1/a$ for all $n\ge 2$, where $u\left(n\right)$ is any solution of system (23) with $u\left(0\right)>0$ according to Lemma 3 in [35]. From Lemma 4 in [35], we have $f\left(u\right)=uexp\left(r_0-d-au\left(n\right)\right)$ is nondecreasing for $u\in (0,1/a]$. Hence, according to Lemma 5 in [35], we have that $x\left(n\right)\le u\left(n\right)$ for all $n\ge 2$, where $u\left(n\right)$ is any solution of (23) satisfying the initial condition $u\left(2\right)=x\left(2\right)$. Furthermore, because $r_0-d<2$, we have $\underset{n\to \infty }{lim}{u}_n=(r_0-d)/a$ according to Lemma 3 in [35]. Therefore, we have

$$S_1=\underset{n\to \infty }{lim\; sup}x\left(n\right)\le \underset{n\to \infty }{lim}u\left(n\right)={\displaystyle \frac{r_0-d}{a}}.$$

Hence, for sufficiently small $\epsilon >0$ there exists an integer $N_1>2$ such that if $n\ge N_1$, then

$$x\left(n\right)\le {\displaystyle \frac{r_0-d}{a}}+\epsilon :=U_1^x.$$

(ii) 

Given by the second equation of system (11), we get

$$\begin{array}{cc}\hfill y(n+1)& =y\left(n\right)exp\left({\displaystyle h-\frac{y\left(n\right)}{m+x\left(n\right)}}\right)\hfill \\ \hfill & \le y\left(n\right)exp\left({\displaystyle h-\frac{y\left(n\right)}{m+U_1^x}}\right),\forall n\ge N_1.\hfill \end{array}$$

Consider the auxiliary equation

$$u(n+1)=u\left(n\right)exp\left({\displaystyle h-\frac{y\left(n\right)}{m+U_1^x}}\right).$$

Since $h\le 1$, we obtain that $u\left(n\right)\le m+U_1^x$ for all $n\ge N_1$ from Lemma 3 in [35]. According to Lemma 4 in [35], we have that $f\left(u\right)=uexp\left({\displaystyle h-\frac{y\left(n\right)}{m+U_1^x}}\right)$ is nondecreasing for $u\in (0,m+U_1^x]$. Hence, according to Lemma 5 in [35], we have that $y\left(n\right)\le u\left(n\right)$ for all $n\ge N_1$. Therefore, as with (I) (i), we get

$$S_2=\underset{n\to \infty }{lim\; sup}y\left(n\right)\le \underset{n\to \infty }{lim}u\left(n\right)=h(m+U_1^x).$$

Hence, for sufficiently small $\epsilon >0$ there exists an integer $N_2>N_1$ such that if $n\ge N_2$, then

$$y\left(n\right)\le h(m+U_1^x)+\epsilon :=U_1^y.$$

(II) 

Through the second iteration, we can prove that there exist $V_1^x,V_1^y>0$, such that $I_1\ge V_1^x$ and $I_2\ge V_1^y$.

(i) 

Given by the first equation of system (11), we obtain that

$$x(n+1)\ge x\left(n\right)exp\left({\displaystyle \frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y-ax\left(n\right)}\right),\forall n\ge N_2.$$

Consider the auxiliary equation

$$u(n+1)=u\left(n\right)exp\left({\displaystyle \frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y-au\left(n\right)}\right).$$

(24)

Since $(1+k{U}_1^y)(d+b{U}_1^y)<r_0\le d+1$, we obtain

$${\displaystyle 0<\frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y\le r_0-d\le 1.}$$

Hence $u\left(n\right)\le 1/a$ for all $n\ge N_2$ according to Lemma 3 in [35]. From Lemma 4 in [35], we know that $f\left(u\right)=uexp\left({\displaystyle \frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y-au}\right)$ is nondecreasing for $u\in (0,1/a]$. Hence, according to Lemma 5 in [35], we have that $x\left(n\right)\ge u\left(n\right)$ for all $n\ge N_2$. According to Lemma 3 in [35], we have

$${\displaystyle I_1=\underset{n\to \infty }{lim\; inf}x\left(n\right)\ge \underset{n\to \infty }{lim}u\left(n\right)={\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y}\right).}$$

Hence, for sufficiently small $\epsilon >0$, there exists an integer $N_3>N_2$ such that if $n\ge N_3$, then

$${\displaystyle x\left(n\right)\ge {\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y}\right)-\epsilon :=V_1^x.}$$

(ii) 

Given by the second equation of system (11), we obtain

$$y(n+1)\ge y\left(n\right)exp\left({\displaystyle h-\frac{y\left(n\right)}{m+V_1^x}}\right),\forall n\ge N_3.$$

Consider the auxiliary equation

$$u(n+1)=u\left(n\right)exp\left(h-\frac{u\left(n\right)}{m+V_1^x}\right).$$

Since $h\le 1$, with a similar argument as above, we can obtain

$$I_2=\underset{n\to \infty }{lim\; inf}y\left(n\right)\ge \underset{n\to \infty }{lim}u\left(n\right)=h(m+V_1^x).$$

Hence, for sufficiently small $\epsilon >0$ there exists an integer $N_4>N_3$ such that if $n\ge N_4$, then

$$y\left(n\right)\ge h(m+V_1^x)-\epsilon :=V_1^y.$$

(III) 

Through the third iteration, we can prove that there exist $U_2^x\le U_1^x,U_2^y\le U_1^y$, such that $S_1\le U_2^x$ and $S_2\le U_2^y$.

(i) 

Given by the first equation of system (11), we obtain that

$$x(n+1)\le x\left(n\right)exp\left({\displaystyle \frac{r_0}{1+k{V}_1^y}-d-b{V}_1^y-ax\left(n\right)}\right),\forall n\ge N_4.$$

Since $U_1^y>V_1^y$, we obtain ${\displaystyle 0<\frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y<\frac{r_0}{1+k{V}_1^y}-d-b{V}_1^y\le r_0-d\le 1}$, according to Lemma 5 in [35], we have that

$$S_1=\underset{n\to \infty }{lim\; sup}x\left(n\right)\le {\displaystyle \frac{1}{a}}\left(\frac{r_0}{1+k{V}_1^y}-d-b{V}_1^y\right).$$

Hence, for sufficiently small $\epsilon >0$, there exists an integer $N_5>N_4$ such that if $n\ge N_5$, then

$$x\left(n\right)\le {\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{V}_1^y}-d-b{V}_1^y}\right)+{\displaystyle \frac{\epsilon }{2}}:=U_2^x\le U_1^x.$$

(ii) 

Given by the second equation of system (11), we obtain

$$y(n+1)\le y\left(n\right)exp\left(h-\frac{y\left(n\right)}{m+U_2^x}\right),\forall n\ge N_5.$$

Since $h\le 1$, with a similar argument as above, we can obtain

$$S_2=\underset{n\to \infty }{lim\; sup}y\left(n\right)\le h(m+U_2^x).$$

Hence, for sufficiently small $\epsilon >0$, there exists an integer $N_6>N_5$ such that if $n\ge N_6$, then

$$y\left(n\right)\le h(m+U_2^x)+\frac{\epsilon }{2}:=U_2^y\le U_1^y.$$

(IV) 

Through the fourth iteration, we can prove that there exist $V_2^x\ge V_1^x,V_2^y\ge V_1^y$, such that $I_1\ge V_2^x$ and $I_2\ge V_2^y$.

(i) 

Similarly, from the first equation of system (11), we obtain

$$x(n+1)\ge x\left(n\right)exp\left({\displaystyle \frac{r_0}{1+k{U}_2^y}-d-b{U}_2^y-ax\left(n\right)}\right),\forall n\ge N_6.$$

Since $U_1^y\ge U_2^y$, we can obtain

$${\displaystyle 0<\frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y\le \frac{r_0}{1+k{U}_2^y}-d-b{U}_2^y\le r_0-d\le 1.}$$

According to Lemma 5 in [35], we have

$$I_1=\underset{n\to \infty }{lim\; inf}x\left(n\right)\ge {\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_2^y}-d-b{U}_2^y}\right).$$

Hence, for sufficiently small $\epsilon >0$ there exists an integer $N_7>N_6$ such that if $n\ge N_7$, then

$$x\left(n\right)\ge {\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_2^y}-d-b{U}_2^y}\right)-{\displaystyle \frac{\epsilon }{2}}:=V_2^x\ge V_1^x.$$

(ii) 

Similarly, from the second equation of system (11), we obtain

$$y(n+1)\ge y\left(n\right)exp\left(h-\frac{y\left(n\right)}{m+V_2^x}\right),\forall n\ge N_7.$$

Since $h\le 1$, with a similar argument as above, we can obtain

$$I_2=\underset{n\to \infty }{lim\; inf}y\left(n\right)\ge h(m+V_2^x).$$

Hence, for sufficiently small $\epsilon >0$ there exists an integer $N_8>N_7$ such that if $n\ge N_8$, then

$$y\left(n\right)\ge h(m+V_2^x)-\frac{\epsilon }{2}:=V_2^y\ge V_1^y.$$

By repeating steps (III)–(IV), sequences $\left\{U_n^x\right\}$, $\left\{U_n^y\right\}$, $\left\{V_n^x\right\}$, $\left\{V_n^y\right\}$, are constructed as follows (define $V_0^y$ as 0):

$$\left\{\begin{array}{cc}\hfill & U_n^x={\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{V}_{n-1}^y}-d-b{V}_{n-1}^y}\right)+{\displaystyle \frac{\epsilon }{n}},\hfill \\ \hfill & U_n^y=h(m+U_n^x)+\frac{\epsilon }{n},\hfill \\ \hfill & V_n^x={\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_n^y}-d-b{U}_n^y}\right)-{\displaystyle \frac{\epsilon }{n}},\hfill \\ \hfill & {\displaystyle V_n^y=h(m+V_n^x)-\frac{\epsilon }{n}.}\hfill \end{array}\right.$$

(25)

These four sequences clearly satisfy inequalities (21), and it is shown below that they also satisfy equalities (22).

We first prove that the limits of all four sequences exist, which can be proved by using the monotone bounded criterion of sequence. First, it is proved by mathematical induction that all four sequences are monotone, and that $\left\{U_n^x\right\}$ and $\left\{U_n^y\right\}$ are monotonically decreasing, while $\left\{V_n^x\right\}$ and $\left\{V_n^y\right\}$ are monotonically increasing.

Obviously, if $n=1$, we have

$$U_1^x\ge U_2^x,V_1^x\le V_2^x,U_1^y\ge U_2^y,V_1^y\le V_2^y.$$

For $n=i-1$ ($i\ge 2$), we assume that $U_{i-1}^x\ge U_i^x$ and $V_{i-1}^x\le V_i^x$ hold. Then, further, we have

$$\begin{array}{cc}\hfill & U_{i-1}^y=h(m+U_{i-1}^x)+\frac{\epsilon }{i-1}\ge h(m+U_i^x)+\frac{\epsilon }{i}=U_i^y,\hfill \\ \hfill & V_{i-1}^y=h(m+V_{i-1}^x)-\frac{\epsilon }{i-1}\le h(m+V_i^x)-\frac{\epsilon }{i}=V_i^y.\hfill \end{array}$$

Consequently,

$$\begin{array}{cc}\hfill U_i^x={\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{V}_{i-1}^y}-d-b{V}_{i-1}^y}\right)+{\displaystyle \frac{\epsilon }{i}}& \ge {\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{V}_i^y}-d-b{V}_i^y}\right)+{\displaystyle \frac{\epsilon }{i+1}}=U_{i+1}^x,\hfill \\ \hfill U_i^y=h(m+U_i^x)+\frac{\epsilon }{i}& \ge h(m+U_{i+1}^x)+\frac{\epsilon }{i+1}=U_{i+1}^y,\hfill \\ \hfill V_i^x={\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_i^y}-d-b{U}_i^y}\right)-{\displaystyle \frac{\epsilon }{i}}& \le {\displaystyle \frac{1}{a}}\left({\displaystyle \frac{r_0}{1+k{U}_{i+1}^y}-d-b{U}_{i+1}^y}\right)-{\displaystyle \frac{\epsilon }{i+1}}=V_{i+1}^x,\hfill \\ \hfill V_i^y=h(m+V_i^x)-\frac{\epsilon }{i}& \le h(m+V_{i+1}^x)-\frac{\epsilon }{i+1}=V_{i+1}^y.\hfill \end{array}$$

Therefore, these four sequences are all monotone by mathematical induction. Combining with inequalities (21), it is easy to obtain that sequences $\left\{U_n^x\right\}$ and $\left\{U_n^y\right\}$ have lower bounds and monotonic decreasing, and sequences $\left\{V_n^x\right\}$ and $\left\{V_n^y\right\}$ have upper bounds and monotonic increasing. Therefore, it can be known from the monotone bounded criterion that the limits of these four sequences all exist. Let

$$\underset{n\to \infty }{lim}{U}_n^x=x_1,\underset{n\to \infty }{lim}{V}_n^x=x_2,\underset{n\to \infty }{lim}{U}_n^y=y_1,\underset{n\to \infty }{lim}{V}_n^y=y_2.$$

(26)

Obviously, they are all positive. Then, it is only necessary to prove that $x_1=x_2=x^{*}$, $y_1=y_2=y^{*}$.

Let us take the limit of both sides of the equality sign of the four equations in (25); taking $n\to \infty$, we obtain

$$\left\{\begin{array}{cc}\hfill & {\displaystyle x_1={\displaystyle \frac{1}{a}}\left(\frac{r_0}{1+k{y}_2}-d-b{y}_2\right),}\hfill \\ \hfill & y_1=h(m+x_1),\hfill \\ \hfill & {\displaystyle x_2={\displaystyle \frac{1}{a}}\left(\frac{r_0}{1+k{y}_1}-d-b{y}_1\right),}\hfill \\ \hfill & y_2=h(m+x_2).\hfill \end{array}\right.$$

(27)

Note that the system of Equation (27) is equivalent to (16) if $x_1=x_2,y_1=y_2$, and $U_1^y=h(m+U_1^x)+\epsilon$, so it can be known from $r_0>(1+k{U}_1^y)(d+b{U}_1^y)$ and from Theorem 2 that the positive equilibrium exists at this time, and therefore the system of Equation (27) has a positive solution $x_1=x_2=x^{*},y_1=y_2=y^{*}$ under the assumption of Theorem 7. Following we will prove that the positive solution of the system (27) is unique.

The following can be obtained by simplification of (27):

$$\left\{\begin{array}{cc}\hfill & \left[a{x}_1+d+bh(m+x_2)\right]\left[1+kh(m+x_2)\right]=r_0,\hfill \\ \hfill & \left[a{x}_2+d+bh(m+x_1)\right]\left[1+kh(m+x_1)\right]=r_0.\hfill \end{array}\right.$$

(28)

Similarly, we know that (28) has at least one positive solution $x_1=x_2=x^{*}$. Below, we prove that Equation (28) has only one positive solution:

Otherwise, Equation (28) has other positive solutions, that is, there are positive numbers $x_1\ne x_2$ satisfying (28). Subtracting the first equation of (28) from the second equation, we obtain

$$(x_1-x_2)\left[-2b{h}^{2}km-b{h}^{2}k(x_1+x_2)+ahkm-dhk+hkm-bh+a+1\right]=0.$$

(29)

Because $x_1\ne x_2$, it can be obtained from (29) that

$$x_1+x_2={\displaystyle \frac{-(d+bhm)hk-(1+hkm)(bh-a-1)}{b{h}^{2}k}}.$$

Since $bh-a-1>0$, we have $x_1+x_2<0$, which contradicts the fact that both $x_1$ and $x_2$ are positive.

Hence, Equation (28) has only one positive solution $x_1=x_2=x^{*}$. Then, we have $y_1=y_2=y^{*}$ from (27). Thus, we prove that these four sequences also satisfy the equalities (22). Therefore, the positive equilibrium $E^{*}$ is globally attractive.

The proof of Theorem 7 is finished. □

4.2. The Global Attractivity of the Boundary Equilibrium $E_1$

In the previous section, we proved that the positive equilibrium is globally attractive under certain conditions, which indicates that predator and prey populations can coexist stably under certain conditions. When the prey population goes extinct, the predator population probably does not go extinct with it, and tends to stabilize because there are other food sources. This section discusses this possibility, and proves that the boundary equilibrium $E_1$ is also globally attractive under certain conditions.

Theorem 8.

If $r_0>d$, any positive solution of system (11) is ultimately bounded.

Proof. 

Let $\left(x\right(n),y(n\left)\right)$ be any positive solution of system (11). Since

$$x(n+1)\le x\left(n\right)exp[r_0-d-ax\left(n\right)]$$

and $r_0-d>0$, then from Lemma 6 in [35] we obtain

$$\underset{n\to \infty }{lim\; sup}x\left(n\right)\le {\displaystyle \frac{exp(r_0-d-1)}{a}}.$$

Therefore, for any sufficiently small $\epsilon >0$ there exists an integer $N_1>0$ such that if $n\ge N_1$, then

$$x\left(n\right)\le {\displaystyle \frac{exp(r_0-d-1)}{a}}+\epsilon :=U.$$

(30)

Then, from the second equation of the system (11) and (30), we obtain

$$y(n+1)\le y\left(n\right)exp{\displaystyle \left(h-\frac{y\left(n\right)}{m+U}\right)},n\ge N_1.$$

According to Lemma 6 in [35], we have

$$\underset{n\to \infty }{lim\; sup}y\left(n\right)\le (m+U)exp(h-1).$$

Therefore, for any sufficiently small $\epsilon >0$ there exists an integer $N_2>N_1$ such that if $n\ge N_2$, then

$$y\left(n\right)\le (m+U)exp(h-1)+\epsilon :=V.$$

(31)

Combining (30) and (31), it can be concluded that any positive solution of system (11) is ultimately bounded.

The proof of Theorem 8 is finished. □

Next we use Lemmas 6 and 7 from [35] and Theorem 8 to discuss the global attractivity of $E_1$.

Theorem 9.

Assuming that

$$r_0<d+bhm$$

(32)

and

$$h<ln2+1$$

(33)

hold, the boundary equilibrium $E_1(0,hm)$ is globally attractive, i.e.,

$$\underset{n\to \infty }{lim}x\left(n\right)=0,\underset{n\to \infty }{lim}y\left(n\right)=hm,$$

where $\left(x\right(n),y(n\left)\right)$ is any positive solution of system (11).

Proof. 

Assume that $\left(x\right(n),y(n\left)\right)$ is any positive solution of system (11). First let us prove that $\underset{n\to \infty }{lim}x\left(n\right)=0$.

For all $j\in \mathbb{N}$, according to system (11), we have

$$\begin{array}{cc}\hfill ln{\displaystyle \frac{x(j+1)}{x\left(j\right)}}& {\displaystyle =\frac{r_0}{1+ky\left(j\right)}-d-ax\left(j\right)-by\left(j\right)}\hfill \\ \hfill & \le r_0-d-ax\left(j\right)-by\left(j\right),\hfill \end{array}$$

(34)

$$\begin{array}{c}\hfill ln{\displaystyle \frac{y(j+1)}{y\left(j\right)}}=h-{\displaystyle \frac{y\left(j\right)}{m+x\left(j\right)}}.\end{array}$$

From (32), we have

$${\displaystyle \frac{r_0-d}{h}}<bm.$$

Hence, there exist the positive constants p and q such that

$${\displaystyle \frac{r_0-d}{h}}<{\displaystyle \frac{q}{p}}<bm.$$

Then,

$$bmp-q>0,$$

(35)

and

$$p(r_0-d)-qh<0.$$

Then, there exists $\delta >0$ such that

$$p(r_0-d)-qh<-\delta <0.$$

(36)

From (34)–(36), we obtain

$$\begin{array}{cc}\hfill pln{\displaystyle \frac{x(j+1)}{x\left(j\right)}}-qln{\displaystyle \frac{y(j+1)}{y\left(j\right)}}& \le \left[p(r_0-d)-hq\right]-\left[apx\left(j\right)+{\displaystyle \frac{bmp-q}{m}}y\left(j\right)\right]\hfill \\ \hfill & \le p(r_0-d)-hq<-\delta <0.\hfill \end{array}$$

Adding both sides from 0 to $n-1$, we have

$$\sum _{j=1}^n(pln{\displaystyle \frac{x(j+1)}{x\left(j\right)}}-qln{\displaystyle \frac{y(j+1)}{y\left(j\right)}})=pln{\displaystyle \frac{x\left(n\right)}{x\left(0\right)}}-qln{\displaystyle \frac{y\left(n\right)}{y\left(0\right)}}<-\delta n,$$

then

$$x\left(n\right)<x\left(0\right){\left({\displaystyle \frac{y\left(n\right)}{y\left(0\right)}}\right)}^{\frac{q}{p}}exp\left(-{\displaystyle \frac{\delta }{p}}n\right).$$

(37)

From (31) and (37), we obtain

$$x\left(n\right)<x\left(0\right){\left({\displaystyle \frac{V}{y\left(0\right)}}\right)}^{\frac{q}{p}}exp\left(-{\displaystyle \frac{\delta }{p}}n\right),\forall n\ge N_2.$$

Hence, $x\left(n\right)\to 0$ with $n\to \infty$, i.e.,

$$\underset{n\to \infty }{lim}x\left(n\right)=0.$$

(38)

Therefore, for any sufficiently small $\epsilon >0$, there exists an integer $N_3>N_2$ such that if $n\ge N_3$, then

$$x\left(n\right)<\epsilon .$$

(39)

Then, we prove that $\underset{n\to \infty }{lim}y\left(n\right)=hm.$

Consider the system:

$$y_1(n+1)=y_1\left(n\right)exp\left(h-{\displaystyle \frac{y_1\left(n\right)}{m}}\right).$$

(40)

The system has a unique positive equilibrium $y_1^{*}=hm$. Let $\left\{y_1\left(n\right)\right\}$ be any positive solution of system (40); $h<2$ from (33), so $\underset{n\to \infty }{lim}{y}_1\left(n\right)=hm$ can be obtained by applying Lemma 3 in [35]. Therefore, to prove that $\underset{n\to \infty }{lim}y\left(n\right)=hm$, we need only prove that

$$\underset{n\to \infty }{lim}\left(y\left(n\right)-y_1\left(n\right)\right)=0.$$

Assume that

$$y\left(n\right)=y_1\left(n\right)exp\left[c\left(n\right)\right];$$

(41)

then, we only need to prove $\underset{n\to \infty }{lim}c\left(n\right)=0$.

Substituting (39)–(41) into the second equation of system (11), we can obtain the following equation by simple calculation:

$$\begin{array}{cc}\hfill c(n+1)& =c\left(n\right)+{\displaystyle \frac{y_1\left(n\right)}{m}}-{\displaystyle \frac{y_1\left(n\right)exp\left[c\left(n\right)\right]}{m+x\left(n\right)}}\hfill \\ \hfill & \le \left[1-{\displaystyle \frac{y_1\left(n\right)exp\left(\theta \right(n\left)c\right(n\left)\right)}{m+\epsilon }}\right]c\left(n\right)+{\displaystyle \frac{y_1\left(n\right)\epsilon }{m^2}},\hfill \end{array}$$

(42)

where $\theta \left(n\right)\in [0,1]$, which means $y_1\left(n\right)exp\left(\theta \right(n\left)c\right(n\left)\right)\in \left[y_1\left(n\right),y\left(n\right)\right]$.

Since $y\left(n\right)exp\left[h-{\displaystyle \frac{y\left(n\right)}{m}}\right]\le y(n+1)\le y\left(n\right)exp\left[h-{\displaystyle \frac{y\left(n\right)}{m+\epsilon }}\right]$, from Lemmas 6 and 7 in [35], we obtain

$$\begin{array}{cc}\hfill & \underset{n\to \infty }{lim\; sup}y\left(n\right)\le (m+\epsilon )exp(h-1):=S,\hfill \\ & \underset{n\to \infty }{lim\; inf}y\left(n\right)\ge hmexp\left(h-{\displaystyle \frac{S}{m}}\right):=I.\hfill \end{array}$$

According to system (40), we obtain

$$\begin{array}{cc}\hfill & \underset{n\to \infty }{lim\; sup}{y}_1\left(n\right)\le mexp(h-1)\le S,\hfill \\ & \underset{n\to \infty }{lim\; inf}{y}_1\left(n\right)\ge hmexp\left(h-{\displaystyle \frac{S}{m}}\right)=I.\hfill \end{array}$$

Therefore, for any sufficiently small $\epsilon >0$, there exists an integer $N_4>N_3$ such that

$$I-\epsilon \le y\left(n\right),y_1\left(n\right)\le S+\epsilon ,\forall n\ge N_4.$$

(43)

Given (33), we obtain $exp(h-1)<2$. In addition, if $\epsilon$ is sufficiently small, we obtain

$$exp(h-1)<{\displaystyle \frac{2m}{m+\epsilon }},$$

from which we have

$$max\left\{\left|1-{\displaystyle \frac{1}{m}}S\right|,\left|1-{\displaystyle \frac{1}{m}}I\right|\right\}<1.$$

Let

$${\alpha }_{\epsilon }=max\left\{\left|1-{\displaystyle \frac{1}{m+\epsilon }}(S+\epsilon )\right|,\right|1-{\displaystyle \frac{1}{m+\epsilon }}(I-\epsilon )\left|\right\},$$

(44)

since $\epsilon$ is sufficiently small, we obtain

$${\alpha }_{\epsilon }<1.$$

(45)

From (42)–(44), we obtain

$$\begin{array}{cc}\hfill \left|c\right(n+1\left)\right|& \le max\left\{\left|1-{\displaystyle \frac{1}{m+\epsilon }}(I-\epsilon )\right|,\left|1-{\displaystyle \frac{1}{m+\epsilon }}(S+\epsilon )\right|\right\}\left|c\left(n\right)\right|+{\displaystyle \frac{exp(h-1)}{m}}\epsilon \hfill \\ \hfill & ={\alpha }_{\epsilon }\left|c\left(n\right)\right|+{\displaystyle \frac{exp(h-1)}{m}}\epsilon ,\forall n\ge N_4.\hfill \end{array}$$

(46)

then we can obtain

$${\displaystyle \left|c\left(n\right)\right|\le {\alpha }_{\epsilon }^{n-N_4}\left|c\left(N_4\right)\right|+\frac{exp(h-1)}{m}\frac{1-{\alpha }_{\epsilon }^{n-N_4}}{1-{\alpha }_{\epsilon }}\epsilon .}$$

Given (45) and because $\epsilon$ is sufficiently small, we can obtain $\underset{n\to \infty }{lim}c\left(n\right)=0$ from the above inequality.

The proof of Theorem 9 is finished. □

5. Bifurcation Analysis

In this section, we discuss bifurcations at the equilibria of the system (11). We can easily derive the following theorems by using the central manifold and bifurcation theory [37,38].

5.1. Flip Bifurcation

First, we discuss the flip bifurcation of the system (11) at $E_1(0,hm)$.

Theorem 10.

System (11) undergoes flip bifurcation at $E_1(0,hm)$ if parameters vary in the small neighborhood of $F_A=\left\{(r_0,k,d,a,b,h,m)\in {\mathbb{R}}_{+}^7:h=2,r_0\ne (d+bhm)(1+hkm)\right\}$.

Theorem 11.

System (11) undergoes flip bifurcation at $E_2\left({\displaystyle \frac{r_0-d}{a}},0\right)$ if parameters vary in the small neighborhood of $F_B=\{(r_0,k,d,a,b,h,m)\in {\mathbb{R}}_{+}^7:r_0=d+2\}$.

Proof. 

The proof of Theorem 11 is similar to the proof of Theorem 10. Therefore, we only give the proof of Theorem 10. First, we can obtain that if $h=2$ and $r_0\ne (d+bhm)(1+hkm)$ hold, the eigenvalues of $J\left(E_1\right)$ are ${\lambda }_2=-1$ and ${\lambda }_1\ne 1,-1$. Since a center manifold of system (11) at $E_1$ is $x=0$ and the system (11) restricted to it is the logistic model:

$$y\to f\left(y\right)=yexp\left(h-{\displaystyle \frac{y}{m}}\right).$$

Its nontrivial fixed point is $y_1=hm$. Then, we can easily obtain that

$$f^{\prime }\left(y_1\right)=1-h=-1$$

when parameters vary in the small neighborhood of $F_A$. Hence, flip bifurcation can occur at $E_1(0,hm)$, as shown in Figure 1.

The proof of Theorem 10 is finished. □

5.2. Transcritical Bifurcation

In this section, our claim is that fixed point $E_0(0,0)$ undergoes transcritical bifurcation.

Theorem 12.

System (11) undergoes transcritical bifurcation at its trivial equilibrium $E_0(0,0)$ if parameters vary in the small neighborhood of $F_C=\left\{(r_0,k,d,a,b,h,m)\in {\mathbb{R}}_{+}^7:r_0=d\right\}$.

Proof. 

First, we can obtain that if $r_0=d$ and $h>0$ hold, the eigenvalues of $J\left(E_0\right)$ are ${\lambda }_1=1$ and ${\lambda }_2=exp\left(h\right)\ne 1,-1$.

As $(r_0,k,d,a,b,h,m)\in F_C$, we can describe system (11) by the following map:

$$\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{c}xexp{\displaystyle \left(\frac{(r_0^{*}+\tilde{r})+x}{1+ky}-d-ax-by\right)}\\ yexp{\displaystyle \left(h-\frac{y}{m+x}\right)}\end{array}\right),$$

(47)

where $r_0^{*}=d$ and parameter $\tilde{r}$ represents a very small perturbation parameter such that $|\tilde{r}|\ll 1$. If we expand the system (47) by Taylor at $(x,y,\tilde{r})=(0,0,0)$, we have:

$$\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{cc}1& 0\\ 0& {\lambda }_2\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}f\left(x,y,\tilde{r}\right)\\ g\left(x,y,\tilde{r}\right)\end{array}\right),$$

(48)

where

$$\begin{array}{ccc}\hfill f\left(x,y,\tilde{r}\right)& =& -a{x}^2+\tilde{r}x+(-dk+b)xy+{\displaystyle \frac{1}{2}}{a}^2{x}^3-a\tilde{r}{x}^2+(adk+ab)x^{2}y+(\frac{1}{2}{d}^2{k}^2+bdk\hfill \\ & & +d{k}^2+\frac{1}{2}{b}^2)x^{2}y+\frac{1}{2}{\tilde{r}}^{2}x-(dk+b+k)\tilde{r}xy+O\left(\left(\right|x|+\left|y\right|+\left|\tilde{r}\right|{)}^4\right),\hfill \\ \hfill g\left(x,y,\tilde{r}\right)& =& {\lambda }_{2}y-\frac{{\lambda }_2}{m}{y}^2+\frac{{\lambda }_2}{m^2}x{y}^2+\frac{{\lambda }_2}{2m^2}{y}^3+O\left(\left(\right|x|+\left|y\right|+\left|\tilde{r}\right|{)}^4\right).\hfill \end{array}$$

Therefore, there exists a center manifold $W^c(0,0,0)$ for map (48) at $(0,0)$ in a small neighborhood of $\tilde{r}=0$:

$$W^c(0,0,0)=\left\{\left(x,y,\tilde{r}\right)\in {\mathbb{R}}^3|y=h\left(x,\tilde{r}\right),h(0,0)=0,Dh(0,0)=0\right\}.$$

Then, h must satisfy

$$h\left(x+f\left(x,h\left(x,\tilde{r}\right),\tilde{r}\right),\tilde{r}\right)={\lambda }_{2}h\left(x,\tilde{r}\right)+g\left(x,h\left(x,\tilde{r}\right),\tilde{r}\right).$$

(49)

h can be written as

$$h\left(x,\tilde{r}\right)=h_1{x}^2+h_{2}x\tilde{r}+h_3{\tilde{r}}^2+O\left(\left(\right|x|+|\tilde{r}{\left|\right)}^3\right).$$

(50)

Substituting (50) into (49), we can obtain $h_1=h_2=h_3=0$. Additionally, the map restricted to the center manifold $W^c(0,0,0)$ is given by:

$$F:x\to x-a{x}^2+x\tilde{r}+\frac{1}{2}{a}^2{x}^3-a{x}^2\tilde{r}+\frac{1}{2}x{\tilde{r}}^2+O\left(\left(\right|x|+|\tilde{r}{\left|\right)}^4\right).$$

Next, we establish ${\mathcal{L}}_1\ne 0$ and ${\mathcal{L}}_2\ne 0$ as follows:

$$\left\{\begin{array}{cc}\hfill & {\displaystyle {\mathcal{L}}_1={\left(\frac{{\partial }^{2}F}{\partial x^2}\right)}_{(0,0)}=-2a\ne 0,}\hfill \\ \hfill & {\displaystyle {\mathcal{L}}_2={\left(\frac{{\partial }^{2}F}{\partial x\partial \tilde{r}}\right)}_{(0,0)}=1\ne 0.}\hfill \end{array}\right.$$

Thus, we can state that system (11) undergoes transcritical bifurcation at $E_0$ if $r_0=d$.

The proof of Theorem 12 is finished. □

6. Numerical Simulations

In this section, we show the feasibility of the main results of this paper. Through numerical simulation, we directly analyze the influences of fear effect on discrete system (11).

Example 1.

The stability of the positive equilibrium $E^{*}$.

(1) 

When

$$(k,d,a,b,h,m,r_0)=(0.7,0.3,3.0,0.5,0.5,1.0,4.0),$$

we take the initial values of system (11) as $(0.6,0.3)$, we have $(1+hkm)(d+bhm)=0.74<r_0$, then the system (11) admits a unique positive equilibrium $E^{*}(0.62,0.81)$ according to Theorem 2. At this point, we have $P=1.18,Q=2.35$, then $E^{*}$ is locally asymptotically stable according to Theorem 6, as shown in Figure 2. In fact, if we take the initial values of system (11) as $(0.6,0.3),(0.9,0.2),(0.7,0.4),(0.4,0.5)(0.2,0.2),(0.8,0.8)$, we can see that $E^{*}$ is also globally stable in Figure 3, but at this point we have $r_0=4>1.3=d+1$. It shows that it is possible for $E^{*}$ to be globally attractive even if the conditions of Theorem 7 are not satisfied. This means that sufficient conditions to ensure the globally asymptotically stable equilibrium $E^{*}$ is too strict, which is probably due to the contraction and expansion of the proof.

(2) 

When

$$(k,d,a,b,h,m,r_0)=(0.0,0.3,3.0,0.5,0.5,1.0,4.0)$$

and $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$, we have $(1+hkm)(d+bhm)=0.55<r_0$, the system (11) has a unique positive equilibrium $E^{*}(1.06,1.03)$ according to Remark 3. At this point, we have $P=1.73,Q=3.68$, then $E^{*}$ is a saddle according to Theorem 6, as shown in Figure 4.

Example 2.

The stability of prey-free equilibrium $E_1$.

(1) 

When

$$(k,d,a,b,h,m,r_0)=(0.7,0.3,3.0,0.5,0.5,1.0,0.4)$$

and $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$, we have $(1+hkm)(d+bhm)=0.74>r_0$ and $h<2$, then $E_1$ is locally asymptotically stable according to Theorem 4, as shown in Figure 5. Further, we can see that $d+bhm=0.55>r_0>d$ and $h<ln2+1$, which means that $E_1$ is globally attractive according to Theorem 9. We take the initial values of the system (11) as

$$(0.6,0.3),(0.9,0.2),(0.7,0.4),(0.4,0.5),(0.2,0.2),(0.8,0.8),$$

respectively, for numerical simulation, and the results verify the accuracy of the conclusion of Theorem 9, as shown in Figure 6.

(2) 

When

$$(k,d,a,b,h,m,r_0)=(0.7,0.3,3.0,0.5,1.9,1.0,2.8)$$

we have

$$(1+hkm)(d+bhm)=2.91>r_0=2.8>1.25=d+bhm$$

and $ln2+1<h=1.9<2$, which satisfy the conditions of Theorem 4, but do not meet the requirements of Theorem 9. However, we take the initial values of the system as

$$(0.5,1.5),(0.3,2.7),(1.1,1.2),(0.7,2.4)(0.2,0.2),(0.8,1.8),$$

we can see that $E_1(0,1.9)$ is also globally attractive at this time from Figure 7.

Example 3.

The impact of fear effect on system stability is studied by numerical simulations.

(1) 

First, we consider the case where there is no fear effect. Set each parameter value of system (11) as $(k,d,a,b,h,m)=(0.0,0.3,3.0,0.5,0.5,1.0)$ and set $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$. Through numerical simulation, we find that when $r_0<2.67$, system (11) has a stable equilibrium point. When $r_0$ continues to increase, this equilibrium point first bifurcates into a two-period orbit ($2.67<r_0<3.25$), then bifurcates into a four-period orbit ($3.25<r_0<3.33$), and finally chaos occurs ($r_0>3.33$), as shown in Figure 8. We can see that with the increasing of grow rate of prey species, the dynamic behaviors of both predator and prey becomes complicated.

(2) 

Then, we consider the fear effect. Set each parameter value of system (11) as:

$$(d,a,b,m,h,r_0)=(0.3,3.0,0.5,1.0,0.5,4.0),$$

set $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$, and plot with k as the abscissa. We can see from Figure 9 that the system (11) changes from chaotic state ($0<k<0.204$) to eight-period orbit ($0.204<k<0.214$), four-period orbit ($0.214<k<0.277$), then to two-period orbit ($0.277<k<0.440$), to stable state ($0.440<k<12.545$), and finally to be the state that the prey is driven to extinction while the predator survive in a stable state ($k>12.545$). This means that the stability of system (11) increases as the fear effect of prey increases within a certain range, which is similar to the result in [35]. In addition, it can be seen from Figure 9 (c)-(d) that the positive equilibrium solution $(x^{*}\left(k\right),y^{*}\left(k\right))$ of the system (11) will decrease with the increasing of k, which is consistent with the conclusion of Remark 2. If k is sufficiently large, the prey will become extinct, while the predator will not become extinct due to the presence of other food sources.

The phase portraits correlated with Figure 9 are displayed in Figure 10, which includes orbits of periods 2, 4, and 8. When $k=0.005$, we can see chaotic sets in Figure 10a.

Example 4.

The influence of fear effect on the positive equilibrium solution of system.

Set each parameter value of system (11) as:

$$(d,a,b,m,h,r_0)=(0.3,3.0,0.5,1.0,0.5,4.0)$$

and $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$. The fear effect values were set as

$$k_1=0.6,k_2=1,k_3=2,k_4=4,k_5=8,k_6=20$$

for numerical simulation, and the number of iterations n was taken as the abscissa for plotting, as shown in Figure 11. It can be seen from Figure 11 that the larger the fear effect is, the smaller the positive equilibrium solution $(x^{*}\left(k\right),y^{*}\left(k\right))$ is. In other words, the density of prey and predator population in the stable state will decrease with the increasing of the fear effect. When the fear effect is sufficiently large, the prey population will become extinct, while the predator population will become stable. In addition, it is noted that when the fear effect is in a certain range, the larger the fear effect is, the fewer iterations it takes for the predator and prey populations to stabilize, that is, the faster the two populations stabilize. Therefore, it is concluded that the fear effect k can enhance the stability of the positive equilibrium solution $(x^{*}\left(k\right),y^{*}\left(k\right))$ and reduce the value of the positive equilibrium solution $(x^{*}\left(k\right),y^{*}\left(k\right))$, which is consistent with the conclusion of Remark 2 and Figure 9c,d.

Example 5.

The impact of fear effect and other food resource on system stability is studied by numerical simulations.

(1) 

First, we consider the case where there is no fear effect and has the variable other food resource. Set each parameter value of system (11) as $(k,d,a,b)=(0.0,0.3,3.0,0.5)$ and set $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$. Depending on $r_0>2$, $r_0<2$, $h>2$ and $h<2$, we give the corresponding numerical simulations (Figure 12, Figure 13, Figure 14 and Figure 15). We found that in any cases, too many other food sources can accelerate the demise of prey populations. The possible reason is that with the increasing of other food sources, the number of predator populations has increased. Therefore, although the influence of a single predator on the prey is reduced, the overall predator population will have an excessive impact on the prey population, which will lead to the decreasing of the prey population, and finally the prey species will be driven to extinction.

(2) 

Then, we consider the influence of fear effect. Set each parameter value of system (11) as:

$$(d,a,b)=(0.3,3.0,0.5),$$

set $\left(x\right(0),y(0\left)\right)=(0.5,0.3)$, and plot with k as the abscissa. Corresponding to the four cases discussed above, we choose suitable m, such that two species could be coexist, no matter in the stable state or chaos state. Figure 16, Figure 17, Figure 18 and Figure 19 show that in any cases, if the fear effect is too large, the prey species will be driven to extinction, while depending on the intrinsic growth rate of predator species, the predator species will state in stable state for the case $h<2$, and in chaos state for the case $h>2$.

7. Discussion

Based on the traditional Lotka-Volterra predator prey model, Wang, Zanette, and Zou [22] proposed system (5), in which the authors argued that the fear effect of predator to prey may reduce the birth rate of the prey species. Their studied indicated that if the system exists a positive equilibrium, it then is globally asymptotically stable. since the existence condition is independent of the fear effect, one could draw the conclusion that fear effect has no influence to the dynamic behaviors of the system (5). One may conjecture that such a property still hold for the discrete counterpart system, however, Kundu, Pal, Samanta, and Sen [34] showed that for system (9), the situations become complicated. If the birth rate of the prey species is enough large, then, without the influence of fear effect, both predator and prey species may have chaotic behaviors, then, for fixed birth rate, with the increasing of fear effect, the system could becomes stable. They then drew the conclusion “fear effect enhances stability in a predator-prey system”. However, an amazing property about system (9) is that the system may have Hopf bifurcation, as was shown in Theorem 4.1 and Figure 5 in [34]. The existence of Hopf bifurcation exclude the global asymptotically stability of the positive equilibrium, since such kind of phenomenon is induced by fear effect, we think the fear effect also have negative effect on the stability of the system.

Zhu, Wu, Lai and Yu [23] proposed system (6), in which the predator species has other food resource, they showed that the fear effect is one of the most important factors that lead to the extinction of the prey species. The authors of [23] only paid attention to the influence of fear effect, the dynamic behaviors of system (5) and (6) are quite different, the reason may rely in the fact that in system (6), the predator species has other food resource, and this leads to the lower bound of the predator species, and predator populations have lasting effects on prey populations. While in system (5), with the reduce of prey species, the predator species also reduced, since it has no enough food resource. Corresponding to system (6), Chen, He, and Chen [35] proposed system (10). Their studied indicated that the prey free equilibrium may be globally stable under some suitable conditions.

Stimulated by the aforementioned works, we proposed a modified Leslie–Gower discrete predator–prey model in which the prey population has a fear effect and the predator population has other food sources, i.e., system (11). The fundamental difference between model (10) and model (11) is that in model (10), the authors assume that the predator population will increase its intrinsic growth rate after capturing the prey, while in model (11), we assume that the captures the prey will increase the environmental capacity of the predator population.

We discussed the existence, local stability, global attractivity and bifurcation phenomenon of the equilibrium points.

By using iterative method and difference equation comparison principle, we obtained the conclusion that the positive equilibrium point is globally attractive under certain conditions. However, due to scaling in the proving process, the conditions of global attraction of positive equilibrium obtained are more stringent. The results of numerical simulation also show that the global attractivity of positive equilibrium may occur when the theorem conditions are not fixed.

Then, we proved that the boundary equilibrium of prey extinction is also globally attractive under certain conditions, and its stability is related to the birth rate of prey and the intrinsic growth rate of predators. This means that the predator population stabilizes while the prey population becomes extinct, the reason relies in the fact that the predator species have other food sources.

Next, we discussed bifurcations at the equilibrium points of the system. At the same time, through numerical simulation, we found that with the increase of growth rate of predator population, a chaotic pattern appears in both prey and predator populations, and the prey population will eventually become extinct.

Finally, we discussed the effect of fear effect on the stability of system (11). We found that fear effect within a certain range can enhance the stability of the system, and this result was presented intuitively via numerical simulation. However, for the case two species could coexist in a stable state or chaotic state without fear effect, only with the increasing of fear effect, the prey species will finally be driven to extinction. (see Example 5 for more detail). Hence, fear effect of the prey to predator species is one of the essential factor that could lead to the extinction of the prey species.

In this paper, we only considered the fear effect on the first species, and did not consider the influence of functional response. It may be interesting to investigate the discrete predator–prey system with fear effect and functional response; for example, to date, there has been no investigation of the discrete type of system (8). In future, we will try to conduct research in this direction.

At the end of the paper, based on the numeric simulations of this paper, we would like to give two conjectures about the system (11).

Conjecture 1. 

The conditions which ensure the local stability of the prey free boundary equilibrium is enough to ensure its global asymptotical stability.

Conjecture 2. 

Assume that

$${\displaystyle \frac{r_0}{1+k{U}_1^y}-d-b{U}_1^y<2}$$

(51)

and

$$h<2$$

(52)

hold, then system (11) admits a unique positive equilibrium which is globally asymptotically stable.