Transcritical bifurcation and Neimark-Sacker bifurcation of a discrete predator-prey model with herd behaviour and square root functional response

ABSTRACT

In this paper, a discrete predator-prey model incorporating herd behaviour and square root response function is deduced from its continuous version by the semi-discretization method. Firstly, the existence and local stability of fixed points of the system are studied by applying a key lemma. Secondly, by employing the centre manifold theorem and bifurcation theory, the conditions for the occurrences of the transcritical bifurcation and Neimark-Sacker bifurcation are obtained. Not only that but the direction and stability conditions of the bifurcated closed orbits are also clearly shown. Finally, numerical simulations are also given to confirm the existence of Neimark-Sacker bifurcation.

KEYWORDS:

• Predator-prey system with herd behaviour
• square root functional response
• semi-discretization method
• transcritical bifurcation
• Neimark-Sacker bifurcation

1. Introduction

Since the pioneering work of Lotka and Volterra in the 1920s, the study of realistic mathematical models in ecology, especially for the relationship between species and environment, has become a hot issue that attracted flocks of famous mathematicians and biologists. Over the past several decades, the predator-prey models have been widely studied [1–11]. Recently, Huang and Ruan [12] considered the classical Gause type predator-prey model in the following form:

(1.1) $\left\{\begin{array}{c}\frac{dx}{dt}=xg(x,k)-yp(x),\\ \frac{dy}{dt}=y(-s+cq(x)),\end{array}\right.$(1.1)

where $x(t)$ and $y(t)$ represent the population densities of prey and predator in time $t$ respectively, $g(x,k)$ describes the specific growth rate of the prey in the absence of predator, the parameters $c$ and $s$ are the positive constants, respectively, indicating the efficiency of predator in converting consumed prey into their growth and predator mortality rate. The functional response $p(x)$ of predator to the prey describes the variation in the density of the prey attacked by per predator per unit time as the prey density changes. The function $q(x)$ describes how the predator converts the consumed prey into the growth of predator [13–16].

Due to the realistic meaning and in order to show the crowding effect [17,18], when the prey is large, the prey growth rate $g(x,k)$ in the model (1.1) is usually a negative value. One can assume that the average growth rate of a typical prey species is the logistic form $xg(x,k)=rx(1-\frac{x}{K})$, among them, the positive constants $r$ and $K$ respectively represent the inherent growth rate of the prey and the carrying capacity of environment to the prey without the predator. Hence, the system (1.1) can be written as

(1.2) $\left\{\begin{array}{c}\frac{dx}{dt}=rx(1-\frac{x}{K})-yp(x),\\ \frac{dy}{dt}=y(-s+cq(x)).\end{array}\right.$(1.2)

Generally speaking, there are many types of the functional response $g(x,y)$ [19,20]. These predator-prey systems with prey groups’ defence ability have more abundant and interesting dynamic characteristics and attract attention of many scholars [21,22]. Recently, an article [23] models the fact that it is the individuals at the edge of the herd that generally suffer the heaviest consequences of the predators’ attacks. Some predator- prey models in which the prey exhibits herd behaviour were considered. Bian et al. proposed a system [22] with the Beddington-DeAngelis functional response; De Assis et al. proposed a system [21] with the square-root functional response and so on. As a mathematical consequence of the herd behaviour, the interaction terms in a system use the square root of the prey population rather than simply the prey population. The use of the square-root properly accounts for the assumption that the interactions occur along the boundary of the population. In this paper, one talks about the following system analysed by Braza [7] where $p(x)=q(x)=m\sqrt{x}$:

(1.3) $\left\{\begin{array}{c}\frac{dx}{dt}=rx(1-\frac{x}{K})-m\sqrt{x}y,\\ \frac{dy}{dt}=y(-s+cm\sqrt{x}),\end{array}\right.$(1.3)

here, the biological significance of the positive parameters $c$, $s$, $r$ and $K$ are as described above. The parameter $m$ is also positive, which describes the search efficiency coefficient of the predator for the prey.

Although the authors obtained some good results for the system (1.3), some problems remain. In an ordinary way, it is impossible to obtain an exact solution for a complex differential equation system. Therefore, many researchers derive the corresponding approximate solution by computer. What the computer considers is the discrete points, so it is very natural and practical to discretize the model. For a given continuous system, there are many different discrete methods, including the Euler forward difference scheme, the Euler backward difference scheme, the semi-discretization method, etc. [24–28]. In this paper, we use the semi-discretization method.

First, to simplify the later calculations, by letting $\frac{x}{K}\to x$, $\frac{my}{rK}\to y$, $rt\to t$, $\sqrt{K}\to a$, $\frac{s}{r}\to d$, $\frac{cm\sqrt{K}}{r}\to b$, we obtain a simpler form of the system (1.3) as follows:

(1.4) $\left\{\begin{array}{c}\frac{dx}{dt}=x(1-x)-a\sqrt{x}y,\\ \frac{dy}{dt}=y(-d+b\sqrt{x}),\end{array}\right.$(1.4)

where $a$, $b$ and $d$ are all positive constants.

Next, we adopt the semi-discretization method, which does not need to consider the step size [20,27], to derive its discrete model. For this, suppose that $[t]$ denotes the greatest integer not exceeding $t$. Consider the average change rate of the system (1.4) at integer number points

(1.5) $\left\{\begin{array}{c}\frac{1}{x(t)}\frac{\mathrm{d}x(t)}{\mathrm{d}t}=1-x([t])-\frac{ay([t])}{\sqrt{x([t])}},\\ \frac{1}{y(t)}\frac{\mathrm{d}y(t)}{\mathrm{d}t}=-d+b\sqrt{x([t])}.\end{array}\right.$(1.5)

It is very straightforward to see that piecewise constant arguments occur in the system (1.5) and that any solution $(x(t),y(t))$ of (1.5) for $t\in [0,+\mathrm{\infty })$ is in possession of the following three characteristics:

1. $x(t)$ and $y(t)$ are continuous on the interval $[0,+\mathrm{\infty })$;

2. when $t\in [0,+\mathrm{\infty })$ possibly except for the points $t\in \{0,1,2,3,\cdots \}$,$\frac{\mathrm{d}x(t)}{\mathrm{d}t}$ and $\frac{\mathrm{d}y(t)}{\mathrm{d}t}$ exist anywhere;

3. the system (1.5) is true in any interval $[n,n+1)$ with $n=0,1,2,3,\cdots$.

The following system can be derived by integrating the system (1.5) in the interval $[n,t]$ for any $t\in [n,n+1)$ and $n=0,1,2,\cdots$:

(1.6) $\left\{\begin{array}{c}x(t)=x_n{e}^{1-x_n-\frac{a{y}_n}{\sqrt{x_n}}}(t-n),\\ y(t)=y_n{e}^{-d+b{x}_n}(t-n),\end{array}\right.$(1.6)

where $x_n=x(n)$ and $y_n=y(n)$.

Letting $t\to (n+1{)}^{-}$ in the system (1.6) produces

(1.7) $\left\{\begin{array}{c}{x}_{n+1}=x_n{e}^{1-x_n-\frac{a{y}_n}{\sqrt{x_n}}},\\ y_{n+1}=y_n{e}^{-d+b\sqrt{x_n}},\end{array}\right.$(1.7)

where all the parameters

$(a,b,d)\in \mathrm{\Omega }=\{(a,b,d)\in R_{+}^3|a>0,b>0,d>0\}.$

In this paper, our main goal is to consider the dynamics of the system (1.7), mainly for its bifurcation problems besides its stability. So as to study the stability and local bifurcation of the fixed points, we need a definition and a key lemma [13–16,27], which are included in the appendix for readers’ convenience.

The rest of this paper is organized as follows. In Section 2, the existence and local stability of the fixed points of the system (1.7) are studied. In Section 3, we derive the sufficient conditions for the occurrences of transcritical bifurcation and Neimark-Sacker bifurcation of the system (1.7). In Section 4, numerical simulations are performed to confirm the above theoretical results. Finally, this paper is ended with some discussions and conclusions in Section 5.

2. The existence and stability of fixed points

In this section, we will not only study the existence of fixed points of the system (1.7) but also obtain the local stability of these fixed points.

The fixed points of the system (1.7) satisfy

$x=x{e}^{1-x-\frac{ay}{\sqrt{x}}},y=y{e}^{-d+b\sqrt{x}}.$

Given the biological meanings of the system (1.7), we only take into consideration its nonnegative fixed points. Thus, the system (1.7) has and only has two nonnegative fixed points $E_1(1,0)$ and $E_2(x_0,y_0)$ for $b>d$, where

$x_0=\frac{d^2}{b^2},y_0=\frac{1}{a}(1-\frac{d^2}{b^2})\frac{d}{b}.$

The Jacobian matrix of the system (1.7) at any fixed point $E(x,y)$ takes the following form:

$J(E)=\left(\begin{array}{cc}[1+x(-1+\frac{1}{2}a{x}^{-\frac{3}{2}}y)]e^{1-x-\frac{ay}{\sqrt{x}}}& -a\sqrt{x}{e}^{1-x-\frac{ay}{\sqrt{x}}}\\ \frac{1}{2}b{x}^{-\frac{1}{2}}y{e}^{-d+b\sqrt{x}}& e^{-d+b\sqrt{x}}\end{array}\right).$

The characteristic polynomial of Jacobian matrix $J(E)$ reads

$F(\lambda )={\lambda }^2-p\lambda +q,$

where

$p=\mathrm{t}\mathrm{r}J(E)=[1+x(-1+\frac{1}{2}a{x}^{-\frac{3}{2}}y)]e^{1-x-\frac{ay}{\sqrt{x}}}+e^{-d+b\sqrt{x}},$

$q=detJ(E)=[1+x(-1+\frac{1}{2}a{x}^{-\frac{3}{2}}y)+\frac{1}{2}aby]e^{1-x-\frac{ay}{\sqrt{x}}-d+b\sqrt{x}}.$

By applying Definition 5.1 and Lemma 5.2 in the appendix, we can easily come to the following Theorems 2.1, 2.2, respectively, for the stability of fixed points $E_1$ and $E_2$.

Theorem 2.1.

The following statements about the fixed point $E_1(1,0)$ of the system (1.7) are true.

1. If $b>d$, then $E_1$ is a saddle.

2. If $b=d$, then $E_1$ is non-hyperbolic.

3. If $b<d$, then $E_1$ is a stable node, i.e. a sink.

Proof. The Jacobian matrix of the system (1.7) at the fixed point $E_1$ can be simplified as follows:

$J(E_1)=\left(\begin{array}{cc}0& -a\\ 0& e^{b-d}\end{array}\right).$

Obviously, ${\lambda }_1=0$ and $|{\lambda }_2|=e^{b-d}$.

Note $|{\lambda }_1|<1$ is always true. If $b>d$, then $|{\lambda }_2|=e^{b-d}>1$, so $E_1$ is a saddle; If $b=d$, then $|{\lambda }_2|=e^{b-d}=1$, implying $E_1$ is non-hyperbolic; If $b<d$, then $|{\lambda }_2|=e^{b-d}<1$, therefore $E_1$ is a stable node, namely, a sink. The proof is complete.

Theorem 2.2.

When $b>d$, $E_2(\frac{d^2}{b^2},\frac{1}{a}(1-\frac{d^2}{b^2})\frac{d}{b})$ is a positive fixed point of the system (1.7). Let $b_1\stackrel{\mathrm{\Delta }}{=}d\sqrt{\frac{d+3}{d+1}}$, then the following statements are true about the fixed point $E_2$ of the system (1.7).

1. If $d<b<b_1$, then $E_2$ is a stable node, i.e. a sink.

2. If $b=b_1$, then $E_2$ is non-hyperbolic.

3. If $b>b_1$, then $E_2$ is an unstable node, i.e. a source.

Proof. The Jacobian matrix of the system (1.7) at the fixed point $E_2$ can be simplified as follows:

$J(E_2)=\left(\begin{array}{cc}\frac{3}{2}(1-\frac{d^2}{b^2})& -\frac{ad}{b}\\ \frac{1}{2}\frac{b}{a}(1-\frac{d^2}{b^2})& 1\end{array}\right).$

Hereout we obtain the characteristic polynomial of Jacobian matrix $J(E_2)$, which reads as

(2.1) $F(\lambda )={\lambda }^2-p\lambda +q,$(2.1)

where

$p=\frac{5}{2}-\frac{3}{2}\frac{d^2}{b^2},q=(\frac{3}{2}+\frac{d}{2})(1-\frac{d^2}{b^2}).$

For $b>d$ we have $p=\frac{5}{2}-\frac{3}{2}\frac{d^2}{b^2}>1$, $q>0$ and

$p<2\iff \frac{5}{2}-\frac{3}{2}\frac{d^2}{b^2}<2\iff b<\sqrt{3}d,$

$d<b_1=d\sqrt{\frac{d+3}{d+1}}<\sqrt{3}d\iff d>0,$

$\left\{\begin{array}{c}q<1\\ q=1\\ q>1\end{array}\iff \begin{array}{c}b<b_1\\ b=b_1\\ b>b_1\end{array}.\right.$

Moreover,

$F(1)=\frac{d}{2}(1-\frac{d^2}{b^2})>0,F(-1)=2+3(1-\frac{d^2}{b^2})+\frac{d}{2}(1-\frac{d^2}{b^2})>0.$

So, by applying Lemma 5.2, one has:

1. if $d<b<b_1$, then $|{\lambda }_1|<1$ and $|{\lambda }_2|<1$, $E_2$ is a stable node, i.e. a sink;

2. if $b=b_1$, then Eq. (2.1) has a pair of conjugate complex roots ${\lambda }_1$ and ${\lambda }_2$ with $|{\lambda }_1|=|{\lambda }_2|=1$, implying $E_2$ is non-hyperbolic;

3. if $b>b_1$, then $|{\lambda }_1|>1$ and $|{\lambda }_2|>1$, meaning $E_2$ is an unstable node, i.e. a source.

The proof is finished.

3. Bifurcation problems and their proofs

In this section, we apply the Center Manifold Theorem and local bifurcation theory to principally analyse the local bifurcation problems of the system (1.7) at the fixed points $E_1(1,0)$ and $E_2(\frac{d^2}{b^2},\frac{1}{a}(1-\frac{d^2}{b^2})\frac{d}{b}$ ($b>d$), considering the practical biological meaning.

3.1. For the fixed point $E_1(1,0)$

We know from Eq. (1.4) that no matter what values the parameters $a$ $b$ and $d$ take, the fixed point $E_1$ always exists. Theorem 2.1 indicates that a bifurcation of the system (1.7) may occur at the fixed points $E_1(1,0)$ in the space of parameters $\mathrm{\Omega }$.

Theorem 3.1.

Consider the system (1.7). Set the parameters $(a,b,d)\in \mathrm{\Omega }$. Take $b_0=d$, then the system (1.7) undergoes a transcritical bifurcation at the fixed point $E_1$ when the parameter $b$ varies in a small neighbourhood of the critical value $b_0$.

Proof. In order to demonstrate the detailed process, we adopt the steps as follows.

Step one. Transform the fixed point $E_1$ to the origin $O(0,0)$ by taking the changes of variables $u_n=x_n-1$, $v_n=y_n-0$. Then, the system (1.7) is changed to

(3.1) $\left\{\begin{array}{c}{u}_{n+1}=(u_n+1)e^{-u_n-\frac{a{v}_n}{\sqrt{u_n+1}}}-1,\\ v_{n+1}=v_n{e}^{-d+b\sqrt{u_n+1}}.\end{array}\right.$(3.1)

Step two. Giving a small perturbation $b^{\ast }=b-b_0$ with $0<|b^{\ast }|\ll 1$, of the parameter $b$ around the critical value $b_0$, the system (3.1) is perturbed into

(3.2) $\left\{\begin{array}{c}{u}_{n+1}=(u_n+1)e^{-u_n-\frac{a{v}_n}{\sqrt{u_n+1}}}-1,\\ v_{n+1}=v_n{e}^{-d+(b^{\ast }+b_0)\sqrt{u_n+1}}.\end{array}\right.$(3.2)

Letting $b_{n+1}^{\ast }=b_n^{\ast }=b^{\ast }$, the system (3.2) can be written as

(3.3) $\left\{\begin{array}{c}{u}_{n+1}=(u_n+1)e^{-u_n-\frac{a{v}_n}{\sqrt{u_n+1}}}-1,\\ v_{n+1}=v_n{e}^{-d+(b_n^{\ast }+b_0)\sqrt{u_n+1}},\\ b_{n+1}^{\ast }=b_n^{\ast }.\end{array}\right.$(3.3)

Step three. Taylor expanding of the system (3.3) at $(u_n,v_n,b_n^{\ast })=(0,0,0)$ to the third-order term can be written as

(3.4) $\left(\begin{array}{c}{u}_{n+1}\\ v_{n+1}\\ b_{n+1}^{\ast }\end{array}\right)=\left(\begin{array}{ccc}0& -a& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{u}_n\\ v_n\\ b_n^{\ast }\end{array}\right)+\left(\begin{array}{c}{g}_1(u_n,v_n,b_n^{\ast })+o({\rho }_1^3)\\ g_2(u_n,v_n,b_n^{\ast })+o({\rho }_1^3)\\ 0\end{array}\right),$(3.4)

where ${\rho }_1=\sqrt{u_n^2+v_n^2+b{{ }_n^{\ast }}^2}$,

$g_1(u_n,v_n,b_n^{\ast })=-\frac{1}{2}{u}_n^2+\frac{a^2}{2}{v}_n^2+\frac{a}{2}{u}_n{v}_n+\frac{1}{3}{u}_n^3$

$-\frac{a^3}{6}{v}_n^3+\frac{a}{8}{u}_n^2{v}_n-\frac{a^2}{2}{u}_n{v}_n^2,$

$g_2(u_n,v_n,b_n^{\ast })=\frac{d}{2}{u}_n{v}_n+v_n{b}_n^{\ast }-\frac{d-d^2}{8}{u}_n^2{v}_n$

$+\frac{1}{2}{v}_{n}b{ }_n^{\ast }2+\frac{d+1}{2}{u}_n{v}_n{b}_n^{\ast }.$

It is easy to derive the three eigenvalues of the matrix

$A=\left(\begin{array}{ccc}0& -a& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$

to be ${\lambda }_1=0$ and ${\lambda }_{2,3}=1$ with corresponding eigenvectors

${\xi }_1=(1,0,0{)}^T,{\xi }_2=(-a,1,0{)}^T,{\xi }_3=(0,0,1{)}^T.$

Step four. Set $T=\left({\xi }_1,{\xi }_2,{\xi }_3\right),$ i.e.

$T=\left(\begin{array}{ccc}1& -a& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$

then

$T^{-1}=\left(\begin{array}{ccc}1& a& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right).$

Using transformation

$(u_n,v_n,b_n^{\ast }{)}^T=T(X_n,Y_n,{\delta }_n{)}^T,$

the system (3.4) is transformed into the following form

(3.5) $\left(\begin{array}{c}{X}_{n+1}\\ Y_{n+1}\\ {\delta }_{n+1}\end{array}\right)=\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{X}_n\\ Y_n\\ {\delta }_n\end{array}\right)+\left(\begin{array}{c}{g}_3(X_n,Y_n,{\delta }_n)+o({\rho }_2^3)\\ g_4(X_n,Y_n,{\delta }_n)+o({\rho }_2^3)\\ 0\end{array}\right),$(3.5)

where ${\rho }_2=\sqrt{X_n^2+Y_n^2+{\delta }_n^2}$,

$g_3(X_n,Y_n,{\delta }_n)=g_1(X_n-a{Y}_n,Y_n,{\delta }_n)+a{g}_2(X_n-a{Y}_n,Y_n,{\delta }_n),$

$g_4(X_n,Y_n,{\delta }_n)=g_2(X_n-a{Y}_n,Y_n,{\delta }_n).$

Step five. Suppose on the centre manifold

$X_n=h(Y_n,{\delta }_n)=a_{20}{Y}_n^2+a_{11}{Y}_n{\delta }_n+a_{02}{\delta }_n^2+o({\rho }_3^2),$

where ${\rho }_3=\sqrt{Y_n^2+{\delta }_n^2}$, then,

$h(Y_{n+1},{\delta }_{n+1})=a_{20}{Y}_{n+1}^2+a_{11}{Y}_{n+1}{\delta }_{n+1}+a_{02}{\delta }_{n+1}^2+o({\rho }_3^2)$

$=a_{20}[Y_n+g_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n){]}^2$

$+a_{11}[Y_n+g_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)]{\delta }_n+a_{02}{\delta }_n^2+o({\rho }_3^2),$

$X_{n+1}=g_3(h(Y_n,{\delta }_n),Y_n,{\delta }_n)+o({\rho }_2^3)$

$=g_1(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)$

$+a{g}_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)+o({\rho }_3^2),$

and $X_{n+1}=h(Y_{n+1},{\delta }_{n+1})$, so we obtain the centre manifold equation

$g_1(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)+a{g}_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)+o({\rho }_3^2)$

$=a_{20}[Y_n+g_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n){]}^2$

$+a_{11}[Y_n+g_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)]{\delta }_n+a_{02}{\delta }_n^2+o({\rho }_3^2).$

Comparing the corresponding coefficients of terms with the same orders in the above centre manifold equation, we get

$a_{20}=-\frac{a^2(d+1)}{2},a_{11}=a,a_{02}=0.$

So, the system (3.5) restricted to the centre manifold takes as

$Y_{n+1}=f_1(Y_n,{\delta }_n):=Y_n+g_2(h(Y_n,{\delta }_n)-a{Y}_n,Y_n,{\delta }_n)+o({\rho }_3^2)$

$=Y_n-\frac{ad}{2}{Y}_n^2+Y_n{\delta }_n+o({\rho }_3^2).$

It is not difficult to calculate

$f_1(0,0)=0,\frac{\mathrm{\partial }{f}_1(0,0)}{\mathrm{\partial }{Y}_n}=1,\frac{\mathrm{\partial }{f}_1(0,0)}{\mathrm{\partial }{\delta }_n}=0,$

$\frac{{\mathrm{\partial }}^2{f}_1(0,0)}{\mathrm{\partial }{Y}_n{\delta }_n}=1\ne 0,\frac{{\mathrm{\partial }}^2{f}_1(0,0)}{\mathrm{\partial }{Y}_n^2}=-ad\ne 0.$

According to (21.1.43)-(21.1.46) in Wiggins [29], all the conditions for the occurrence of a transcritical bifurcation are established. There is no doubt that a transcritical bifurcation occurs in the fixed point $E_1$. The proof is complete.

3.2. For the fixed point $E_2(\frac{d^2}{b^2},\frac{1}{a}(1-\frac{d^2}{b^2})\frac{d}{b})$

When $b=b_1=d\sqrt{\frac{d+3}{d+1}}$ and $b>d$, Theorem 2.2 and Lemma 5.2 show that $F(1)>0$, $F(-1)>0$, $q=1$ and $-2<p<2$, therefore, ${\lambda }_1$ and ${\lambda }_2$ are a pair of conjugate complex roots with $|{\lambda }_1|=|{\lambda }_2|=1$. In this case, we derive that the system (1.7) at the fixed point $E_2$ can undergo a Neimark-Sacker bifurcation in the space of parameters $(a,b,d)\in S_2=\{(a,b,d)\in R_{+}^3|a>0,b>d>0\}$.

To express this process more clearly, let us take the following steps.

The first step. Transform the fixed point $E_2$ to the origin $O(0,0)$ by taking the changes of variables $u_n=x_n-x_0$, $v_n=y_n-y_0$. Then, the system (1.7) is changed into the following form

(3.6) $\left\{\begin{array}{c}{u}_{n+1}=(u_n+x_0)e^{1-u_n-x_0-\frac{a(v_n+y_0)}{\sqrt{u_n+x_0}}}-x_0,\\ v_{n+1}=(v_n+y_0)e^{-d+b\sqrt{u_n+x_0}}-y_0.\end{array}\right.$(3.6)

The second step. Give a small perturbation $b^{\ast }=b-b_1$ with $0<|b^{\ast }|\ll 1$, of the parameter $b$ around the critical value $b_1$, then the system (3.6) is perturbed into

(3.7) $\left\{\begin{array}{c}{u}_{n+1}=(u_n+x_0)e^{1-u_n-x_0-\frac{a(v_n+y_0)}{\sqrt{u_n+x_0}}}-x_0,\\ v_{n+1}=(v_n+y_0)e^{-d+(b^{\ast }+b_1)\sqrt{u_n+x_0}}-y_0.\end{array}\right.$(3.7)

Denote the characteristic equation of the linearized equation of the system (3.7) at the coordinate origin point (0,0) as

$F(\lambda )={\lambda }^2-p(b^{\ast })\lambda +q(b^{\ast })=0,$

where

$p(b^{\ast })=\frac{5}{2}-\frac{3}{2}\frac{d^2}{{(b^{\ast }+b_1)}^2},$

and

$q(b^{\ast })=(\frac{3}{2}+\frac{d}{2})(1-\frac{d^2}{{(b^{\ast }+b_1)}^2}).$

It is not hard to find $p^2(0)-4q(0)<0$, then one can obtain that the two roots of $F(\lambda )=0$ are

${\lambda }_{1,2}(b^{\ast })=\omega \pm \mu i=\frac{p(b^{\ast })\pm \sqrt{p^2(b^{\ast })-4q(b^{\ast })}}{2}=\frac{p(b^{\ast })\pm i\sqrt{4q(b^{\ast })-p^2(b^{\ast })}}{2},$

where $\omega =\frac{p(b^{\ast })}{2}$, $\mu =\frac{1}{2}\sqrt{4q(b^{\ast })-p^2(b^{\ast })}.$

Not only that, there are

${\lambda }_{1,2}(0)=\frac{(d+6)\pm i\sqrt{3d^2+12d}}{2(d+3)}.$

Moreover,

$(|{\lambda }_{1,2}(b^{\ast })|){|}_{b^{\ast }=0}=\sqrt{q(b^{\ast })}{|}_{b^{\ast }=0}$

$=\sqrt{(\frac{3}{2}+\frac{d}{2})(1-\frac{d^2}{b_1^2})}=1,$

which means

${\left.(\frac{d|{\lambda }_{1,2}(b^{\ast })|}{d{b}^{\ast }})\right|}_{b^{\ast }=0}=-(\frac{3}{2}+\frac{d}{2})(\frac{d+1}{d+3}{)}^{\frac{3}{2}}<0.$

Since $p(b^{\ast }){|}_{b^{\ast }=0}=\frac{5}{2}-\frac{3}{2}\frac{d^2}{b_1^2}$ and $q(b^{\ast }){|}_{b^{\ast }=0}=1$, we have

${\lambda }_{1,2}(0)=\frac{(d+6)\pm i\sqrt{3d^2+12d}}{2(d+3)},$

hence it is easy to derive ${\lambda }_{1,2}^m(0)\ne 1$ for all $m=1,2,3,4$.

Thus, we know that the following two conditions for the occurrence of Neimark-Sacker bifurcation are satisfied:

$(A1){\left.\frac{\mathrm{d}|{\lambda }_{1,2}(b^{\ast })|}{\mathrm{d}{b}^{\ast }}\right|}_{b^{\ast }=0}\ne 0;$

$(A2){\lambda }_{1,2}^i(0)\ne 1,i=1,2,3,4.$

The third step. In order to discover the normal form of the system (3.7) when $b^{\ast }=0$, we reshape the system (3.7) in Taylor expansion to the third-order form around the origin $(0,0)$ as follows

(3.8) $\left(\begin{array}{c}{u}_{n+1}\\ v_{n+1}\end{array}\right)=\left(\begin{array}{cc}\frac{3}{d+3}& -a\sqrt{\frac{d+1}{d+3}}\\ \frac{d}{a\sqrt{(d+1)(d+3)}}& 1\end{array}\right)\left(\begin{array}{c}{u}_n\\ v_n\end{array}\right)+\left(\begin{array}{c}{g}_5(u_n,v_n)+o({\rho }_4^3)\\ g_6(u_n,v_n)+o({\rho }_4^3)\end{array}\right),$(3.8)

where ${\rho }_4=\sqrt{u_n^2+v_n^2}$,

$g_5(u_n,v_n)=a_{20}{u}_n^2+a_{11}{u}_n{v}_n+a_{02}{v}_n^2$

$+a_{30}{u}_n^3+a_{21}{u}_n^2{v}_n+a_{12}{u}_n{v}_n^2+a_{03}{v}_n^3,$

$g_6(u_n,v_n)=b_{20}{u}_n^2+b_{11}{u}_n{v}_n+b_{02}{v}_n^2$

$+b_{30}{u}_n^3+b_{21}{u}_n^2{v}_n+b_{12}{u}_n{v}_n^2+b_{03}{v}_n^3,$

$a_{20}=-\frac{2d^2+15d+9}{4(d+1)(d+3)},a_{11}=-a(\frac{1}{2}\sqrt{\frac{d+3}{d+1}}-\frac{d}{\sqrt{(d+1)(d+3)}}),$

$a_{02}=\frac{a^2}{2},a_{30}=\frac{8d^3+51d^2+36d-27}{24(d+{1)}^2(d+3)},$

$a_{21}=a(-\frac{{(d+3)}^{\frac{3}{2}}}{8(d+{1)}^{\frac{3}{2}}}+\frac{3d{(d+3)}^{\frac{1}{2}}}{4(d+{1)}^{\frac{3}{2}}}-\frac{4d^2-3d-9}{8(d+{1)}^{\frac{3}{2}}{(d+3)}^{\frac{1}{2}}}+\frac{3(d+{1)}^{\frac{1}{2}}}{8d{(d+3)}^{\frac{3}{2}}}),$

$a_{12}=-\frac{3}{4}{a}^2,a_{03}=-\frac{a^3}{6}\sqrt{\frac{d+3}{d+1}},$

$b_{20}=-\frac{(d-d^2)(d+{3)}^{\frac{1}{2}}}{4a{(1+d)}^{\frac{3}{2}}},b_{11}=\frac{d(d+3)}{2(d+1)},b_{02}=b_{12}=b_{03}=0,$

$b_{30}=\frac{1}{a}\frac{{(d+3)}^{\frac{3}{2}}}{{(d+1)}^{\frac{5}{2}}}(\frac{d-d^2}{8}+\frac{d^3}{24}),b_{21}=\frac{d^2-d}{8}(\frac{d+3}{d+1}{)}^2.$

It is easy to figure out the two eigenvalues of the matrix

$B=\left(\begin{array}{cc}\frac{3}{d+3}& -a\sqrt{\frac{d+1}{d+3}}\\ \frac{d}{a\sqrt{(d+1)(d+3)}}& 1\end{array}\right)=\left(\begin{array}{cc}{a}_{10}& a_{01}\\ b_{10}& b_{01}\end{array}\right)$

to be ${\lambda }_1(0)$ and ${\lambda }_2(0)$.

The fourth step. Take matrix

$T=\left(\begin{array}{cc}0& a_{01}\\ \mu & 1-\omega \end{array}\right),then{T}^{-1}=\left(\begin{array}{cc}\frac{\omega -1}{\mu a_{01}}& \frac{1}{\mu }\\ \frac{1}{a_{01}}& 0\end{array}\right).$

Using the transformation

$(u,v{)}^T=T(X,Y{)}^T,$

the system (3.8) is transformed into the following form

(3.9) $\left(\begin{array}{c}X\\ Y\end{array}\right)\to \left(\begin{array}{cc}\omega & -\mu \\ \mu & \omega \end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)+\left(\begin{array}{c}F(X,Y)+o({\rho }_5^4)\\ G(X,Y)+o({\rho }_5^4)\end{array}\right),$(3.9)

where ${\rho }_5=\sqrt{X^2+Y^2}$,

$F(X,Y)=c_{20}{u}^2+c_{11}uv+c_{02}{v}^2+c_{30}{u}^3+c_{21}{u}^{2}v+c_{12}u{v}^2+c_{03}{v}^3,$

$G(X,Y)=d_{20}{u}^2+d_{11}uv+d_{02}{v}^2+d_{30}{u}^3+d_{21}{u}^{2}v+d_{12}u{v}^2+d_{03}{v}^3,$

of which $u=a_{01}Y$, and $v=\mu X+(1-\omega )Y$,

$c_{20}=\frac{a_{20}(\omega -1)}{\mu a_{01}}+\frac{b_{20}}{\mu },c_{11}=\frac{a_{11}(\omega -1)}{\mu a_{01}}+\frac{b_{11}}{\mu },c_{02}=\frac{a_{02}(\omega -1)}{\mu a_{01}}+\frac{b_{02}}{\mu },$

$c_{30}=\frac{a_{30}(\omega -1)}{\mu a_{01}}+\frac{b_{30}}{\mu },c_{21}=\frac{a_{21}(\omega -1)}{\mu a_{01}}+\frac{b_{21}}{\mu },c_{12}=\frac{a_{12}(\omega -1)}{\mu a_{01}}+\frac{b_{12}}{\mu },$

$c_{03}=\frac{a_{03}(\omega -1)}{\mu a_{01}}+\frac{b_{03}}{\mu },d_{20}=\frac{a_{20}}{a_{01}},d_{11}=\frac{a_{11}}{a_{01}},d_{02}=\frac{a_{02}}{a_{01}},d_{30}=\frac{a_{30}}{a_{01}},$

$d_{21}=\frac{a_{21}}{a_{01}},d_{12}=\frac{a_{12}}{a_{01}},d_{03}=\frac{a_{03}}{a_{01}}.$

Moreover,

$F_{XX}{|}_{(0,0)}=2c_{02}{\mu }^2,F_{XY}{|}_{(0,0)}=c_{11}{a}_{01}\mu +2c_{02}\mu (1-\omega ),$

$F_{YY}{|}_{(0,0)}=2c_{20}{a}_{01}^2+2c_{11}{a}_{01}(1-\omega ),F_{XXX}{|}_{(0,0)}=6c_{03}{\mu }^3,$

$F_{XXY}{|}_{(0,0)}=2c_{21}{a}_{01}{\mu }^2+6c_{03}{\mu }^2(1-\omega ),$

$F_{XYY}{|}_{(0,0)}=2c_{21}{a}_{01}^2\mu +4c_{12}{a}_{01}\mu (1-\omega )+6c_{03}\mu (1-\omega {)}^2,$

$F_{YYY}{|}_{(0,0)}=4(1-\omega {)}^3+6c_{30}{a}_{01}^3+4c_{21}{a}_{01}^2(1-\omega )+6c_{12}{a}_{01}(1-\omega {)}^2,$

$G_{XX}{|}_{(0,0)}=2d_{02}{\mu }^2,G_{XY}{|}_{(0,0)}=d_{11}{a}_{01}\mu +2d_{02}\mu (1-\omega ),$

$G_{YY}{|}_{(0,0)}=2d_{20}{a}_{01}^2+2d_{11}{a}_{01}(1-\omega ),G_{XXX}{|}_{(0,0)}=6c_{03}{\mu }^3,$

$G_{XXY}{|}_{(0,0)}=2d_{21}{a}_{01}{\mu }^2+6d_{03}{\mu }^2(1-\omega ),$

$G_{XYY}{|}_{(0,0)}=2d_{21}{a}_{01}^2\mu +4d_{12}{a}_{01}\mu (1-\omega )+6d_{03}\mu (1-\omega {)}^2,$

$G_{YYY}{|}_{(0,0)}=4(1-\omega {)}^3+6d_{30}{a}_{01}^3+4d_{21}{a}_{01}^2(1-\omega )+6d_{12}{a}_{01}(1-\omega {)}^2.$

The fifth step. In order to make sure that the system (3.8) undergoes a Neimark-Sacker bifurcation and to determine the stability and direction of the bifurcation curve, the discriminating quantity $L$ should be calculated and required not to be zero, where

(3.10) $L=-Re(\frac{(1-2{\lambda }_1){\lambda }_2^2}{1-{\lambda }_1}{\zeta }_{20}{\zeta }_{11})-\frac{1}{2}|{\zeta }_{11}{|}^2-|{\zeta }_{02}{|}^2+Re({\lambda }_2{\zeta }_{21}),$(3.10)

${\zeta }_{20}=\frac{1}{8}[F_{XX}-F_{YY}+2G_{XY}+i(G_{XX}-G_{YY}-2F_{XY})]{|}_{(0,0)},$

${\zeta }_{11}=\frac{1}{4}[F_{XX}+F_{YY}+i(G_{XX}+G_{YY})]{|}_{(0,0)},$

${\zeta }_{02}=\frac{1}{8}[F_{XX}-F_{YY}-2G_{XY}+i(G_{XX}-G_{YY}+2F_{XY})]{|}_{(0,0)},$

${\zeta }_{21}=\frac{1}{16}[F_{XXX}+F_{XYY}+G_{XXY}+G_{YYY}$

$+i(G_{XXX}+G_{XYY}-F_{XXY}-F_{YYY})]{|}_{(0,0)}.$

According to [29], one can see that if $L<(>)0$, then an attracting invariant closed curve bifurcates from the fixed point $E_2$ for $b>(<)b_1$. Based on the above analysis, the following result can be summarized.

Theorem 3.2.

Suppose the parameters $a,b,d$ in the space $S_2$. When the parameter $b$ varies in a small neighbourhood of the critical value $b_1=d\sqrt{\frac{d+3}{d+1}}$, a Neimark-Sacker bifurcation occurs at the fixed point $E_2$ of the system (1.7) if $L\ne 0$. Furthermore, if $L<(>)0$, then an attracting (repelling) invariant closed curve bifurcates from the fixed point $E_2$ for $b>(<)b_1$.

4. Numerical simulation

In this section, by using Matlab software, we can obtain the bifurcation diagrams, Lyapunov exponents and phase portraits of the system (1.7) at the fixed points $E_1$ and $E_2$ with some peculiar parameter values to corroborate the theoretical results previously derived.

First consider the fixed point $E_1$. Choose the parameter $a=1$, $d=2$, the initial values $(x_0,y_0)=(0.996,0.994)$ and $b\in (0,4)$, then we can obtain Figure 1(a) and observe the existence of transcritical bifurcation when $b\approx 2$, which is in accordance with the result in Theorem 3.1. Figure 1(b) means the spectrum of maximum Lyapunov exponent with respect to the parameter $b$.

Figure 1. Bifurcation of the system (7) in $(b,x)$-plane and maximal lyapunov exponents.

Fix the parameter values $a=0.05$, $d=0.55$, let $b\in (0.7,0.93)$ and take the initial values $(0.44,4.5)$, $(0.44,7.47)$ in Figures 2 and 3, respectively. Figure 1(c) illustrates the bifurcation diagram of $(b,x)$-plane, and shows the existence of a Neimark-Sacker bifurcation at the fixed point $E_2=(0.45486,7.35321)$ when $b=b_1=0.8155$, whose multipliers are ${\lambda }_{1,2}=0.908855\pm 0.376306$ with $|{\lambda }_{1,2}|=1$. We can also clearly see that the fixed point $E_2$ is unstable when $b>b_1$, stable when $b<b_1$ and that a Neimark-Sacker bifurcation occurs when $b=b_1$.

Figure 2. Phase portraits for the system (7) with $\mathrm{a}=0.05$, $\mathrm{d}=0.55$ and different $b$ with the initial value $(x_0,y_0)=(0.44,4.5)$ outside the closed orbit.

Figure 3. Phase portraits for the system (7) with $\mathrm{a}=0.05$, $\mathrm{d}=0.55$ and different $b$ with the initial value $(x_0,y_0)=(0.44,7.47)$ inside the closed orbit.

Figure 1(d) depicts the corresponding maximal Lyapunov exponents, which are subtractive for the parameter $b\in \left(0.7,0.93\right)$ which means the occurrence of chaos [24] in the system (7). Figures 2(a-i) and 3(a-d) show that the dynamical properties of the fixed point $E_2$ change from stable to unstable as the value of the parameter $b$ increases and there is an occurrence of invariant closed curve around $E_2$ when $b=b_1$, which fits the conclusion of Theorem 3.2 very well. Not only that, we also give the dynamic behaviour of the system (7) when $d=1.5$, as shown in Figure 4. Figure 4 describes such a fact that the predator death rate parameter $s$ is too big, bigger than its consumption rate $r$ of the prey, which is equivalent to $d>1$.

Figure 4. Phase portraits for the system (7) with $\mathrm{a}=0.05$, $\mathrm{d}=1.5$ and different $b$ with the initial value $(x_0,y_0)=(0.46,6.5)$.

5. Conclusion

In this paper, we study the dynamical behaviour of a discrete predator-prey system with herd behaviour and square root functional response by using the semi-discretization method. Given the parameter conditions, we completely formulate the existence and stability of two nonnegative equilibria $E_1(1,0)$, and $E_2(\frac{d^2}{b^2},\frac{1}{a}(1-\frac{d^2}{b^2})\frac{d}{b})$ for $b>d$. We also derive the sufficient conditions for the transcritical bifurcation and Neimark-Sacker bifurcation to occur at the fixed points $E_1$ and $E_2$ respectively in certain parameter spaces. The most key is that our model has more simple form, namely, with the fewer number of parameters, and more clear expressions of bifurcation results. Unlike the known work that only presents the existence results for local bifurcation of this system, we also state the stability and direction of the closed orbits bifurcated. To our surprise, the system (3) cannot undergo the Turing-Hopf bifurcation, whereas the known result has, which is due to the different discrete method we adopt. Our results elaborate this fact: For the same one continuous system, different discrete ways probably deduce different results. Finally, some interesting dynamical properties for Neimark-Sacker bifurcation are obtained through applying numerical simulations.

By comparing these figures, one finds that, the dynamical behaviour of this system is entirely reasonable from an ecological point of view:

1. When $d>1$, the predator naturally dies off, leaving the prey to maintain itself over time at the carrying capacity.

2. When $d<1$, the predator and prey coexist in a stable equilibrium.

The numerical simulation results confirm that the prey in the system exhibits strong herd structure when the square root is used as a functional response in the system, implying that the predator generally interacts with the prey along the outer corridor of the herd. This behaviour makes the predator hard to obtain food (such as tigers on a vast grassland). If the situation is not resolved, the predator will be in peril of extinction. In order to prevent the predator from dying out, one must take effective measures to ensure that such an outcome does not occur.

Authors’ contributions

All authors contributed equally and significantly in writing this paper. All authors read and approved the final manuscript.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is partly supported by the National Natural Science Foundation of China (grant: 61473340), the Distinguished Professor Foundation of Qianjiang Scholar in Zhejiang Province (grant: F703108L02), and the Natural Science Foundation of Zhejiang University of Science and Technology (grant: F701108G14).

There are no applicable data associated with this manuscript.

Disclosure statement

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

Additional information

Funding

The work was supported by the National Natural Science Foundation of China [61473340]; Distinguished Professor Foundation of Qianjiang Scholar in Zhejiang Province [F703108L02].

appendix

For readers’ convenience, we have included a Definition and a Lemma in the appendix.

Definition 5.1.

Let $E(x,y)$ be a fixed point of the system (1.7) with multipliers ${\lambda }_1$ and ${\lambda }_2$.

$(i)$ If $|{\lambda }_1|<1$ and $|{\lambda }_2|<1$, $E(x,y)$ is called sink, so a sink is locally

asymptotically stable.

$(ii)$ If $|{\lambda }_1|>1$ and $|{\lambda }_2|>1$, $E(x,y)$ is called source, so a source is

locally asymptotically unstable.

$(iii)$ If $|{\lambda }_1|<1$ and $|{\lambda }_2|>1$ $($or $|{\lambda }_1|>1$ and $|{\lambda }_2|<1$$)$, $E(x,y)$ is called

saddle.

$(iv)$ If either $|{\lambda }_1|=1$ or $|{\lambda }_2|=1$, $E(x,y)$ is called to be non-hyperbolic.

Lemma 5.2.

Let $F(\lambda )={\lambda }^2+B\lambda +C$, where $B$ and $C$ are two real constants. Suppose ${\lambda }_1$ and ${\lambda }_2$ are two roots of $F(\lambda )=0$. Then, the following statements hold.

$(i)$ If $F(1)>0$, then

$(i.1)$ $|{\lambda }_1|<1$ and $|{\lambda }_2|<1$ if and only if $F(-1)>0$ and $C<1$;

$(i.2)$ ${\lambda }_1=-1$ and ${\lambda }_2\ne -1$ if and only if $F(-1)=0$ and $B\ne 2$;

$(i.3)$ $|{\lambda }_1|<1$ and $|{\lambda }_2|>1$ if and only if $F(-1)<0$;

$(i.4)$ $|{\lambda }_1|>1$ and $|{\lambda }_2|>1$ if and only if $F(-1)>0$ and $C>1$;

$(i.5)$ ${\lambda }_1$ and ${\lambda }_2$ are a pair of conjugate complex roots with $|{\lambda }_1|=|{\lambda }_2|=1$

if and only if $-2<B<2$ and $C=1$;

$(i.6)$ ${\lambda }_1={\lambda }_2=-1$ if and only if $F(-1)=0$ and $B=2$.

$(ii)$ If $F(1)=0$, namely, $1$ is one root of $F(\lambda )=0$, then the other root $\lambda$

satisfies $|\lambda |=(<,>)1$ if and only if $|C|=(<,>)1$.

$(iii)$ If $F(1)<0$, then $F(\lambda )=0$ has one root lying in $(1,\mathrm{\infty })$. Moreover,

$(iii.1)$ the other root $\lambda$ satisfies $\lambda <(=)-1$ if and only if $F(-1)<(=)0$;

$(iii.2)$ the other root $-1<\lambda <1$ if and only if $F(-1)>0$.