Effects of prey refuge and predator cooperation on a predator–prey system

ABSTRACT

We develop and investigate a discrete-time predator–prey model with cooperative hunting among predators and a spatial prey refuge. The system can exhibit two positive equilibria if the magnitude of cooperation is large and the maximal reproduction number of predators is less than one. In such a scenario, the predator population may exhibit a strong Allee effect, and therefore the predator may survive if its density is above the threshold. When the positive equilibrium is unique, we prove that hunting cooperation can destabilize the equilibrium through a Neimark–Sacker bifurcation. Numerical findings indicate that a high degree of predator hunting cooperation can help rescue the predator population if the proportion of prey refuge is large, while hunting cooperation becomes destabilizing when the proportion of refuge is small. Despite hunting cooperation's destabilizing effect, it can facilitate predator persistence, particularly with a large proportion of prey refuge. Furthermore, there exists a wide parameter space where the predator–prey interaction may exhibit chaotic behaviour.

KEYWORDS:

• Prey refuge
• cooperative hunting
• Neimark–Sacker bifurcation
• maximal Lyapunov exponents

2020 MSC:

• 92D25
• 92B05

1. Introduction

Cooperative hunting among predators is an important and common biological phenomenon in many ecological communities [16]. The phenomenon is frequently observed in carnivores such as wolves, African wild dogs, lions [13,28,37,38], as well as in other organisms such as spiders, birds and ants [6]. Cooperative hunting has many biological advantages for predator evolution. The mechanism not only increases predation success but also reduces chasing distance and enhances the probabilities of capturing larger prey [16,28]. During recent years, a variety of mathematical models have been developed to study the effects of hunting cooperation. In Ref. [7,8,10,11,33], discrete-time models of cooperative hunting are studied, while models of ordinary differential equations are investigated in Ref. [2,20,23,25,36,39,40,42,48], and reaction diffusion models are explored in Ref. [14,15,21,32,34,44]. Different biological conclusions are reached from these mathematical models depending on modelling assumptions.

On the other hand, prey populations often have access to areas where they are safe from predators [4]. It is widely believed that prey refuge plays a critical role in the survival of prey populations as well as on species composition and diversity [4]. To incorporate spatial refuges into mathematical modelling, one uses either a fixed number or a fixed fraction of prey that is free from predators [22,29]. Undoubtedly, the existence of prey refuges clearly has important consequences on the coexistence of predators and prey. Various mathematical models have been proposed to study the effects of prey refuge on the population dynamics of predator–prey interactions. See [3,24,26,29–31,35,48] and the references cited therein for a list of continuous-time models on prey refuges.

There are many populations in ecosystems that do not reproduce continuously but rather seasonally. Discrete-time models of difference equations are therefore important tools to study populations with non-overlapping generations. For example, Elaydi and Yakubu [17] and Yakubu [45,47] investigate cycles and dispersal in discrete population models, while Friedman and Yakubu [19], Franke and Yakubu [18], and Yakubu [46] use discrete-time models to understand infectious diseases.

To study spatial refuges for populations with non-overlapping generations, Hassell [22] adopts a modified Nicholson-Bailey model with either a fixed number or a fixed proportionality of prey refuge to examine the stabilizing/destabilizing effects of prey refuge. It is concluded in Ref. [22] that proportionate refuges can lead to stability, while constant number of prey refuge cannot. In addition to Ref. [22], several other discrete-time models of prey refuge have been constructed and studied subsequently. See for example [5,9,27,43].

In this work, the goal is to investigate the combined effects of prey refuge and predator hunting cooperation on predator–prey interactions. Since there are many populations that reproduce seasonally, we will use models of difference equations as a framework to explore the dynamical effects of prey refuge and hunting cooperation of predators. In particular, we will verify that the model undergoes a Neimark–Sacker bifurcation at the unique positive equilibrium, and numerical simulations will be carried out to study dynamical consequences.

The remainder of this manuscript is organized as follows. The model is presented in the following section, with mathematical properties of the model provided in the same section. In particular, the results of the model of prey refuge with no hunting cooperation from [9] and the model with no refuge from [10,11] are briefly summarized. Section 3 studies the model numerically, and biological conclusions are given in the final section.

2. Model and analysis

In this section, we present the mathematical model and its analysis. The section is separated into three subsections depending on the magnitude of hunting cooperation.

2.1. Model development and preliminary results

The well-known Beverton-Holt equation is one of the most simple discrete-time models exhibiting only equilibrium dynamics [1]. The model is considered as the discrete analogue of the continuous-time logistic equation. It models contest competition between individuals within the same population when the intra-specific competition involves resources such as mates or food that are stable [28]. Since many mathematical models of ordinary differential equations use logistic growth for the prey population, we adopt the Beverton-Holt equation as a framework to study prey refuge and predator cooperation. As a consequence, any dynamical destabilization would result from interactions with predators.

Let $x_n$ and $y_n$ denote the prey and predator densities at time $n=0,1,2,\dots$, respectively. The parameter $\lambda >0$ is the prey population's intrinsic growth rate whereas k>0 relates to its carrying capacity. In the absence of predators and prey refuge, the evolution of the prey population is governed by the Beverton-Holt equation $x_{n+1}=\frac{\lambda x_n}{1+k{x}_n}$, $n=0,1,2,\dots$ Let γ, $0<\gamma <1$, denote the constant fraction of the prey population that is available to predators. Then $1-\gamma$ is the constant proportion in the refuge patch. All the prey population is in the refuge patch if $\gamma =0$ and there is no prey refuge if $\gamma =1$.

At time n and when there is no predator cooperation, the number of encounters between the prey and predator populations is assumed to follow the law of mass action, $a\left(\gamma x_n\right)y_n$, where a>0 is the constant searching efficiency of the predators and $\gamma x_n$ is the prey population that is available to the predators. It is further assumed that the number of encounters is distributed randomly and follows a Poisson distribution with probability $p\left(m\right)=\frac{{\mathrm{e}}^{-\mu }{\mu }^m}{m!}$, where $m=0,1,2,\dots$ is the number of encounters and µ is the average encounters per prey per time unit. By the given assumptions, $\mu =a\left(\gamma x_n\right)y_n/\left(\gamma x_n\right)=a{y}_n$ and thus $1-p\left(0\right)=1-{\mathrm{e}}^{-a{y}_n}$ is the probability of an individual prey being preyed upon. Let $\text{β}>0$ denote the predator conversion from each prey consumed. The model with no hunting cooperation studied in [9] is given by $\begin{array}{rl}& x_{n+1}=\frac{\lambda (1-\gamma )x_n}{1+k{x}_n}+\frac{\lambda \gamma x_n}{1+k{x}_n}{\mathrm{e}}^{-a{y}_n},\\ & y_{n+1}=\text{β}\gamma x_n(1-{\mathrm{e}}^{-a{y}_n}).\end{array}$Similar to Ref. [10,11], to incorporate predator hunting cooperation, the number of encounters between prey and predators at time n is modified as $a\left(\gamma x_n\right)y_n(1+c{y}_n)$ where $c\ge 0$ is the degree of hunting cooperation. There is no cooperation among predators if c = 0 and the magnitude of cooperation is stronger if c is larger. Applying a similar approach as above, $1-{\mathrm{e}}^{-a{y}_n(1+c{y}_n)}$ becomes the probability of an individual prey being preyed upon. As a result, the interaction between predators and prey is given by the following first-order difference equations (1) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+k{x}_n}+\frac{\lambda \gamma x_n}{1+k{x}_n}{\mathrm{e}}^{-a{y}_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-a{y}_n(1+c{y}_n)}).\end{array}$(1) There are six parameters in (1(1) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+k{x}_n}+\frac{\lambda \gamma x_n}{1+k{x}_n}{\mathrm{e}}^{-a{y}_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-a{y}_n(1+c{y}_n)}).\end{array}$(1) ). To simplify the number of parameters, we let $\hat{x}=kx$, $\hat{y}=ay$, $\widehat{\text{β}}=\frac{\text{β}a}{k}$, $\hat{c}=\frac{c}{a}$, and by ignoring the hats, system (1(1) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+k{x}_n}+\frac{\lambda \gamma x_n}{1+k{x}_n}{\mathrm{e}}^{-a{y}_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-a{y}_n(1+c{y}_n)}).\end{array}$(1) ) is converted to (2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) System (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) in the case of no cooperative hunting, c = 0, is studied in Ref. [9] and the results are summarized below. In particular, using β as the bifurcation parameter, it can be shown that the system exhibits a Neimark–Sacker bifurcation.

Proposition 2.1

Let c = 0. The following statements hold true for (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ).

1. If $\lambda <1$, the extinction equilibrium $E_0=(0,0)$ is globally asymptotically stable.

2. Let $\lambda >1$. Then (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has another equilibrium $E_1=(\lambda -1,0)$.

   1. If $\text{β}\gamma (\lambda -1)<1$, then $E_1$ is globally asymptotically stable.

   2. If $\text{β}\gamma (\lambda -1)>1$, then (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) is uniformly persistent and has a unique positive equilibrium $E_0^{\ast }=(x_0^{\ast },y_0^{\ast })$. Moreover, there exists a unique positive ${\text{β}}_1$ such that $E_0^{\ast }$ undergoes a Neimark–Sacker bifurcation at $\text{β}={\text{β}}_1$.

In the case that there is no prey refuge, $\gamma =1$, system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) is investigated in Ref. [10,11] and we collect several of the results in Proposition 2.2. The analysis of a Neimark–Sacker bifurcation is also centred at the parameter β.

Proposition 2.2

Let $\gamma =1$. The following statements hold for system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ).

1. If $\lambda <1$, then the extinction equilibrium $E_0=(0,0)$ is globally asymptotically stable.

2. Let $\lambda >1$ and $\text{β}\overline{x}>1$. System (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has a unique positive equilibrium and is uniformly persistent. There exists a unique ${\text{β}}_d>0$ such that the unique positive equilibrium undergoes a Neimark–Sacker bifurcation at $\text{β}={\text{β}}_d$.

3. Let $\lambda >1$ and $\text{β}\overline{x}\le 1$. If $0<2c\le \frac{3\lambda -1}{\lambda -1}$, then (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has no positive equilibrium and there are at most two positive equilibria if $2c>\frac{3\lambda -1}{\lambda -1}$.

For model (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) in the general setting, it is clear that its solutions exist, remain nonnegative and are bounded for $n=0,1,2,\dots$ Moreover, the extinction equilibrium $E_0=(0,0)$ always exists and from the Jacobian matrix of (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) at $E_0$, $E_0$ is asymptotically stable if $\lambda <1$ and a saddle point if $\lambda >1$. For the latter, the stable manifold lies on the nonnegative y-axis. By a simple comparison method, it can be proven that $E_0$ is globally asymptotically stable if $\lambda <1$ and globally attracting if $\lambda =1$. In addition, system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has the following properties listed in Proposition 2.3.

Proposition 2.3

The following statements hold for system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ).

1. If $\lambda <1$, then $E_0=(0,0)$ is globally asymptotically stable for (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) and $E_0$ is globally attracting if $\lambda =1$.

2. If $\lambda >1$, then $E_0$ is unstable and there exists another boundary equilibrium $E_1=(\overline{x},0)$, $\overline{x}=\lambda -1$, where $E_1$ is locally asymptotically stable if $\text{β}\gamma \overline{x}<1$ and unstable if $\text{β}\gamma \overline{x}>1$.

3. If $\lambda >1$ and $\text{β}\gamma \overline{x}>1$, then (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has a unique positive equilibrium if $2c\le 1$. There are no positive equilibria if $2c\le 1$ and $\text{β}\gamma \overline{x}<1$.

Proof.

The proof of (a) is trivial by a comparison argument. To verify (b) and (c), notice that $E_1=(\overline{x},0)$ is biologically feasible, where $\overline{x}=\lambda -1$. The Jacobian matrix evaluated at $E_1$ has the form $\left(\begin{array}{cc}1/\lambda & \ast \\ 0& \text{β}\gamma \overline{x}\end{array}\right)$, where * is an unimportant entry. Therefore, $E_1$ is asymptotically stable if $\text{β}\gamma \overline{x}<1$ and a saddle point with stable manifold on the positive x-axis if $\text{β}\gamma \overline{x}>1$.

For the existence of a positive equilibrium, notice that the nontrivial x- and y-isoclines are given respectively by (3) $\begin{array}{rl}& x=\lambda (1-\gamma )-1+\lambda \gamma {\mathrm{e}}^{-y(1+cy)}:=f\left(y\right),\\ & x=\frac{y}{\text{β}\gamma (1-{\mathrm{e}}^{-y(1+cy)})}:=h\left(y\right),\end{array}$(3) where $f\left(0\right)=\overline{x},f\left(\mathrm{\infty }\right)=\lambda (1-\gamma )-1,f^{\mathrm{\prime }}\left(y\right)=-\lambda \gamma {\mathrm{e}}^{-y(1+cy)}(1+2cy)<0\text{for}y\ge 0.$If $\lambda (1-\gamma )-1\ge 0$, then $f\left(y\right)>0$ for $y\ge 0$ while if $\lambda (1-\gamma )-1<0$ there exists (4) $y_r=\frac{-1+\sqrt{1+4c\ln (\frac{\lambda \gamma }{1-\lambda (1-\gamma )}})}{2c}>0$(4) such that $f\left(y_r\right)=0$, $f\left(y\right)>0$ on $[0,y_r)$ and $f\left(y\right)<0$ for $y>y_r$. For $h\left(y\right)$, we take advantage of a similar function studied in [9]. In particular, $h\left(0\right)=\frac{1}{\text{β}\gamma },h\left(\mathrm{\infty }\right)=\mathrm{\infty },h^{\mathrm{\prime }}\left(y\right)=\frac{{\mathrm{e}}^{-y(1+cy)}g\left(y\right)}{\text{β}\gamma (1-{\mathrm{e}}^{-y(1+cy)}{)}^2},$where $g\left(y\right)={\mathrm{e}}^{y(1+2cy)}-1-y(1+2cy)$. One can conclude that $h^{\mathrm{\prime }}\left(y\right)>0$ for $y\ge 0$ if $2c\le 1$, and there exists a unique critical point $\overline{y}>0$ such that $h^{\mathrm{\prime }}\left(y\right)<0$ on $[0,\overline{y})$ while $h^{\mathrm{\prime }}\left(y\right)>0$ on $[\overline{y},\mathrm{\infty })$ if 2c>1. Consequently, if $2c\le 1$, then (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has a unique positive equilibrium $E^{\ast }=(x^{\ast },y^{\ast })$ if $\text{β}\gamma \overline{x}>1$ and there is no positive equilibrium if $\text{β}\gamma \overline{x}<1$.

In view of Proposition 2.3, we assume $\lambda >1$ for the remainder of this section and separate it into the cases of $2c\le 1$ and 2c>1.

2.2. The case of $2c\le 1$

Let $2c\le 1$. It can be shown as in the proof of Theorem 4.2 of [11] that $E_1$ is globally asymptotically stable if $\text{β}\gamma \overline{x}<1$. Let $\text{β}\gamma \overline{x}>1$ so that the unique positive equilibrium $E^{\ast }=(x^{\ast },y^{\ast })$ exists by Proposition 2.3(c). The Jacobian matrix evaluated at $E^{\ast }$ is given by $J\left(E^{\ast }\right)=\left(\begin{array}{cc}\hat{a}& \hat{b}\\ \hat{c}& \hat{d}\end{array}\right)$, where (5) $\begin{array}{rl}& \hat{a}=\frac{1}{1+x^{\ast }},\hat{b}=-x^{\ast }(1-\frac{\lambda (1-\gamma )}{1+x^{\ast }})(1+2c{y}^{\ast }),\\ & \hat{c}=\frac{y^{\ast }}{x^{\ast }},\hat{d}=\frac{y^{\ast }{\mathrm{e}}^{-y^{\ast }(1+c{y}^{\ast })}(1+2c{y}^{\ast })}{1-{\mathrm{e}}^{-y^{\ast }(1+c{y}^{\ast })}}.\end{array}$(5) The following theorem summarizes the stability of $E^{\ast }$.

Theorem 2.1

Let $\lambda >1$ and $2c\le 1$. If $\text{β}\gamma \overline{x}<1$ then $E_1=(\overline{x},0)$ is globally asymptotically stable in Γ. If $\text{β}\gamma \overline{x}>1$, then $E_1$ is unstable (a saddle point) and (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has a unique positive equilibrium $E^{\ast }=(x^{\ast },y^{\ast })$ with $0<tr\left(J\right(E^{\ast }\left)\right)<1+det\left(J\right(E^{\ast }\left)\right)$ and $det\left(J\right(E^{\ast }\left)\right)>0$. Moreover, $E^{\ast }$ is asymptotically stable (stable focus) if $det\left(J\right(E^{\ast }\left)\right)<1$ and unstable (unstable focus) if $det\left(J\right(E^{\ast }\left)\right)>1$.

Proof.

The global stability of $E_1$ can be proven similarly as in [11, Theorem 4.2] when $\text{β}\gamma \overline{x}<1$. To understand the stability of $E^{\ast }$, we apply the Jury conditions [1]. Specifically, $E^{\ast }$ is locally asymptotically stable if $|tr\left(J\right(E^{\ast }\left)\right)|<1+det\left(J\right(E^{\ast }\left)\right)<2$, where $tr\left(J\right)$ and $det\left(J\right)$ are the trace and determinant of the matrix J, respectively. Clearly, $0<\hat{a}<1$, $\hat{c}>0$ and $\hat{d}>0$. Also, $\hat{b}<0$ since $1-\frac{\lambda (1-\gamma )}{1+x^{\ast }}=\frac{\lambda \gamma }{1+x^{\ast }}{\mathrm{e}}^{-y^{\ast }(1+c{y}^{\ast })}>0$ from the first equation of (3(3) $\begin{array}{rl}& x=\lambda (1-\gamma )-1+\lambda \gamma {\mathrm{e}}^{-y(1+cy)}:=f\left(y\right),\\ & x=\frac{y}{\text{β}\gamma (1-{\mathrm{e}}^{-y(1+cy)})}:=h\left(y\right),\end{array}$(3) ). Hence, $det\left(J\right(E^{\ast }\left)\right)>0$ and $tr\left(J\right(E^{\ast }\left)\right)>0$. Moreover, $tr\left(J\right(E^{\ast }\left)\right)<1+det\left(J\right(E^{\ast }\left)\right)$ is equivalent to $(1-\hat{a})(1-\hat{d})-\hat{b}\hat{c}>0$, where $\hat{d}<1$ if and only if $c(y^{\ast }{)}^2<\frac{\left(y^{\ast }{)}^2\right(1+c{y}^{\ast }{)}^2}{2}+\frac{\left(y^{\ast }{)}^3\right(1+c{y}^{\ast }{)}^3}{6}+\cdots$, which is trivially true since $2c\le 1$, i.e. $tr\left(J\right(E^{\ast }\left)\right)<1+det\left(J\right(E^{\ast }\left)\right)$ holds. Therefore, $E^{\ast }$ is asymptotically stable if $det\left(J\right(E^{\ast }\left)\right)<1$ and unstable if $det\left(J\right(E^{\ast }\left)\right)>1$.

Figure 1 presents components of the positive equilibrium $E^{\ast }$ in (a), (c) and (e) whereas (b), (d) and (f) plot $det\left(J\right(E^{\ast }\left)\right)$ against c for $0\le c\le 0.5$. In this example, we see that a large proportion of refuge, $1-\gamma$, not only increases the equilibrium level of both populations but also stabilizes the positive equilibrium $E^{\ast }$.

Figure 1. The left and right panels plot components of the unique positive equilibrium and the determinant of $J\left(E^{\ast }\right)$, respectively. Fixed parameter values are $\lambda =10$ and $\text{β}=3$ with c varying between 0 and 0.5. The solid and dashed curves in (a), (c), and (e) denote the x and y components of the positive equilibrium, respectively.

Figure 2 illustrates the destabilization effect of hunting cooperation when there is no prey refuge, $\gamma =1$. Other parameter values are $\lambda =10$, $\text{β}=0.3$ with c given by 0.4 and 0.5 in (a) and (b) respectively. Here, we use a smaller β value since Figure 1(f) shows that $E^{\ast }$ is unstable for $0\le c\le 0.5$ if β is large. Although not presented, there is only equilibrium dynamics when c = 0. As c is increased to c = 0.4, the equilibrium $E^{\ast }$ becomes unstable and is marked as a black circle in Figure 2(a,b). The unstable $E^{\ast }$ is enclosed by the invariant circle, where the closed invariant curves are generated using the initial condition $(4,0.6)$ by eliminating the first 5000 iterations and plotting the next 1000 iterations.

Figure 2. Phase portrait along with the unstable positive equilibrium for system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) are plotted with fixed parameter values of $\lambda =10$, $\gamma =1$, and $\text{β}=0.3$. The degree of hunting cooperation is c = 0.4 in (a) and c = 0.5 in (b). These demonstrate that hunting cooperation destabilizes the predator–prey interaction when there is no prey refuge.

It is expected that $E^{\ast }$ undergoes a Neimark–Sacker bifurcation at $det\left(J\right(E^{\ast }\left)\right)=1$ by Theorem 2.1. We now verify the Neimark–Sacker bifurcation of model (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) at the unique positive equilibrium $E^{\ast }$, where c is assumed to be the bifurcation parameter with $c\le 0.5$ under the condition $\text{β}\gamma \overline{x}>1$. The characteristic equation of the Jacobian matrix $J\left(E^{\ast }\right)$ has two complex conjugate roots with modulus one if $det\left(J\right(E^{\ast }\left)\right)=1$ which is equivalent to $\hat{a}\hat{d}-\hat{b}\hat{c}=1$. Therefore, $E^{\ast }$ may undergo a Neimark–Sacker bifurcation if c obeys the following equation: $\frac{\text{β}\gamma x^{\ast }{\mathrm{e}}^{-y^{\ast }(1+c{y}^{\ast })}(1+2c{y}^{\ast })}{(1+x^{\ast })}+y^{\ast }(1-\frac{\lambda (1-\gamma )}{1+x^{\ast }})(1+2c{y}^{\ast })=1.$

Theorem 2.2

Let $\text{β}\gamma \overline{x}>1$ and $2c\le 1$. If (9(9) $p\left(0\right)=\frac{1}{1+x^{\ast }}+[\text{β}\gamma x^{\ast }-y^{\ast }][1+2c_1{y}^{\ast }]\ne 0,1.$(9) ) is fulfilled and $\theta \ne 0$, then model (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) exhibits a Neimark–Sacker bifurcation at the unique positive equilibrium point $E^{\ast }=(x^{\ast },y^{\ast })$ when the parameter c is varied within a small neighbourhood of HB. In addition, an attracting (resp., repelling) invariant closed curve bifurcates from $E^{\ast }$ for $c>c_1$ (resp., $c<c_1$) if $\theta <0$ (resp., $\theta >0$).

Proof.

We define the following bifurcation set: $\begin{array}{rl}HB& =\{(c,\lambda ,\text{β},\gamma )\in {\mathbb{R}}_{+}^4:c\le 0.5,\frac{\text{β}\gamma x^{\ast }{\mathrm{e}}^{-y^{\ast }(1+c{y}^{\ast })}(1+2c{y}^{\ast })}{(1+x^{\ast })}\\ & +y^{\ast }(1-\frac{\lambda (1-\gamma )}{1+x^{\ast }})(1+2c{y}^{\ast })=1\}.\end{array}$A Neimark–Sacker bifurcation occurs at the unique positive equilibrium point $E^{\ast }$ when c is varied around a small neighbourhood of HB. We consider model (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) with parameters $(c_1,\lambda ,\text{β},\gamma )$, where the parameters $(c_1,\lambda ,\text{β},\gamma )$ are chosen randomly from the set HB, and introduce a perturbation of $c_1$ as follows: (6) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n[1+(c_1+c^{\ast }\left)y_n\right]},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n[1+(c_1+c^{\ast }\left)y_n\right]}),\end{array}$(6) where $|c^{\ast }|\ll 1$ is a small perturbation parameter.

Let $u_n=x_n-x^{\ast }$, and $v_n=y_n-y^{\ast }$. Then model (6(6) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n[1+(c_1+c^{\ast }\left)y_n\right]},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n[1+(c_1+c^{\ast }\left)y_n\right]}),\end{array}$(6) ) is turned into (7) $\begin{array}{rl}{u}_{n+1}& =a_{10}{u}_n+a_{01}{v}_n+a_{20}{u}_n^2+a_{11}{u}_n{v}_n+a_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\\ v_{n+1}& =b_{10}{u}_n+b_{01}{v}_n+b_{11}{u}_n{v}_n+b_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\end{array}$(7) where $\{\begin{array}{rl}{a}_{10}& =\frac{1}{1+x^{\ast }},a_{01}=-\frac{\lambda (\text{β}\gamma x^{\ast }-y^{\ast })}{\text{β}(1+x^{\ast })}[1+2(c_1+c^{\ast }\left)y^{\ast }\right],a_{20}=-\frac{1}{(1+x^{\ast }{)}^2},\\ a_{11}& =-\frac{\lambda (\text{β}\gamma x^{\ast }-y^{\ast })}{\text{β}{x}^{\ast }(1+x^{\ast }{)}^2}[1+2(c_1+c^{\ast }\left)y^{\ast }\right],\\ a_{02}& =\frac{\lambda (\text{β}\gamma x^{\ast }-y^{\ast })}{2\text{β}(1+x^{\ast })}[1-2(c_1+c^{\ast })+4(c_1+c^{\ast })y^{\ast }+4(c_1+c^{\ast }{)}^2{y}^2],\\ b_{10}& =\frac{y^{\ast }}{x^{\ast }},b_{01}=[\text{β}\gamma x^{\ast }-y^{\ast }][1+2(c_1+c^{\ast }\left)y^{\ast }\right],\\ b_{11}& =\frac{\text{β}\gamma x^{\ast }-y^{\ast }}{x^{\ast }}[1+2(c_1+c^{\ast }\left)y^{\ast }\right],\\ b_{02}& =-(\text{β}\gamma x^{\ast }-y^{\ast })[1-2(c_1+c^{\ast })+4(c_1+c^{\ast })y^{\ast }+4(c_1+c^{\ast }{)}^2{y}^2].\end{array}$The linearized model of (7(7) $\begin{array}{rl}{u}_{n+1}& =a_{10}{u}_n+a_{01}{v}_n+a_{20}{u}_n^2+a_{11}{u}_n{v}_n+a_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\\ v_{n+1}& =b_{10}{u}_n+b_{01}{v}_n+b_{11}{u}_n{v}_n+b_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\end{array}$(7) ) at $(u,v)=(0,0)$ has the following characteristic equation: (8) ${\eta }^2-p\left(c^{\ast }\right)\eta +q\left(c^{\ast }\right)=0,$(8) where $\begin{array}{rl}p\left(c^{\ast }\right)& =\frac{1}{1+x^{\ast }}+[\text{β}\gamma x^{\ast }-y^{\ast }][1+2(c_1+c^{\ast }\left)y^{\ast }\right],\\ q\left(c^{\ast }\right)& =\frac{(\text{β}{x}^{\ast }+\lambda y^{\ast })}{\text{β}{x}^{\ast }(1+x^{\ast })}[\text{β}\gamma x^{\ast }-y^{\ast }][1+2(c_1+c^{\ast }\left)y^{\ast }\right].\end{array}$The roots of (8(8) ${\eta }^2-p\left(c^{\ast }\right)\eta +q\left(c^{\ast }\right)=0,$(8) ) at $c^{\ast }=0$ are complex conjugate numbers η and $\overline{\eta }$ with modulus 1 as the parameters $(c_1,\lambda ,\text{β},\gamma )\in HB$, where $\eta ,\overline{\eta }=\frac{p\left(c^{\ast }\right)}{2}\pm \frac{i}{2}\sqrt{4q\left(c^{\ast }\right)-p^2\left(c^{\ast }\right)}.$Therefore, we have $|\eta {|}_{c^{\ast }=0}=\sqrt{q\left(0\right)}=1,\text{and}l=\frac{d|\eta |}{d{c}^{\ast }}{|}_{c^{\ast }=0}=y^{\ast }\sqrt{\frac{(\text{β}{x}^{\ast }+\lambda y^{\ast })(\text{β}\gamma x^{\ast }-y^{\ast })}{\left[\text{β}{x}^{\ast }\right(1+x^{\ast }\left)\right](1+2c_1{y}^{\ast })}}>0.$Moreover, for a Neimark–Sacker bifurcation to occur at $c^{\ast }=0$, it is necessary that ${\eta }^m,{\overline{\eta }}^m\ne 1(m=1,2,3,4)$, which is the same as $p\left(0\right)\ne -2,0,1,2$. Note that $p\left(0\right)\ne \pm 2$ is inferred from the fact that $(c_1,\lambda ,\text{β},\gamma )\in HB$. Thus, it is sufficient to require that $p\left(0\right)\ne 0,1$, which results in (9) $p\left(0\right)=\frac{1}{1+x^{\ast }}+[\text{β}\gamma x^{\ast }-y^{\ast }][1+2c_1{y}^{\ast }]\ne 0,1.$(9) Hence, when $c^{\ast }=0$ and the criteria (9(9) $p\left(0\right)=\frac{1}{1+x^{\ast }}+[\text{β}\gamma x^{\ast }-y^{\ast }][1+2c_1{y}^{\ast }]\ne 0,1.$(9) ) hold, the eigenvalues $\eta ,\overline{\eta }$ do not lie near the point where the unit circle and the coordinate axes cross.

Let $\mu =\frac{p\left(0\right)}{2}$ and $\omega =\frac{1}{2}\sqrt{4q\left(0\right)-p^2\left(0\right)}$. Then the following transformation yields the normal form of (7(7) $\begin{array}{rl}{u}_{n+1}& =a_{10}{u}_n+a_{01}{v}_n+a_{20}{u}_n^2+a_{11}{u}_n{v}_n+a_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\\ v_{n+1}& =b_{10}{u}_n+b_{01}{v}_n+b_{11}{u}_n{v}_n+b_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\end{array}$(7) ) at $c^{\ast }=0$: $\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{cc}{a}_{01}& 0\\ \mu -a_{10}& -\omega \end{array}\right)\left(\begin{array}{c}\tilde{x}\\ \tilde{y}\end{array}\right).$Hence, model (7(7) $\begin{array}{rl}{u}_{n+1}& =a_{10}{u}_n+a_{01}{v}_n+a_{20}{u}_n^2+a_{11}{u}_n{v}_n+a_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\\ v_{n+1}& =b_{10}{u}_n+b_{01}{v}_n+b_{11}{u}_n{v}_n+b_{02}{v}_n^2+O\left(\right(|u_n|+|v_n|{)}^3),\end{array}$(7) ) changes to $\left(\begin{array}{c}{\overline{x}}_{n+1}\\ {\overline{y}}_{n+1}\end{array}\right)=\left(\begin{array}{cc}\mu & -\omega \\ \omega & \mu \end{array}\right)\left(\begin{array}{c}{\tilde{x}}_n\\ {\tilde{y}}_n\end{array}\right)+\left(\begin{array}{c}{\tilde{f}}_1({\tilde{x}}_n,{\tilde{y}}_n)\\ {\tilde{f}}_2({\tilde{x}}_n,{\tilde{y}}_n)\end{array}\right),$where $\begin{array}{rl}{\tilde{f}}_1({\tilde{x}}_n,{\tilde{y}}_n)& =\frac{a_{20}}{a_{01}}{u}_n^2+\frac{a_{11}}{a_{01}}{u}_n{v}_n+\frac{a_{02}}{a_{01}}{v}^2+O\left(\right(|{\tilde{x}}_n|+|{\tilde{y}}_n|{)}^3),\\ {\tilde{f}}_2({\tilde{x}}_n,{\tilde{y}}_n)& =\frac{(\mu -a_{10})a_{20}}{a_{01}\omega }{u}_n^2+\frac{(\mu -a_{10})a_{11}-a_{01}{b}_{11}}{a_{01}\omega }{u}_n{v}_n\\ & +\frac{(\mu -a_{10})a_{02}-a_{01}{b}_{02}}{a_{01}\omega }{v}_n^2+O\left(\right(|{\tilde{x}}_n|+|{\tilde{y}}_n|{)}^3),\end{array}$and $\begin{array}{rl}{u}_n^2& =a_{01}^2{\tilde{x}}_n^2,u_n{v}_n=a_{01}(\mu -a_{10}){\tilde{x}}_n^2-a_{01}\omega {\tilde{x}}_n{\tilde{y}}_n,\\ v_n^2& =(\mu -a_{10}{)}^2{\tilde{x}}_n^2-2\omega (\mu -a_{10}){\tilde{x}}_n{\tilde{y}}_n+{\omega }^2{\tilde{y}}_n^2.\end{array}$To continue, let us define the next nonzero real number: (10) $\theta ={[-\mathrm{\Re }\left(\frac{(1-2\eta ){\overline{\eta }}^2}{1-\eta }{L}_{11}{L}_{12}\right)-\frac{1}{2}|L_{11}{|}^2-|L_{21}^2|+\mathrm{\Re }(\overline{\eta }{L}_{22}\left)\right]}_{c^{\ast }=0},$(10) where $\begin{array}{rl}{L}_{11}& =\frac{1}{4}\left[\right({\tilde{f}}_{1{\tilde{x}}_n{\tilde{x}}_n}+{\tilde{f}}_{1{\tilde{y}}_n{\tilde{y}}_n})+i({\tilde{f}}_{2{\tilde{x}}_n{\tilde{x}}_n}+{\tilde{f}}_{2{\tilde{y}}_n{\tilde{y}}_n}\left)\right],\\ L_{12}& =\frac{1}{8}\left[\right({\tilde{f}}_{1{\tilde{x}}_n{\tilde{x}}_n}-{\tilde{f}}_{1{\tilde{y}}_n{\tilde{y}}_n}+2{\tilde{f}}_{2{\tilde{x}}_n{\tilde{y}}_n})+i({\tilde{f}}_{2{\tilde{x}}_n{\tilde{x}}_n}-{\tilde{f}}_{2{\tilde{y}}_n{\tilde{y}}_n}-2{\tilde{f}}_{1{\tilde{x}}_n{\tilde{y}}_n}\left)\right],\\ L_{21}& =\frac{1}{8}\left[\right({\tilde{f}}_{1{\tilde{x}}_n{\tilde{x}}_n}-{\tilde{f}}_{1{\tilde{y}}_n{\tilde{y}}_n}-2{\tilde{f}}_{2{\tilde{x}}_n{\tilde{y}}_n})+i({\tilde{f}}_{2{\tilde{x}}_n{\tilde{x}}_n}-{\tilde{f}}_{2{\tilde{y}}_n{\tilde{y}}_n}+2{\tilde{f}}_{1{\tilde{x}}_n{\tilde{y}}_n}\left)\right],\\ L_{22}& =\frac{1}{16}\left[\right({\tilde{f}}_{1{\tilde{x}}_n{\tilde{x}}_n{\tilde{x}}_n}+{\tilde{f}}_{1{\tilde{x}}_n{\tilde{y}}_n{\tilde{y}}_n}+{\tilde{f}}_{2{\tilde{x}}_n{\tilde{x}}_n{\tilde{y}}_n}+{\tilde{f}}_{2{\tilde{y}}_n{\tilde{y}}_n{\tilde{y}}_n})\\ & +i({\tilde{f}}_{2{\tilde{x}}_n{\tilde{x}}_n{\tilde{x}}_n}+{\tilde{f}}_{2{\tilde{x}}_n{\tilde{y}}_n{\tilde{y}}_n}-{\tilde{f}}_{1{\tilde{x}}_n{\tilde{x}}_n{\tilde{y}}_n}-{\tilde{f}}_{1{\tilde{y}}_n{\tilde{y}}_n{\tilde{y}}_n})].\end{array}$

2.3. The case of 2c>1

It is clear that (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has a unique positive equilibrium $E^{\ast }=(x^{\ast },y^{\ast })$ if $\text{β}\gamma \overline{x}>1$ and 2c>1. In this scenario, since $tr\left(J\right(E^{\ast }\left)\right)<1+det\left(J\right(E^{\ast }\left)\right)$ may or may not hold true, the proof of Theorem 2.2 with some modifications can be adopted to show the existence of a Neimark–Sacker bifurcation at $E^{\ast }$. It is suspected that $E^{\ast }$ is globally asymptotically stable whenever it is locally asymptotically stable. However, this conjecture is not easy to verify. Instead, we show that system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) is uniformly persistent using a similar argument as in the proof of Ref. [11, Theorem 2.2]. It follows that both prey and predator populations can coexist if $\text{β}\gamma \overline{x}>1$ even with a large fraction of prey refuge as long as β or λ is large.

Theorem 2.3

Let $\lambda >1$, $\text{β}\gamma \overline{x}>1$ and 2c>1. Then system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has a unique positive equilibrium and is uniformly persistent.

Suppose now $\text{β}\gamma \overline{x}<1$. Notice $\overline{x}$ is the prey's carrying capacity and hence $\text{β}\gamma \overline{x}$ can be interpreted as the maximal reproduction number of the predator. When this reproductive number is less than one, the extinction equilibrium $E_0=(0,0)$ is locally asymptotically stable and the predator goes extinct if both populations have sufficiently small densities. Therefore, (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) is not uniformly persistent. Moreover, the number of positive equilibria can not be determined from the graphs of the isoclines directly. There can be either zero, one or two positive equilibria. When $\gamma =1$, it is shown in Ref. [11, Theorem 3.1] that there exists no positive equilibrium if $2c\le \frac{3\lambda -1}{\lambda -1}$ and $\text{β}\gamma \overline{x}\le 1$ whereas there are at most two positive equilibria if $2c>\frac{3\lambda -1}{\lambda -1}$ and $\text{β}\gamma \overline{x}\le 1$. See Proposition 2.2(c). For $0<\gamma <1$ we can apply a similar argument as in Ref. [11, Theorem 3.1]. Indeed, setting $f\left(y\right)=h\left(y\right)$, results in (11) $\text{β}\gamma (1-{\mathrm{e}}^{-y(1+cy)})\left(\lambda \right(1-2\gamma )-1+2\lambda \gamma {\mathrm{e}}^{-y(1+cy)})=y.$(11) Let the left hand sides of (11(11) $\text{β}\gamma (1-{\mathrm{e}}^{-y(1+cy)})\left(\lambda \right(1-2\gamma )-1+2\lambda \gamma {\mathrm{e}}^{-y(1+cy)})=y.$(11) ) be denoted by $w\left(y\right)$. Then $w^{\mathrm{\prime }}\left(y\right)=\text{β}\gamma {\mathrm{e}}^{=y(1+cy)}(1+2cy)\left(\lambda \right(1-2\gamma )-1+2\lambda \gamma {\mathrm{e}}^{-y(1+cy)}),$and $w^{\mathrm{\prime }\mathrm{\prime }}\left(y\right)=\text{β}\gamma {\mathrm{e}}^{-2y(1+cy)}u\left(y\right),$where $u\left(y\right)=(2c-(1+2cy{)}^2\left)\right(\lambda (1-2\gamma -1){\mathrm{e}}^{-y(1+cy)}+4\lambda \gamma c-4\lambda \gamma (1+2cy{)}^2.$Observe that $u\left(0\right)=2c(\lambda -1)+1-\lambda -2\lambda \gamma >0$ if and only if $2c>\frac{\lambda (2\gamma +1)-1}{\lambda -1}$. By similar computations and arguments as in [11, Theorem 3.1], we obtain the following.

Theorem 2.4

Let $\lambda >1$ and $\text{β}\gamma \overline{x}\le 1$. Then system (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has no positive equilibrium if $1<2c\le \frac{\lambda (2\gamma +1)-1}{\lambda -1}$ and there are at most two positive equilibria if $2c>\frac{\lambda (2\gamma +1)-1}{\lambda -1}$.

The statement of Theorem 2.4 reduces to Proposition 2.2(c) when $\gamma =1$. Figure 3(a) provides a bifurcation diagram with respect to c when $\lambda =3$, $\text{β}=4$ and $\gamma =0.1$ so that $\text{β}\gamma \overline{x}=0.8<1$. In this example, $\frac{\lambda (2\gamma +1)-1}{\lambda -1}=1.3$ and there is no positive equilibrium if $c\le 0.65$ by Theorem 2.4. Even though we only present the simulation result for up to c = 80, the following observation remains true for $c\in [0,1000]$ simulated. There always exists an asymptotically stable equilibrium $E_0=(0,0)$ for $c\ge 0$, which is marked by the black solid line in Figure 3(a). A saddle node bifurcation occurs at $c\approx 2.241$ (not proven), beyond which there are two positive equilibria $E_i^{\ast }=(x_i^{\ast },y_i^{\ast })$, i = 1, 2, $y_1^{\ast }<y_2^{\ast }$, and below which there are no positive equilibria. The upper branch equilibrium $E_2^{\ast }$ (the red solid curve) is asymptotically stable whereas the lower branch equilibrium $E_1^{\ast }$, marked as the blue dashed curve, is a saddle point for c>2.241. As a result, bistability occurs for c>2.241. The predator density in the stable equilibrium increases with increasing c and there is no observed Neimark–Sacker bifurcation as c is increased to 1000. In this parameter space, the predator exhibits a strong Allee effect [12] due to its large degree of cooperation. The predator population can either persist or go to extinction depending on its density.

Figure 3. Fixed parameter values are $\lambda =3$, $\text{β}=4$ and $\gamma =0.1$. A bifurcation diagram with respect to c is given in (a) and (b)–(c) present basins of attraction of $E_2^{\ast }$ (in red) and of $E_0$ (in cyan). The degree of hunting cooperation is c = 5 in (b) and c = 8 in (c).

With the same parameter values of λ, β and γ, Figure 3(b,c) presents basins of attraction of $E_2^{\ast }$ marked in red and of $E_0$ shown in cyan colour with c = 5 in (b) and c = 8 in (c). Here, 50, 000 randomly chosen initial conditions are used to simulate the basin. The equilibrium points are $E_1^{\ast }=(1.976,0.064)$ and $E_2^{\ast }=(1.720,0.642)$ in Figure 3(b) whereas they are $E_1^{\ast }=(1.986,0.036)$ and $E_2^{\ast }=(1.704,0.672)$ for c = 8 in Figure 3(c). It is apparent that as c is increased from 5 to 8 the basin of attraction of $E_2^{\ast }$ is enlarged, i.e. the Allee threshold of the predator is smaller as the degree of cooperation increases. Therefore, cooperative hunting in this parameter space has the effect of promoting predator persistence.

Let (12) $\mathrm{\Gamma }=\left\{\right(x,y)\in {\mathbb{R}}_{+}^2:x>0\}.$(12) In the case that $\text{β}\gamma \overline{x}$ is small enough, it can be shown that $E_1=(\overline{x},0)$ is globally asymptotically stable in Γ. The proof is similar to the proof of Ref. [11, Theorem 4.4] and is therefore omitted.

Theorem 2.5

Let $\lambda >1$, 2c>1 and $\text{β}\gamma \overline{x}\sqrt{2c}\left({\mathrm{e}}^{\frac{1-2c}{4c}}\right)<1$. Then $E_1=(\overline{x},0)$ is globally asymptotically stable in Γ.

Since $\sqrt{2c}\left({\mathrm{e}}^{\frac{1-2c}{4c}}\right)>1$ for 2c>1, the assumption $\text{β}\gamma \overline{x}\sqrt{2c}\left({\mathrm{e}}^{\frac{1-2c}{4c}}\right)<1$ in Theorem 2.5 requires that $\text{β}\gamma \overline{x}$ is sufficiently smaller than one, and Theorem 2.5 implies that the model (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ) has no positive equilibrium in this parameter regime.

3. Numerical simulations

In this section, the goal is to use numerical techniques to investigate the predator–prey interactions encapsulated by model (2(2) $\begin{array}{rl}{x}_{n+1}& =\frac{\lambda (1-\gamma )x_n}{1+x_n}+\frac{\lambda \gamma x_n}{1+x_n}{\mathrm{e}}^{-y_n(1+c{y}_n)},\\ y_{n+1}& =\text{β}\gamma x_n(1-{\mathrm{e}}^{-y_n(1+c{y}_n)}),\\ & x_0,y_0\ge 0.\end{array}$(2) ). In particular, we wish to understand the effects of prey refuge and hunting cooperation of predators. Unless otherwise stated, parameter values $\lambda =10$ and $\text{β}=0.3$ are fixed. In this circumstance, the interaction only possesses equilibrium dynamics when there is no cooperative hunting c = 0 and no prey refuge $\gamma =1$. See Figure 5(a,b) to be presented below.

Figure 4. Bifurcation diagrams with respect to c for $0\le c\le 5$ are plotted. Fixed parameter values are $\lambda =10$ and $\text{β}=0.3$. The value of γ is 0.9 in (a)–(b), 0.6 in (c)–(d) and 0.3 in (e)–(f).

Figure 5. Bifurcation diagrams with respect to γ for $0\le \gamma \le 1$ are plotted. Fixed parameter values are $\lambda =10$ and $\text{β}=0.3$. The parameter value of c is 0 in (a)–(b), 0.9 in (c)–(d), and 25 in (e)–(f).

We first study the effects of cooperative hunting relative to population interactions. Bifurcation diagrams with respect to c for $0\le c\le 5$ are presented in Figure 4 with different γ values. Specifically, $\gamma =0.9$ in Figure 4(a,b), $\gamma =0.6$ in Figure 4(c,d), and $\gamma =0.3$ in Figure 4(e,f). From the figure we see that prey refuge can stabilize the predator–prey interaction since as γ is decreased to 0.6 or 0.3 the system converges to equilibria for $0\le c\le 5$. This conclusion is consistent with that given for Figure 1. Moreover, when there is a large proportion of prey population in the refuge patch, large hunting cooperation can promote predator persistence as shown in Figure 4(e,f).

Next, bifurcation diagrams with respect to γ for $0\le \gamma \le 1$ are presented in Figure 5 with $\lambda =10$ and $\text{β}=0.3$, where c = 0 in Figure 5(a,b), c = 0.9 in (c)–(d), and c = 25 in (e)–(f). Notice that if $\gamma <0.374$ then $\text{β}\gamma \overline{x}<1$ for the given λ and β values. We see from this figure that if the proportion of prey refuge $1-\gamma$ is large, i.e. γ is small, the predator population cannot survive for all of the different values of c simulated. Moreover, a small proportion of prey refuge is destabilizing to the interaction, shown in Figure 5(c–f). Also, as c is increased to c = 25, the predator population can persist under a smaller γ value. Therefore, hunting cooperation can promote predator persistence. Further, comparing the plots in Figure 5 we reach the conclusion that hunting cooperation is destabilizing to predator–prey interaction, which is also concluded in Figures 4.

We now investigate the effects of parameter β using bifurcation diagrams presented in Figure 6 with fixed $\lambda =2$ and c = 0.2. The other parameter γ is chosen as 0.9, 0.7 and 0.6 in Figure 6(a–f), respectively. It can be seen that large predator conversion is destabilizing in Figure 6(a–d). However, when the proportion $1-\gamma$ of prey refuge is increased to 0.4, the predator–prey interaction exhibits only equilibrium dynamics and increasing the predator conversion can increase the predator population level at the equilibrium, shown in Figure 6(e,f). As γ is decreased to 0.4, although not presented, the predator is extinct if β is small while the interaction is at a positive equilibrium if β is larger.

Figure 6. Bifurcation diagrams with respect to β are plotted. Fixed parameter values are $\lambda =2$ and c = 0.2. The other parameter is $\gamma =0.9$ in (a)–(b), $\gamma =0.7$ in (c)–(d), and $\gamma =0.6$ in (e)–(f).

Finally, we approximate maximal Lyapunov exponents with respect to different parameters to detect possible chaotic interactions. Our numerical method is adopted from [41] with the aid of Matlab. Figure 7 presents the Lyapunov exponents against parameter c in (a)–(b) and versus γ in (c)–(d), where $\lambda =10$ and $\text{β}=0.3$. When $\gamma =0.6$ in (b) the Lyapunov exponents are below 0 for $0\le c\le 5$ implying there is only equilibrium behaviour. On the other hand, when $\gamma =0.9$ in (a) the Lyapunov exponents are positive for some c values between 2 and 3, and between 4 and 5 indicating sensitive dependence on initial conditions and therefore possible chaotic dynamics occur in this parameter space. As a consequence, increasing hunting cooperation can destabilize the dynamics and lead to chaotic interaction when the proportion of prey refuge is small, $1-\gamma =0.1$. Figure 7(c,d) simulates Lyapunov exponents against γ where c = 2.5 in (c) and c = 0 in (d). In particular, Figure 7(d) suggests that there is likely a transcritical bifurcation at γ approximately equal to 0.4. Furthermore, Figure 7(c,d) implies that the dynamics are very likely to be chaotic when γ is large, i.e. when the proportion of prey refuge is small. Therefore, prey refuge has a stabilizing effect on the predator–prey interaction.

Figure 7. Maximal Lyapunov exponents with respect to c in (a)–(b) and γ in (c)–(d) are approximated. Fixed parameter values are $\lambda =10$ and $\text{β}=0.3$. Plots (e) and (f) present the maximal Lyapunov exponents with respect to β with $\gamma =0.9$ in (e) and $\gamma =0.7$ in (f). The other parameter values are $\lambda =2$ and c = 0.2.

Figure 7(e,f) presents the maximal Lyapunov exponents with respect to β for fixed $\lambda =2$ and c = 0.2. The parameter γ is 0.9 in (e) and 0.7 in (f). When the proportion of prey refuge is small in (e), the system is chaotic in a large range of β values between 5 and 8. However, as the proportion of prey refuge is increased in (f), the dynamics are not chaotic for all β between 2 and 8. We again conclude that prey refuge can be stabilizing to the predator–prey interaction.

4. Conclusions

In this work, we utilized a simple model of nonlinear difference equations to examine the impact of prey refuge and predator hunting cooperation on the dynamics of predator–prey interactions. In the absence of predators, the prey exhibits only equilibrium dynamics. Therefore, any non-equilibrium behaviour arises as a result of interactions with predators.

We have proven that the unique positive equilibrium exhibits a Neimark–Sacker bifurcation when the magnitude of hunting cooperation is varied at a small critical value, provided that the maximal reproduction number of predators exceeds one. We conducted numerical simulations to investigate a scenario in which the predator–prey system is at the positive equilibrium in the absence of prey refuge and hunting cooperation. In this situation, a small magnitude of cooperative hunting destabilizes the system when the proportion of prey refuge is small. Conversely, a large degree of cooperation can promote predator persistence when the proportion of prey refuge is large. When the proportion of prey refuge is small, the populations remain at the positive equilibrium without hunting cooperation. However, increasing the degree of hunting cooperation leads to a Neimark–Sacker bifurcation in the system. As the proportion of prey refuge increases, the predator cannot survive without hunting cooperation or with only a small degree of cooperation. On the other hand, large cooperation in this case can promote predator persistence. These conclusions are supported by Figures 4, 5, and 6. Regarding predator conversion, a high turnover of predators destabilizes the interaction in one scenario, but can stabilize the interaction when the proportion of prey refuge is large.

To detect chaotic interactions, we numerically approximated maximal Lyapunov exponents with respect to several model parameters. It was found that the model does show sensitive dependence on initial conditions, indicating the presence of potential chaos in the interactions. Consequently, utilizing predators as biological control agents to constrain the prey population, especially if the prey is a pest, would pose greater challenges.

We will now compare the model results obtained here with those given in Ref. [9–11]. First, both populations go extinct if $\lambda <1$ for all $c\ge 0$ and $0\le \gamma \le 1$. If $\lambda >1$ and the maximal reproduction number of predators is greater than one, i.e. $\text{β}\gamma \overline{x}=\text{β}\gamma (\lambda -1)>1$, then both populations can coexist for all positive initial conditions, given $0<\gamma \le 1$ satisfying the inequality and $c\ge 0$. Additionally, the model has a unique positive equilibrium, and the destabilization of the equilibrium occurs via a Neimark–Sacker bifurcation. This local bifurcation has been verified based on the parameter β in Ref. [9,11], whereas in the present study, it is varied with respect to c. If $\text{β}\gamma \overline{x}<1$, then the predator-extinction equilibrium $(\overline{x},0)$ is globally asymptotically stable if there is no hunting cooperation (c = 0) or if the degree of cooperation is small with $c\le 0.5$. With a high degree of hunting cooperation, the model can exhibit two positive equilibria, and the predator may display a strong Allee effect, where the predator population can persist if its densities are above the Allee threshold.

Acknowledgments

We thank one of the reviewers for the valuable suggestions and comments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.