Bifurcation analysis in a discrete predator–prey model with herd behaviour and group defense

1.   Introduction and preliminaries

Over the past several decades, the predator-prey interaction has become a hot point of studies in biomathematics ^{[1,2,3,4,5,6,7,8,9,10]}. Because differential equations can assume that generations overlap and that populations vary continously in time, the general model for predator-prey interaction may be written as

where $x$ and $y$ are expressed as prey and predator population sizes (or densities), respectively, $f(x)$ denotes the growth rate of prey with the absence of predator, $g(x, y)$ represents the amount of prey consumed per predator per unit time, $h(x, y)$ is on behalf of per capita predator production, and $m$ is the intrinsic death rate of predator. See also ^[1].

Due to the realistic meaning of $f(x)$, one can assume that the prey grows logistically with growth rate $r$ and carrying capacity $k$ in the absence of predator (i.e., $f(x) = r(1-\frac{x}{k})$). Hence the system (1.1) can be written as

where $e$ is the conversion effciency.

As for the functional response $g(x, y)$, there are many different kinds of forms. Bian et al. proposed a system with the Beddington-DeAngelis funcional response ^[5]; De Assis et al. proposed a system with the square-root functional response ^[7] and so on. Notice the fact that in the natural ecosystem, many species may gather together and form herds to either search for food resources or to defend the predators, which means that all members of a group do not interact at one time. This behaviour is often called herd behaviour. In this paper, one talks about the following system ^[6,7]:

Here, the funcional response $\frac{ax}{\sqrt{x+\tilde{h}}}$ can be expressed as the function of the ratio of prey to predator, where $\tilde{h}$ is a threshold for the transition between herd grouping and solitary behaviour and $a$ is the maximum value of prey consumed by each predator per unit time. In this system, all parameters are positive. The biological meanings for the parameters $r$, $k$, $e$, and $m$ are the same as in (1.2).

For the sake of simplicity of mathematical analysis, let $\frac{x}{k} \to x, mt \to t, \frac{{y}}{ek} \to y, \frac{r}{m} \to \gamma, \frac{ae\sqrt{k}}{m} \to \beta, \frac{\tilde{h}}{k} \to h$, then one can derive an equivalent form of the system (1.3) as follows:

This continuous system has been discussed in ^[6,7], but its discrete version has not been investgated as of yet. To be honest, it is very difficult to solve a complicate continuous equation or system without using computer. Therefore, one natuarally wishes to consider the corresponding discrete version of a continuous model. One tries to use various methods to derive the discrete model of the system (1.4) to make it easily studied ^{[8,9,10,11,12,13,14,15,16]}. In this paper, we adopt a semi-discretazation method, which does not need to consider the step size, 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

It is easy to see that the system (1.5) has piecewise constant arguments, and that a solution $(x(t), y(t))$ of the system (1.5) for $t\in [0, +\infty)$ has the following characteristics:

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

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

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

where $x_n = x(n)$ and $y_n = y(n)$. Letting $t\to (n+1)^-$ in the system (1.6) produces

where the parameters $h, \beta, \gamma > 0$, and their biological meanings are the same as in (1.4). The system (1.7) will be considered in the sequel.

The rest of the paper is organized as follows. In Section 2, we investigate the existence and stability of the fixed points of the system (1.7) in detail. In Section 3, we derive the sufficient conditions for transcritical bifurcation, period-doubling bifurcation, and Neimark-Sacker bifurcation of the system (1.7) to occur. In Section 4, numerical simulations are performed to illustrate the above theoretical results. In the end, some brief conclusions are stated in Section 5.

2.   Existence and stability of fixed points

Considering the biological meaning of the system (1.7), we discuss the existence and stability of non-negative fixed points of the system (1.7) in this section. By solving the equations of fixed points of system (1.7)

it's easy to find that there are three nonnegative fixed points $E_{0} = (0, 0)$, $E_{1} = (1, 0)$, and $E_{2} = (x_{0}, y_{0})$ for $\beta > \sqrt{ h+1}$, where

The Jacobian matrix of the system (1.7) at a fixed point $E(x, y)$ is

whose charactertistic polynomial reads as

where

In order to analyze the properties of the fixed points of the system (1.7), we utilize the Appendix definition and Lemma ^[17,18,19].

By using Definition 5.1 and Lemma 5.2 in the Appendix, the following conclusions can be obtained.

Theorem 2.1. The fixed point $E_0 = (0, 0)$ of the system (1.7) is a saddle.

The proof for this theorem is simple and omitted here.

Theorem 2.2. The type of the fixed point $E_1 = (1, 0)$ of the system (1.7) complies with the following results:

Proof. The Jacobian matrix $J({E_1})$ of the system (1.7) at the fixed point $E_1$ reads

Obviously, $\lambda_1 = 1-\gamma$ and $\lambda_2 = e^{\frac{\beta}{\sqrt{h+1}}-1}$.

When $0 < \gamma < 2$, $|\lambda_1| < 1$. If $0 < \beta < \sqrt{h+1}$, then $|\lambda_2| < 1$, so $E_1$is a sink; if $\beta = \sqrt{h+1}$, then $|\lambda_2| = 1$, therefore $E_1$ is non-hyperbolic; if $\beta > \sqrt{h+1}$, meaning $|\lambda_2| > 1$, then $E_1$ is a saddle.

When $\gamma = 2$, which reads $|\lambda_1| = 1$, $E_1$ is non-hyperbolic.

When $\gamma > 2$, $|\lambda_1| > 1$. If $0 < \beta < \sqrt{h+1}$, then $|\lambda_2| < 1$, so $E_1$ is a saddle; if $\beta = \sqrt{h+1}$, then $|\lambda_2| = 1$, therefore $E_1$ is non-hyperbolic; if $\beta > \sqrt{h+1}$, implying $|\lambda_2| > 1$, then $E_1$ is a source. The proof is complete.

We can easily derive the following result.

Lemma 2.3. Consider the function $f(x) = 4x^2-4x+7+(2x-7)\sqrt{4x ^{2}+20x+1}$ with $x \in (1, \infty)$. Then $f(x)$ is strictly increasing for $x \in (1, \infty)$, Furthermore, $f(x)$ has a unique positive root $X_0$ in (2, 2.5).

Proof. Evidently, $f'(x) = 4(2x-1)+\frac{16x^2+32x-68}{\sqrt{4x ^{2}+20x+1}}$ and $f''(x) = 8+\frac{64x^3+520x^2+912x712}{(\sqrt{4x ^{2}+20x+1})^3} > 0$, so, for $x > 1$, $f'(x) > f'(1) = 0$. Hence, $f(x)$ is strictly increasing for $x \in (1, \infty)$. Again, $f(2) = 15-3\sqrt{57} < 0$ and $f(2.5) = 22-2\sqrt{76} > 0$. Therefore, $f(x)$ has a unique positive root $X_0$ in (2, 2.5).

Now consider the stability of the fixed point $E_2$.

Theorem 2.4. For $\beta > \sqrt{h+1}$, $E_2 = (x_0, y_0) = (\frac{1+\sqrt{1+4h\beta^2}}{2\beta^2}, \gamma \frac{1+\sqrt{1+4h\beta^2}}{2\beta^2} (1-\frac{1+\sqrt{1+4h\beta^2}}{2\beta^2}))$ is a positive fixed point of the system (1.7).

Let $X_0$ be the unique positive root of the function $f(x) = 4x^2-4x+7+(2x-7)\sqrt{4x ^{2}+20x+1}$ in (2, 2.5). Put $\beta_{0} = \sqrt{4h+2}$ and $h_0 = \frac{4\beta ^{4}-4\beta ^{2}+7+(2\beta ^{2}-7)\sqrt{4\beta ^{4}+20\beta ^{2}+1}}{72\beta ^{2}}$. Denote $\gamma_{0} = \frac{8 \beta ^2(1+\sqrt{1+4h \beta ^{2}})}{3(1+4h \beta^2)+(7-2\beta^2) \sqrt{1+4h \beta ^{2} }+4(1-\beta ^2)}$, where $h > h_0$. Then the following consequences hold about the fixed point $E_2$ illustrated in the Table 2.

Proof. The Jacobian matrix $J({E_2})$ of the system (1.7) at the fixed point $E_2$ is

whose characteristic polynomial can be written as

where

Note that $x_{0} = \frac{1+\sqrt{1+4h \beta ^{2} }}{2 \beta ^{2} }, y_{0} = \gamma x_{0} (1-x_{0})$, and $E_2 = (x_0, y_0)$ is a positive fixed point, so $0 < x_{0} < 1$. It's easy to calculate that

If $6\beta ^{2}{x^*}^{2}-(2\beta ^{2}-1)x^{*}-1 = 0$ and $x^{*} > 0$, then $x^{*} = \frac{2\beta ^{2}-1+\sqrt{4\beta ^{4}+20\beta ^{2}+1}}{12\beta ^{2}}$. Simultaneously, it is easy to prove $x^{*} < \frac{1}{2}$.

Notice that $0 < h < \beta^2-1$. Moreover, $x_0 > ( = , < )x^* \Leftrightarrow h > ( = , < ) h_0.$ Additionally, $\beta ^{2}-1-h_0 = \frac{68\beta ^{4}-68\beta ^{2}-7-(2\beta ^{2}-7)\sqrt{4\beta ^{4}+20\beta ^{2}+1}}{72\beta ^{2}}.$ Set $x = \beta ^{2} > 1$. Denote

and

Lemma 2.3 tells us that $f(x)$ is strictly increasing for $x > 1$ and has a unique positive root $X_0$ in (2, 2.5). From this one can see

Obiviously, $g(1) = 18 > 0$, and $g'(x) = \frac{68(2x-1)\sqrt{4x ^{2}+20x+1}-16x^2-32x+68}{\sqrt{4x ^{2}+20x+1}} > \frac{68(2x-1)(2x+1)-16x^2-32x+68}{\sqrt{4x ^{2}+20x+1}} = \frac{256x^2-32x}{\sqrt{4x ^{2}+20x+1}} > 0$. So, $g(x) > g(1) > 0$ for $x > 1$. This implies that $h_0 < \beta^2-1$ always holds.

It is easy to see $x_0 > ( = , < )x^{*} \Leftrightarrow 6\beta ^{2}x_{0}^{2}-(2\beta ^{2}-1)x_{0}-1 > ( = , < )0$. From $F(-1) = 4-\frac{\gamma(6\beta ^{2}x_{0}^{2}+2\beta ^{2}x_{0}+x_{0}-1)}{2\beta ^{2}x_{0}} = 0$, one has

Again, $\beta > ( = , < )\beta _{0} = \sqrt{4h+2} \Leftrightarrow x_{0} < ( = , > )\frac{1+\sqrt{1+4h \beta_{0} ^{2} }}{2 \beta_{0} ^{2} } = \frac{1}{2}$.

Now, one considers the following two cases:

1) Case Ⅰ: $0 < h\leqslant \frac{X_0-2}{4}$. Then $\beta_{0} = \sqrt{4h+2}\leq \sqrt{X_0}$.

(a) Subcase 1: $\sqrt{h+1} < \beta \leq \sqrt{X_0}$. Then $h_0 \leq 0 < h$, implying $x^* < x_0$ and $\gamma_0 > 0$.

ⅰ. If $0 < \gamma < \gamma_{0}$, then $F(-1) > 0$.

● For $\beta < \beta _{0}$, $q < 1$, which reads $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| < 1$ by Lemma 6.2$(i.1)$. So, $E_2$ is a sink;

● For $\beta = {\beta _0}$, $q = 1$. Lemma 6.2$(i.5)$ shows that $\left| {{\lambda _1}} \right| = \left| {{\lambda _2}} \right| = 1$, so $E_2$ is non-hyperbolic;

● For $\beta > {\beta _0}$, $q > 1$. Lemma 6.2$(i.4)$ shows that $\left| {{\lambda _1}} \right| > 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a source.

ⅱ. If $\gamma = \gamma_{0}$, then $F(-1) = 0$. In other words, $-1$ is a root of the characteristic polynomial, namely $E_2$ is non-hyperbolic.

ⅲ. If $\gamma > \gamma_{0}$, then $F(-1) < 0$. Using Lemma 6.2$(i.3)$, we conclude that $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a saddle.

(b) Subcase 2: $\beta > \sqrt{X_0}$. Then $h_0 > 0$.

ⅰ. If $0 < h\leq h_0$, then $0 < x_0 \leq x^*$, implying that $6\beta ^{2}x_{0}^{2}-(2\beta ^{2}-1)x_{0}-1 \le 0$. So, $F(-1)\geq 4 > 0$. From $\beta > \sqrt{X_0} \geq \beta _0$, we see $q > 1$. Lemma 6.2$(i.4)$ shows that $\left| {{\lambda _1}} \right| > 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a source.

ⅱ. If $h > h_0$, then $x^* < x_0 < 1$, implying that $\gamma_0 > 0$.

A. If $0 < \gamma < \gamma_{0}$, then $F(-1) > 0$. For $\beta > \sqrt{X_0} \geq \beta _{0}$, $q > 1$. Lemma 6.2$(i.4)$ shows that $\left| {{\lambda _1}} \right| > 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a source.

B. If $\gamma = \gamma_{0}$, then $F(-1) = 0$. In other words, $-1$ is one root of the characteristic polynomial, namely, $E_2$ is non-hyperbolic.

C. If $\gamma > \gamma_{0}$, then $F(-1) < 0$. Lemma 6.2$(i.3)$ shows that $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a saddle.

2) Case Ⅱ: $h > \frac{X_0-2}{4}$. Then $\beta_{0} = \sqrt{4h+2} > \sqrt{X_0}$.

(a) Subcase 1: $\sqrt{h+1} < \beta \leq \sqrt{X_0}$. Then $h_0 \leq 0 < h$, so, $x^* < x_0$ and hence $r_0 > 0$.

ⅰ. If $0 < \gamma < \gamma_{0}$, then $F(-1) > 0$. For $\sqrt{h+1} < \beta \leq \sqrt{X_0} < \beta _{0}$, $q < 1$, which reads $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| < 1$ by Lemma 6.2$(i.1)$. Therefore, $E_2$ is a sink.

ⅱ. If $\gamma = \gamma_{0}$, then $F(-1) = 0$. Hence, $E_2$ is non-hyperbolic.

ⅲ. If $\gamma > \gamma_{0}$, then $F(-1) < 0$. Lemma 6.2$(i.3)$ shows that $E_2$ is a saddle.

(b) Subcase 2: $\beta > \sqrt{X_0}$. Then $h_0 > 0$.

ⅰ. If $\frac{X_0-2}{4} < h\leq h_0$, then $x_0 \leq x^*$, so, $F(-1)\geq 4 > 0$.

● For $\sqrt{X_0} < \beta < \beta _{0}$, $q < 1$, which reads $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| < 1$ by Lemma 6.2$(i.1)$, thus, $E_2$ is a sink;

● For $\beta = {\beta _0}$, $q = 1$. Lemma 6.2$(i.5)$ shows that $\left| {{\lambda _1}} \right| = \left| {{\lambda _2}} \right| = 1$, so $E_2$ is non-hyperbolic;

● For $\beta > {\beta _0}$, $q > 1$. It follows from Lemma 6.2$(i.4)$ that $\left| {{\lambda _1}} \right| > 1$ and $\left| {{\lambda _2}} \right| > 1$, hence $E_2$ is a source.

ⅱ. If $h > h_0$, then $x^* < x_0$, so, $\gamma_0 > 0$.

A. If $0 < \gamma < \gamma_{0}$, then $F(-1) > 0$.

● For $\sqrt{X_0} < \beta < \beta _{0}$, $q < 1$. Lemma 6.2$(i.1)$ tells us $E_2$ is a sink;

● For $\beta = {\beta _0}$, $q = 1$. Therefore, $E_2$ is non-hyperbolic;

● For $\beta > {\beta _0}$, $q > 1$. Lemma 6.2$(i.4)$ shows that $E_2$ is a source.

B. If $\gamma = \gamma_{0}$, then $F(-1) = 0$, which shows $E_2$ is non-hyperbolic.

C. If $\gamma > \gamma_{0}$, then $F(-1) < 0$. Using Lemma 6.2$(i.3)$, we conclude that $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a saddle.

Summarizing the above analysis, the proof is complete.

3.   Bifurcation analysis

In this section, we apply the center manifold theorem and bifurcation theory to investigate the local bifurcation problems of the system at the fixed points $E_1$ and $E_2$.

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

It follows from Eq (1.4) that the fixed point $E_1$ always exists, regardless of what values the parameters $\beta$ and $\gamma$ take. One can see from Theorem 2.2 that the fixed point $E_1$ is a non-hyperbolic fixed point when $\beta = \sqrt{h+1}$ or $\gamma = 2$. As soon as the parameters $\beta$ or $\gamma$ goes through corresponding critical values, the dimensional numbers for the stable manifold and the unstable manifold of the fixed point $E_1$ vary. Therefore, a bifurcation probably occurs. Now, the considered parameter case is divided into the following three subcases:

Case Ⅰ: $\beta = \sqrt{h+1}$, $\gamma\neq 2$;

Case Ⅱ: $\beta\neq \sqrt{h+1}$, $\gamma = 2$;

Case Ⅲ: $\beta = \sqrt{h+1}$, $\gamma = 2$.

First we consider Case Ⅰ: $\beta = \sqrt{h+1}$, $\gamma\neq 2$, i.e., the parameters $\left (h, \beta, \gamma \right)\in \Omega_1 = \left \{\left (h, \beta, \gamma \right)\in R_{+}^{3}\mid h > 0, \beta > 0, \gamma > 0, \gamma \neq 2 \right \}$. Then, the following result is obtained.

Theorem 3.1. Suppose the paramenters $(h, \beta, \gamma) \in \Omega_1$. Let $\beta_1 = \sqrt{h+1}$. If the parameter $\beta$ varies in a small neighborhood of the critical value ${\beta _1}$, then the system (1.7) experiences a transcritical bifurcation at the fixed point $E_1$ when the parameter ${\beta}$ goes through the critical value ${\beta _1}$.

Proof. First, assume that $u_n = x_n-1, v_n = y_n-0$, which transforms the fixed point $E_1$ to the origin, and the system (1.7) to

Second, giving a small perturbation $\beta^*$ of the parameter $\beta$ around $\beta_1$, i.e., $\beta^* = \beta-\beta_1$ with $0 < | \beta^* |\ll 1$, and letting $\beta_{n+1}^* = \beta_n^* = \beta^*$, the system (3.1) is perturbed into

By the Taylor expansion, the system (3.2) at $(u_n, v_n, \beta_n^*) = (0, 0, 0)$ can be written as

where $\rho_1 = \sqrt {u_n^2+v_n^2+{\beta_n^*}^2},$

It is easy to derive the three eigenvalues of the matrix

to be $\lambda _1 = 1-\gamma$ and $\lambda _2 = \lambda _3 = 1$ with corresponding eigenvectors

Notice $0 < \gamma \neq 2$ implies that $|\lambda_1|\neq 1$.

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

then,

Taking the following transformation

the system (3.3) is changed into the following form

where $\rho_2 = \sqrt{X_{n}^{2}+Y_{n}^{2}+\delta _n^{2}}$,

Assume on the center manifold

where $\rho_3 = \sqrt{Y_{n}^{2}+\delta _n^{2}}$, then, from

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

By comparing the corresponding coefficients of terms with the same order in the above center manifold equation, it is easy to derive that

Therefore, the system (3.4) restricted to the center manifold is given by

It is not difficult to calculate

According to (21.1.43)–(21.1.46) in [24, p507], for a transcritical bifurication to occur, all conditions hold, hence, the system (1.7) undergoes a transcritical bifurcation at the fixed point $E_1$. The proof is over.

Next we consider Case Ⅱ: $\beta\neq \sqrt{h+1}$, $\gamma = 2$. By Theorem 2.2, one can see that $\lambda_1 = -1$ and $|\lambda_2|\neq 1$when $\beta\neq \sqrt{h+1}$, $\gamma = 2$. Thereout, the following result can be derived.

Theorem 3.2. Let $\gamma_{1} = 2.$ Suppose the paramenters $\left (h, \beta, \gamma \right)\in \Omega_2 = \left \{\left (h, \beta, \gamma \right)\in R_{+}^{3}\mid h > 0, \beta > 0, \beta \neq \sqrt{h+1}, \gamma > 0 \right \}$. If the parameter $\gamma$ varies in a small neighborhood of the critical value ${\gamma _1}$, then the system (1.7) undergoes a period-doubling bifurcation at the fixed point $E_1$ when the parameter ${\gamma}$ goes through the critical value ${\gamma _1}$.

Proof. Shifting $E_1 = (1, 0)$ to the origin $O(0, 0)$ and giving a small perturbation ${\gamma^*}$ of the parameter ${\gamma}$ at the critical value ${\gamma_1}$ with ${0 < |\gamma^*|\ll 1}$, the system (3.1) is changed into the following form:

Set $\gamma_{n+1}^* = \gamma_n^* = \gamma^*$, then (3.5) can be seen as

By the Taylor expansion, the system (3.6) at $(u_n, v_n, \beta_n^*) = (0, 0, 0)$ can be expended into

where $\rho_4 = \sqrt {u_n^2+v_n^2+{\gamma_n^*}^2},$

It is not difficult to derive the three eigenvalues of the matrix

to be $\lambda _1 = -1$, $\lambda _2 = {\mathrm{e}}^{(\frac{\beta}{\sqrt{h+1}}-1)}$ and $\lambda _3 = 1$ with corresponding eigenvectors

Notice $\beta \neq \sqrt{h+1}$ implies $|\lambda_2|\neq 1$.

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

then,

Taking the following transformation

the system (3.7) is changed into the following form

where $\rho_5 = \sqrt{X_{n}^{2}+Y_{n}^{2}+\delta _n^{2}}$,

Suppose on this center manifold

where $\rho_6 = \sqrt{X_{n}^{2}+\delta _n^{2}}$, which must satisfy

Similar to Case Ⅰ, one can establish the corresponding center manifold equation. Comparing the corresponding coefficients of terms with the same type in the equation produces

That is to say, $Y_{n} = h(X_{n}, \delta_{n}) = o(\rho_6^2).$ Therefore, the center manifold equation is given by

Thereout, one has

Therefore, the following results are derived:

which, according to (21.2.17)–(21.2.22) in [24, p516], satisfy all conditions for a period-doubling bifurcation to occur. Therefore, the system (1.7) undergoes a period-doubling bifurcation at $E_1$. Again,

Therefore, the period-two orbit bifurcated from $E_1$ lies on the right of $\gamma_{1} = 2$.

Of course, one can also compute the following two quantities, which are the transversal condition and non-degenerate condition for judging the occurrence and stability of a period-doubling bifurcation, respectively (see ^{[3,15,16,17,18]}),

It is easy to say $\alpha_1 = -1$ and $\alpha_2 = \frac{2}{3}$. Due to $\alpha_{2} > 0$, the period-two orbit bifurcated from $E_1$ is stable. The proof is complete.

Finally, we consider Case Ⅲ: $\beta = \sqrt{h+1}$, $\gamma = 2$. At this time, the two eigenvalues of the linearized matrix evaluated at this fixed point $E_1$ are $\lambda_1 = -1$ and $\lambda_2 = 1$. The bifurcation problem in this case is very complicated and will be considered future work.

3.2.   For fixed point $E_2= (\frac{1+\sqrt{1+4h\beta^2}}{2\beta^2}, \gamma \frac{1+\sqrt{1+4h\beta^2}}{2\beta^2} (1-\frac{1+\sqrt{1+4h\beta^2}}{2\beta^2}))$

Consider the bifurcation of the system (1.7) at the fixed point $E_2$. The parameters are divided into the following three cases:

Case Ⅰ: $\beta = \sqrt{4h+2}$, $\gamma \neq \frac{8 (4h+2)}{4h +1}$;

Case Ⅱ: $\beta \neq \sqrt{4h+2}$, $\gamma = \frac{8 \beta ^2(1+\sqrt{1+4h \beta ^{2}})}{3(\sqrt{1+4h \beta ^{2}})^2+(7-2\beta^2) \sqrt{1+4h \beta ^{2} }+4(1-\beta ^2)}$;

Case Ⅲ: $\beta = \sqrt{4h+2}$, $\gamma = \frac{8 (4h+2)}{4h +1}$.

According to our calculations, there is no bifurcation under Case Ⅱ. Additionally, the bifurcation problem in case Ⅲ is very complicated and will be considered future work. Therefore, we only consider Case Ⅰ.

Suppose the paramenters

Then the following result may be obtained.

Theorem 3.3. Suppose the paramenters $(h, \beta, \delta) \in \Omega_3$ and meet $\gamma < \frac{8 (4h+2)}{4h +1}$. Let $\beta_2 = \sqrt{4h+2}$. Then the system (1.7) undergoes a Neimark-Sacker bifurcation at the fixed point $E_2$ when the parament $\beta$ varies in a small neighborhood of the critical value ${\beta _2}$. Moreover, if $L < (>)0$ in (3.13), then a (an) stable (unstable) invariant closed orbit is bifurcated out from the fixed point $E_2$ of system (1.7) when $\beta > (<) \beta_2$.

Proof. First, give a small perturbation $\beta^{**}$ of the parameter $\beta$ around $\beta_2$ in the system (3.1), i.e., $\beta^{**} = \beta-\beta_2$ with $0 < |\beta^{**}|\ll1$, and set $x_{01} = x_{01}(\beta^{**}) = \frac{1+\sqrt{1+4h(\beta^{**}+\beta_2)^2}}{2(\beta^{**}+\beta_2)^2}$ and $y_{01} = \gamma x_{01} (1-x_{01})$. Under the perturbation, the system (3.1) reads

The characteristic equation of the linearized equation of the system (3.9) at the origin (0, 0) is

where

Notice $\beta_2 = \sqrt{4h+2}$. For $\gamma < \frac{8 (4h+2)}{4h +1}$, $-2 < p(0) < 2$, $q(0) = 1$, so $p^2(0)-4q(0) < 0$, and hence the two roots of $F(\lambda) = 0$ are

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

It is easy to obseve that $|\lambda_{1, 2}(\beta^{**})| = \sqrt {q(\beta^{**})}$ and ${\left. {(|{\lambda _{1, 2}}({\beta ^{**}})|)} \right|_{{\beta ^{**}} = 0}} = \sqrt {q(0)} = 1$. Therefore, a Neimark-Sacker bifurcation probably occurs.

The occurrence of the Neimark-Sacker bifurcation requires the following two conditions to be satisfied:

1) $\quad \Big(\frac{d|\lambda_{1, 2}(\beta^{**})|}{d\beta^{**}}\Big)\bigg|_{\beta^{**} = 0}\neq 0;$

2) $\quad \lambda_{1, 2}^{i}(0)\ne1, i = 1, 2, 3, 4.$

Notice

Obviously $\lambda_{1, 2}^{i}(0)\ne1$ for $i = 1, 2, 3, 4,$ so the two conditions are satisfied.

Second, in order to derive the normal form of the system (3.9), one expands (3.9) in power series up to the third-order term around the origin to get

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

Take matrix

Make a change of variables

then the system (3.11) is changed to the following form:

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

Furthermore,

To determine the stability and direction of the bifurcation curve (closed orbit) for the system (1.7), the discriminating quantity $L$ should be calculated and not to be zero, where

Based on ^[24,25,26], we see that if $L < (>)0,$ then an attracting (a repelling) invariant closed curve bifurcates from the fixed point for $\beta > (<) \beta_{2}$.

The proof of this theorem is complete.

4.   Numerical simulation

In this section, we utilize Matlab to perform numerical simulations to validate the above theoretical analysis through utilizing bifurcation diagrams, phase portraits, maximum Lyapunov expoents, and fractal dimensions of the system (1.7) at the fixed point $E_2$.

Consider the fixed point $E_2$. Vary $\beta$ in the range $(1.4, 1.85)$, and fix $\gamma = 2, h = 0.2$ with the initial value $(x_0, y_0) = (0.4, 0.5)$. Figure 1(a) shows that the existence of a Neimark-Sacker bifurcation at the fixed point $E_{2} = (0.5, 0.5)$ when $\beta = \beta_2 = \sqrt{2.8}\approx 1.6733$. Figure 1(b) describes the spectrum of maximum Lyapunov exponents, which are positive for the parameter $\beta \in (1.4, 1.85)$, which leads to chaos in system (1.7). For this, the interested readers may refer to ^[28] to create an electronic emulator to get immediate results.

The phase portraits associated with Figure 1(a) are drawn in Figure 2. When $\beta$ increases, a circular curve enclosing the fixed point $E_2$ appears.

By choosing a different initial value $(x_0, y_0)$ = (0.52, 0.48) and three same values of $\beta$, the correspending phase portraits are plotted in Figure 3. Figure 2 implies that the closed curve is stable outside, while Figure 3 indicates that the closed curve is stable inside. That is to say, a stable invariant closed curve around the fixed point $E_2$ occurs. This agrees with the conclusion in Theorem 3.3.

5.   Conclusions

In this paper, we consider a predator–prey model with the prey individual behaviour and herd behaviour. By using the semi-discretization method, the continuous system (1.4) is transformed to the discrete system (1.7). Under the given parametric conditions, we demonstrate the existence and stability of three nonnegative fixed points $E_0 = (0, 0)$, $E_{1} = (1, 0)$ and $E_2 = (\frac{1+\sqrt{1+4h\beta^2}}{2\beta^2}, \gamma \frac{1+\sqrt{1+4h\beta^2}}{2\beta^2} (1-\frac{1+\sqrt{1+4h\beta^2}}{2\beta^2}))$. By using the center manifold theory, we determine the existence conditions of transcritical bifurcation and period-doubling bifurcation in the fixed point $E_1$ and the Neimark-Sacker bifurcation at the fixed point $E_2$ of system (1.7). we also derive that $E_2$ is asymptotically stable when $\beta > {\beta _2} = \sqrt{4h+2}$ and unstable when $\beta < \beta_2$. Additionally, the system (1.7) undergoes a Neimark-Sacker bifurcation when the parameter $\beta$ goes through the critical value $\beta_2$. The occurrence for this phenomenon of Neimark-Sacker bifurcation indicates the coexistence of prey and predator when the parameter $\beta = \beta_2$.

Our findings indicate that the proposed discrete model shows a behaviour similar to the one found in the corresponding continuous model ^[27]. In particular, it gives rise to stable populations limit cycles. Ecologically, this means that the suggested response function may be adequate if we want to model the prey herd behaviour that takes place only for a sizable population, namely when the population level settles in a certain threshold (critical value).

Use of AI tools declaration

The authors declare 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 (61473340), the Distinguished Professor Foundation of Qianjiang Scholar in Zhejiang Province (F703108L02), and the Natural Science Foundation of Zhejiang University of Science and Technology (F701108G14). The authors are specially thankful to the anonymous referee for his/her presenting us the electric version of the monograph ^[28], which is really helpful to us.

Conflict of interest

The authors declare that they have no competing interests. All authors contributed equally and significantly in writing this paper. All authors read and approved the final manuscript.

Appendix

We here give a definition and a key Lemma.

Definition 5.1. Let $E(x, y)$ be a fixed piont 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\ne2$;

$(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 and, $|\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 another 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, \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$.