Bifurcation of a discrete predator-prey model with increasing functional response and constant-yield prey harvesting

1.   Introduction and preliminaries

With the continuous development of human society and the continuous progress of civilization, resource consumption and environmental pollution are being increased day by day, and human beings have also been punished by nature, such as frequent occurrences of natural disasters, viruses wreak havoc, etc. So it is very important to find strategies to deal with environmental problems. Mathematical modelling is a force tool to reveal the changing trend of natural environment. More and more scholars use mathematical methods to study ecological balance problem.

Generally speaking, the classical predator-prey model has the following structure:

where $x(t)$ and $y(t)$ represent the population densities of prey and predator in time $t$ respectively, $f(x)$ is the net growth rate of prey without predator, $g(x, y)$ is the consumption rate of prey by predator, $\epsilon$ and $\mu$ are the positive constants respectively representing the conversion rate of captured prey into predator and the mortality of predator. In order to show the crowding effect, when the prey is large, the prey growth rate $f(x)$ in model (1.1) is usually a negative value. The most famous example of $xf(x)$ is the logistic form:

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. In this paper, we assume that $xf(x)$ takes the logistic form given by above (1.2). Consequently, model (1.1) reads as

The behavioral characteristics of the predator species can be reflected by the key element $g(x, y)$, called functional response or nutritional function. Ultimately, the functional response plays an important role in determining different dynamical behaviors, such as steady state, oscillation, bifurcation and chaos phenomenon ^[1]. The functional response $g(x, y)$ in population dynamics (and other disciplines) has several traditional forms:

(ⅰ) $g(x, y)$ depends on $x$ only (meaning $g(x, y) = g(x)$).

$\diamond$ Holling type Ⅰ ^[2,3,4]:

$g(x) = mx;$

$\diamond$ Holling type Ⅱ ^[5,6,7,8]:

$g(x) = \frac{mx}{a+x};$

$\diamond$ Holling type Ⅲ ^[9,10,11,12]:

$g(x) = \frac{mx^2}{a+x^{2}};$

$\diamond$ Holling type Ⅳ ^{[13,14,15,16]}:

$g(x) = \frac{mx}{a+x^{2}}.$

(ⅱ) $p(x, y)$ depends on both $x$ and $y$.

$\diamond$ Ratio-dependent type ^[17]:

$g(x, y) = \frac{mx}{x+ay}$;

$\diamond$ Beddington-DeAngelis type ^[18,19]:

$g(x, y) = \frac{mx}{a+bx+cy}$;

$\diamond$ Hassell-Varley type ^[20,21]:

$g(x, y) = \frac{mx}{y^{\gamma}+ax}, \gamma = \frac{1}{2}, \frac{1}{3}$.

The parameters $m$, $a$, $b$ and $c$ in the above formulas are all positive constants, and they have different biological meanings in different formulas. In order to propose a functional response to show how a group of predators (for example, a group of tuna) search, contact and then hunt a prey or a group of preys, several biological hypotheses were proposed. Based on these assumptions and the logic of Holling ^[22], the hunting cooperation proposed by Cosner, DeAngelis, Ault and Olson ^[23] has a special functional response, as shown below:

here, $C$ is the score of the prey killed in each encountering each predator, $e_{0}$ is the total encountering coefficient between the predator and the prey, and $h$ is the processing time of each prey. It's monotonous. Ryu, Ko and Haque ^[24] introduced this reaction into model (1.3) and obtained the following system:

A common phenomenon in the predator-prey model is called cooperative hunting between predators. This phenomenon makes the encountering rate between the predator and the prey change with the number of predators ^{[25,26,27,28,29,30]}. However, when encountering a gathering of prey, there may be extreme phenomena leading to the eventual extinction of the predator. Therefore, Shang, Qiao, Duan and Miao ^[31] added the constant yield harvest $H$ to the first equation of the model (1.5) to study the arrangement of renewable resources that ensures the coexistence of two species. By the transformations $\overline{t} = rt$, $\overline{x} = \frac{x}{K}$, $\overline{y} = hCe_{0}Ky$, $a = \frac{1}{Ce_{0}h^{2}K^{2}r }$, $b = \frac{\epsilon}{rh}$, $c = \frac{\mu}{r}$ and $\overline{h} = \frac{H}{rK}$, and dropping the bars in the above alphabets, we get the following predator-prey system:

In the system (1.6), we assume that the initial values $(x_0, y_0)$ are positive to ensure that its solution is positive. Obviously, it is very difficult and complicated to directly find an exact solution of the system (1.6), so we consider to find its approximate solution. This motivates us to study the dynamical properties for the discretization version of the system (1.6).

For a given system, there are many discretization methods, including the forward Euler method, the backward Euler method, semi-discretization, and so on. In this paper, we use the forward Euler method to derive the discrete form of the system (1.6). Applying the forward Euler method to the system (1.6), one has

which is written as a map to get the followng system

where $\delta$ is the step size, and $a$, $b$, $c$, $h$ all are positive constants.

The outline of this paper is as follows: In Section 2, we investigate the existence and stability of fixed points of the system (1.8). In Section 3, we use the central manifold theorem and the bifurcation theory to derive some sufficient conditions that ensure the flip bifurcation and saddle-node bifurcation of the system (1.8) to occur. In Section 4, numerical simulation results are provided to not only support theoretical analysis derived but also illustrate new and rich dynamical behaviors of this system.

Before we analyze the fixed points of the system (1.8), we recall the following lemma ^[32,33].

Lemma 1.1. 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 the 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$.

2.   Existence and stability of fixed points

In this section, we first consider the existence of fixed points of the system (1.8) and then analyze the local stability of these fixed points.

The fixed points of the system (1.8) satisfy the following equations

namely,

Considering the biological meanings of the system (1.8), one only takes into account its nonnegative fixed points. Corresponding analysis is as follows:

(ⅰ) if $h > \frac{1}{4}$, then $x(1-x)-\frac{axy^2}{1+xy}-h < 0$ for all nonnegative $x$ and $y$, hence the system (1.8) has no equilibria;

(ⅱ) if $h = \frac{1}{4}$, then $x(1-x)-\frac{axy^2}{1+xy}-h = 0$ if and only if $x = \frac{1}{2}$ and $y = 0$, so, the system (1.8) has a unique predator free equilibrium $A(\frac{1}{2}, 0)$;

(ⅲ) if $0 < h < \frac{1}{4}$, then the system (1.8) has two boundary equilibria $B(\frac{1-\sqrt{1-4h}}{2}, 0)$ and $C(\frac{1+\sqrt{1-4h}}{2}, 0)$, and some positive equilibria may take place. Next, we further analyse this case.

If the system (1.8) has a positive equilibrium, denoted as $(x, y)$, then following (2.2) we have

where $0 < x < 1$, $b > c$ and $0 < h < \frac{1}{4}$.

Let

then

Since $0 < h < \frac{1}{4}$, $f^\prime(x)$ has two unequal real roots $x_1 = \frac{1-\sqrt{1-3h}}{3}$ and $x_2 = \frac{1+\sqrt{1-3h}}{3}$ in the interval $(0, 1)$. On the other hand, we can see that $0 < f(0) < f(1)$ and $0 < x_1 < x_2 < 1$, and that $f(x)$ is increasing for $x\in(0, x_1)\cup(x_2, 1)$ and decreasing for $x\in(x_1, x_2)$.

For the sake of convenient discussion later, let

It is easy to compute $f(x_2) = \frac{c^2(a-a_0)}{b(b-c)}$. So, we have the following results about the positive real roots $x\in(0, 1)$ of $f(x) = 0$:

(ⅰ) if $a > a_0$, then $f(x_2) > 0$, hence $f(x)$ has no positive real root in $(0, 1)$$\Rightarrow$ the system (1.8) has no positive equilibria;

(ⅱ) if $a = a_0$, then $f(x_2) = 0$, hence $f(x)$ has one real root $x_2$ in $(0, 1)$, and it is a double root$\Rightarrow$ the system (1.8) possesses a unique positive equilibrium $E_1(\frac{1+\sqrt{1-3h}}{3}, \frac{c(1-\sqrt{1-3h})}{h(b-c)})$;

(ⅲ) if $a < a_0$, then $f(x_2) < 0$, hence $f(x)$ has two positive roots $x_A$ and $x_B$ in $(0, 1)$, and $0 < x_1 < x_A < x_2 < x_B < 1$$\Rightarrow$ the system (1.8) has two positive equilibria $E_2(x_A, \frac{c}{(b-c)x_A})$ and $E_3(x_B, \frac{c}{(b-c)x_B})$.

Summarizing the above discussions, we obtain the following result.

Theorem 2.1. Consider the system (1.8). Suppose $a_0$ is defined in (2.4). The existence conditions for all nonnegative fixed points ofthe system (1.8) are summarized in the Table 1.

Now we begin to analyze the stability of these fixed points. The Jacobian matrix $J$ of the system (1.8) at a fixed point $E(x, y)$ is presented as follows:

and the characteristic equation of Jacobian matrix $J(E)$ can be written as

where

For the stability of fixed points $A(\frac{1}{2}, 0)$, $B(\frac{1-\sqrt{1-4h}}{2}, 0)$ and $C(\frac{1+\sqrt{1-4h}}{2}, 0)$, we can easily get the following Theorems 2.2–2.4, respectively.

Theorem 2.2. The fixed point $A = (\frac{1}{2}, 0)$ of the system (1.8) is non-hyperbolic.

Theorem 2.3. For $0 < h < \frac{1}{4}$, the boundary fixed point $B = (\frac{1-\sqrt{1-4h}}{2}, 0)$ ofthe system (1.8) occurs. Moreover, the following statements about the fixed point $B$ are true.

1) If $0 < \delta < \frac{2}{c}$, $B$ is a saddle;

2) if $\delta = \frac{2}{c}$, $B$ is non-hyperbolic;

3) if $\delta > \frac{2}{c}$, $B$ is a source.

Theorem 2.4. For $0 < h < \frac{1}{4}$, the boundary fixed point $C = (\frac{1+\sqrt{1-4h}}{2}, 0)$ ofthe system (1.8) occurs. In addition, the following results in the Table 2 arevalid about the fixed point $C$.

For the stability of the positive equilibrium point $E_1(\frac{1+\sqrt{1-3h}}{3}, \frac{c(1-\sqrt{1-3h})}{h(b-c)})$, one will discuss it in the next section.

3.   Bifurcation analysis

In this section, we use the central manifold theorem and bifurcation theory to discuss the flip bifurcation and saddle-node bifurcation at the boundary fixed point $B$ and the positive equilibrium point $E_1$ of the system (1.8).

3.1.   Flip bifurcation

From Theorem (2.3) one can see that, when the parameter $\delta$ goes through the critical value $\delta_0 = \frac{2}{c}$, the dimension numbers of stable and unstable manifolds of the system (1.8) at the fixed point $B$ change. A bifurcation will occur. Again, for $\delta = \delta_0$, one eigenvalue $-1$ appears. So, at this time, the system may produce a flip bifurcation, which is considered in the following, and $\delta$ is chosen as bifurcation parameter. Remember the parameters

Let $X = x-x_{B}, Y = y-y_{B}, \delta^{*} = \delta-\delta_{0}$. We transform the fixed point $B(x_{B}, y_{B})$ to the origin and consider the parameter $\delta^{*}$ as a new independent variable. Thus, the system (1.8) becomes

Taylor expanding of the system (3.1) at $(X, Y, \delta^*) = (0, 0, 0)$ takes the form:

where $\rho_1 = \sqrt {X_n^2+Y_n^2+{\delta_n^*}^2},$

Namely, the system (3.2) is equivalent to the following form:

where

By the center manifold theorem, the stability of $(X, Y) = (0, 0)$ near $\delta^* = 0$ can be determined by studying a one-parameter family of map on a center manifold, which can be written as:

Assume that $h_1^*(Y, \delta^*)$ has the following form:

where $\rho_3 = \sqrt{Y^2+\delta^{*^2}}$. Then the center manifold equation must satisfy

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

Thus the system (3.3) restricted to the center manifold is given by

and

Therefore, one has

According to ^[34], if the nondegenerency conditions $\frac{\partial^3 F^2}{\partial Y^3}|_{(0, 0)}\neq0$ and $\frac{\partial^2 F^2}{\partial Y \partial \delta^*}\bigg|_{(0, 0)} \ne 0$ hold, then the system (1.8) undergoes a flip bifurcation. Obviously, they hold. Therefore, the following result may be derived.

Theorem 3.1. Assume the parameters $(a, b, c, h, \delta) \in S_{E_{+}} = \{(a, b, c, h, \delta) \in R^5_{+}| 0 < a, 0 < h < \frac{1}{4}, 0 < c < b, \delta > 0 \}$. Let $\delta_0 = \frac{2}{c}$, then the system (1.8) undergoes a flip bifurcation at $B(\frac{1-\sqrt{1-4h}}{2}, 0)$ when the parameter $\delta$ varies in a small neighborhood of the critical value $\delta_0$.

3.2.   Saddle-node bifurcation

In the next one considers the saddle-node bifurcation of the system (1.8) at the positive fixed point $E_1(x_0, y_0)$, where $a$ is chosen as bifurcation parameter. The characteristic equation of Jacobian matrix J of the system (1.8) at the positive fixed point $E_1(x_0, y_0)$ is presented as

where $x_0 = \frac{1+\sqrt{1-3h}}{3}$, $y_{0} = \frac{c}{(b-c)x_{0}}$, $0 < h < \frac{1}{4}$ and $b > c$, $p(x_{0})$ and $q(x_{0})$ are given by

Notice $f(1) = 0$ always holds. So, $\lambda_1 = 1$ is a root of $f(\lambda) = 0$. If

then another eigenvalue of the fixed point $E_1(x_{0}, y_{0})$ satisfies

As this time, the system may produce a fold bifurcation, which is considered in the following.

Let $X = x-x_{0}, Y = y-y_{0}, a^{*} = a-a_{0}$, which transform the fixed point $(x_{0}, y_{0})$ to the origin. Consider the parameter $a^{*}$ as a new independent variable, then the system (1.8) becomes

Taylor expanding of the system (3.6) at $(X, a^{*}, Y) = (0, 0, 0)$ obtains

where $\rho_1^* = \sqrt {X_n^2+a_n^{*^{2}}+Y_n^2}$,

Then the system (3.7) is equivalent to the following form:

where

Assume that

Take

then $T^{-1}$ exists.

Under the transformation

the system (3.8) becomes

where $\rho_2^* = \sqrt{X^{2}+{a^*}^{2}+Y^{2}},$

By the center manifold theorem, the stability of $(U, V) = (0, 0)$ near $a_1^* = 0$ can be determined by studying a one-parameter family of map on a center manifold, which can be written as:

Assume that $h_2^*(U, a_1^*)$ has the following form:

where $\rho_3^* = \sqrt{U^2+a_1^{*^2}}$. Then the center manifold equation must satisfy

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

Thus the system (3.10) restricted to the center manifold is given by

where

Therefore, one has

If the condition $\frac{\partial^2 G^*}{\partial U^2}|_{(0, 0)}\neq0$ is true, then the system (1.8) undergoes a saddle-node bifurcation ^[34]. Therefore, we need assume

And the following result may be derived.

Theorem 3.2. Consider the system (1.8). Let $a_0$ be defined in (2.4). Set the parameters $(a, b, c, h, \delta) \in S_{E_{+}} = \{(a, b, c, h, \delta) \in R^5_{+}| 0 < h < \frac{1}{4}, 0 < c < b$,$\sqrt{1-3h}\neq\frac{\delta(3c(b-c)+2b-c)-3nb}{2\delta(2b-c)}, n = 0, -2\}$.

If the conditions (3.9) and (3.11) hold, then the system (1.8) undergoes a saddle-node bifurcation at $E_1(\frac{1+\sqrt{1-3h}}{3}, \frac{c(1-\sqrt{1-3h})}{h(b-c)})$ when the parameter $a$ varies in a small neighborhood of the critical value $a_0$.

Remark. Since the characteristic equation corresponding to the system (3.10) contains double roots $\lambda_1 = \lambda_3 = 1$, the normal form can not be obtained by known routine method. Here we use a special mathematical skill to find the invertible matrix $T$.

4.   Numerical simulation

In this section, we give the bifurcation diagrams of the system (1.8) to illustrate the above theoretical analyses and further reveal some new dynamical behaviors to occur by Matlab software.

First fix the parameter values $a = 3$, $b = 1.1$, $c = 0.6$, $h = 0.16$, let $\delta \in (2, 5)$ and take the initial values $(x_{0}, y_{0}) = (0.2, 0)$ in Figure 1. We can see that there is a stable fixed point for $\delta\in(2, 3.35)$, and a flip bifurcation occurs at $\delta_0 = 3.35$, eventually, period-double bifurcation to chaos. The fixed point $E$ is unstable when $\delta > \delta_{0}$. This agrees to the results stated in Theorem 3.1.

Then fix the parameter values $b = 0.56$, $c = 0.25$, $\delta = 0.187$, $h = 0.12$, and vary $a$ in the range $(0.18, 0.23)$ with the initial value $(x_0, y_0) = (0.6, 1.3)$ in Figure 2. One can see that there is a stable fixed point for $a\in(0.195, 0.205)$, and that a saddle-node bifurcation occurs at $a_0 = 0.2$. When $a < a_{0}$ and is increasing to $a_0$, the fixed point $E_1$ is gradually stable. When $a > a_{0}$, the fixed point $E_1$ is unstable. This agrees to the results stated in Theorem 3.2.

5.   Discussion and conclusions

In this paper, toward a discrete-time predator-prey system of Gause type with constant-yield prey harvesting and a monotonically increasing functional response in $R^2$, we investigate its flip bifurcation and saddle-node bifurcation problems. By using the center manifold theorem and the bifurcation theory, one shows that the flip bifurcation and saddle-node bifurcation of the discrete-time system take place.

We finally present numerical simulations, which not only illustrate the theoretical analysis results, but also find some new properties of the system (1.8)-chaos occurring.

One of the highlights in this paper is to skillfully find an invertible transform to derive the normal form of the flip (fold) bifurcation of the system (1.8), and determine the stability of the closed orbit bifurcated, while it is impossible for one to use routine methods because its two characteristic roots are double so that corresponding invertible matrix does not exist.

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, and the Natural Science Foundation of Zhejiang University of Science and Technology (F701108G14).

Conflict of interest

The authors declare there is no conflicts of interest.