Flip bifurcation and Neimark-Sacker bifurcation in a discrete predator-prey model with Michaelis-Menten functional response

1.   Introduction and preliminaries

In the past several decades, predator-prey interaction has become a hot point of study in biomathematics ^{[1,2,3,4,5,6,7,8,9,10,11,12]}. The differential equation is the main tool to be used in modeling predator-prey interaction when the populations have generation overlap or the numbers (densities) of populations are regarded as varying continuously in time. It can help us understand the interactions of different species within a fluctuating natural environment. Generally speaking, the classical predator-prey model may be written as

where $x$ and $y$ can be expressed as prey and predator population sizes respectively. The function $f(x)$ denotes the growth rate of prey with the absence of predator. $g(x, y)$ represents the amount of prey consumed by per predator per unit time (also called functional response). $h(x, y)$ is on behalf of predator production per capita, and $d$ is the intrinsic death rate of predator ^[1].

Due to the realistic meaning of $f(x)$, one can assume the prey grows logistically with carrying capacity $k$ and growth rate $r$ in the absence of predator, i.e., $f(x) = r(1-\frac{x}{k})$. Besides a good approximation, many scholars ^[2,3] reduce the system (1.1) as

where $a$ is the conversion effciency.

As for the functional response $g(x, y)$, there are many different types. The dynamical complexity of predator-prey system depends on functional response. Liu and Cheng ^[4] proposed a system with square-root functional response, while Bian et al. ^[5] proposed a system with Beddington-DeAngelis funcional response, and so on. In this paper, we discuss the following system

which is called Michaelis-Menten type predator-prey system (or ratio-dependent predator-prey system). The funcional response $\frac{cx}{{my + x}}$ can be expressed as the function of the ratio of prey to predator, where $m$ is the number of prey necessary to achieve one-half of the maximum rate $c$. In this system all of the parameters are positive.

Now nondimensionalize the system (1.3). Let $\frac{x}{k} \to x, rt \to t, \frac{{my}}{k} \to y, \frac{c}{{mr}} \to \alpha, \frac{f}{r} \to \beta, \frac{d}{f} \to \gamma$. Then, we can derive a simplier form of the system (1.3) as follows

This continuous system has been discussed in ^{[6,7,8,9,10,11,12]}.

To be honest, although many methods for continuous systems are mature and have been used to get some interesting results ^[13,14], it is very difficult to solve a complicated differential equation (system) without using a computer. So, we try to use discretization method to derive and study the discrete model of a complicate differential equation (system) so that we can understand the properties of corresponding continuous systems ^{[15,16,17,18,19,20]}. The discrete system related to the system (1.4) has not been investigated yet. In this paper, we select the semi-discretization method, which has better accuracy, to get the discrete version of the system (1.4). To this end, let $[t]$ represent the greatest integer not exceeding $t$. Now, we explore the average change rate of the system (1.4) at integer number points

We know 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 two characteristics as follows:

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

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

One can derive the following system 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)$. Assuming $t\to (n+1)^-$ in the system (1.6) produces

where $\alpha, \beta, \gamma > 0$.

In the sequel, one considers the dynamical properties of the system (1.7). The rest of the paper is organized as follows: in Section 2, one investigates in detail the existence and local stability of the nonnegative fixed points of the system (1.7). Then, the sufficient conditions are formulated for the occurences of flip bifurcation and Neimark-Sacker bifurcation of the system (1.7) in Section 3. Next, numerical simulations are shown to illustrate the results obtained above in Section 4. In the end, some brief conclusions and discussions are stated in Section 5.

2.   Local stability of fixed points

Because of the biological meaning of the system (1.7), we discuss the local stability of its nonnegative fixed points in this section. By letting

it's easy to find that there are two nonnegative fixed points $E_{1} = (1, 0)$ and $E_{2} = (x_{0}, y_{0})$ for $\max \{\frac{\alpha-1}{\alpha}, 0\} < \gamma < 1$ where

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

whose charactertistic polynomial reads as

In order to analyze the prosperties of the fixed points of the system (1.7), one needs the following lemma and definition ^[21,22].

Definition 2.1. Let $E(x, y)$ be a fixed piont of a 2D discrete system with multipliers $\lambda_1$ and $\lambda_2$.

$(i)$ $E(x, y)$ is called sink if $|\lambda_1| < 1$ and $|\lambda_2| < 1$, so, a sink is locally asymptotically stable.

$(ii)$ $E(x, y)$ is called source if $|\lambda_1| > 1$ and $|\lambda_2| > 1$, so, a source is locally asymptotically unstable.

$(iii)$ $E(x, y)$ is called saddle if $|\lambda_1| < 1$ and $|\lambda_2| > 1$ $($or $|\lambda_1| > 1$ and $|\lambda_2| < 1$$)$, .

$(iv)$ $E(x, y)$ is called to be non-hyperbolic if either $|\lambda_1| = 1$ or $|\lambda_2| = 1$.

Lemma 2.2. Let $F(\lambda) = \lambda^2+M\lambda+N$, where $M$ and $N$ 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 $N < 1$;

  $(i.2)$ $\lambda_1 = -1$ and $\lambda_2\ne-1$ if and only if $F(-1) = 0$ and $M\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 $N > 1$;

  $(i.5)$ $\lambda_1$ and $\lambda_2$ are a pair of conjugate complex roots and, $|\lambda_1| = |\lambda_2| = 1$ \setlength\parindent{4em}if and only if $-2 < M < 2$ and $N = 1$;

  $(i.6)$ $\lambda_1 = \lambda_2 = -1$ if and only if $F(-1) = 0$ and $M = 2$.

$(ii)$ If $F(1) = 0,$ in other words $1$ is one root of $F(\lambda) = 0$, then, another root $\lambda$ satisfies $|\lambda| = (<, > )1$ if and only if $|N| = (<, > )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$.

By using Definition 2.1 and Lemma 2.2, the following conclusions can be obtained.

Theorem 2.3. The fixed point $E_1 = (1, 0)$ of the system (1.7) is a saddle.

The proof for Theorem 2.3 is easy and omitted here. Now consider the fixed point $E_2$. For $\max \{\frac{\alpha-1}{\alpha}, 0\} < \gamma < 1$, denote ${\beta _0} = \frac{{2+2\alpha(1-\gamma^2)}}{{\gamma(1-\gamma)[1+\alpha(1-\gamma)]}}$ and ${\beta _1} = \frac{{\alpha (1 - {\gamma ^2}) - 1}}{{\alpha \gamma {{(1 - \gamma)}^2}}}$. Obviously, $\beta _0 > \beta _1$.

Theorem 2.4. Assume $\max \{\frac{\alpha-1}{\alpha}, 0\} < \gamma < 1$, then, $E_2 = (x_0, y_0)$ is a positive fixed point of the system (1.7). Furthermore, the following statements about the fixed point $E_2$ are true.

1) When $\beta < {\beta _0}$,

(a) if $0 < \alpha \le 1$ or $\alpha > 1$ and $\gamma \le \sqrt{\frac{\alpha -1}{\alpha }}$, then, $E_2$ is a stable node;

(b) if $\alpha > 1$ and $\gamma > \sqrt{\frac{\alpha -1}{\alpha }}$, then, for $0 < \beta < \beta_1$, $E_2$ is an unstable node; for $\beta = \beta_1$, $E_2$ is nonhyperbolic; for $\beta_1 < \beta < \beta_0$, $E_2$ is a stable node.

2) When $\beta = {\beta _0}$, then, $E_2$ is non-hyperbolic.

3) When $\beta > {\beta _0}$, then, $E_2$ is a saddle.

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

whose characteristic polynomial can be written as

where

Obviously,

So, when $\beta < {\beta _0}$, $F(-1) > 0$. If $0 < \alpha \le 1$ or $\alpha > 1$ and $\gamma \le \sqrt{\frac{\alpha -1}{\alpha }}$, then $\beta_1 \le 0 < \beta$, so $q < 1$, which reads $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| < 1$ by Lemma 2.2$(i.1)$, therefore, $E_2$ is a stable node.

If $\alpha > 1$ and $\gamma > \sqrt{\frac{\alpha -1}{\alpha }}$, then, for $0 < \beta < \beta_1$, $q > 1$, which reads $\left| {{\lambda _1}} \right| > 1$ and $\left| {{\lambda _2}} \right| > 1$ by Lemma 2.2$(i.4)$, therefore, $E_2$ is an unstable node; for $\beta = \beta_1$, $q = 1$, so, Lemma 2.2$(i.5)$ says $E_2$ is nonhyperbolic; for $\beta_1 < \beta < \beta_0$, $q < 1$, which indicates $\left| {{\lambda _{1, 2}}} \right| < 1$ by Lemma 2.2$(i.1)$, and so $E_2$ is a stable node.

When $\beta = {\beta _0}$, $F(-1) = 0$. Namely, $-1$ is a root of the characteristic polynomial, namely, $E_2$ is non-hyperbolic.

When $\beta > {\beta _0}$, then, $F(-1) < 0$. Lemma 2.2$(i.3)$ says that $\left| {{\lambda _1}} \right| < 1$ and $\left| {{\lambda _2}} \right| > 1$, so $E_2$ is a saddle. The proof is over.

3.   Bifurcation analysis at the fixed point $E_2$

In this section, by using the bifurcation theory and center manifold theorem in ^{[23,24,25,26,27]}, we must pay attention to the flip bifurcation and Neimark-Sacker bifurcation of the system (1.7) at the fixed point $E_2$.

Theorem 2.4 shows that when $\beta = \beta_0$ or $\beta = \beta_1$, the fixed point $E_2$ is non-hyperbolic. Moreover, the dimensional numbers for the stable manifold and the unstable manifold of the fixed point $E_2$ vary when the parament $\beta$ goes through these values, which indicates a bifurcation may occur at each case. In the following analyzation, the paraments comply with $(\alpha, \beta, \gamma) \in {S_{{E_ + }}} = \{ (\alpha, \beta, \gamma) \in R_ + ^3\left| {\alpha > 0, \beta > 0, \max \{\frac{\alpha-1}{\alpha}, 0\} < \gamma < 1 \} } \right.$.

3.1.   Flip bifurcation

When $\beta = \beta_0$, $F(-1) = 0$, which is an indispensable condition for a flip bifurcation to occur. One now explores whether it really exists at the fixed point $E_2$. In fact, the answer is positive.

Theorem 3.1. Suppose the parameters $(\alpha, \beta, \gamma) \in {S_{{E_{_ + }}}}$. If the parameter $\beta$ varies in a small neighborhood of ${\beta _0}$, then, the system (1.7) experiences a flip bifurcation at the fixed point $E_2$.

Proof. Firstly, let $u_n = x_n-x_0, v_n = y_n-y_0$, which transform the fixed point $E_2$ to the origin, and the system (1.7) to

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

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

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

Let

i.e.,

whose three eigenvalues are

with corresponding eigenvectors

where

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

then,

Let

the system (3.2) is changed as follows

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

where $A = \frac{{ \mu - K}}{{4 \alpha\mu\gamma^2 [1 + \alpha (1 - {\gamma })]}}, B = \frac{{ \mu + K}}{{4 \alpha\mu\gamma^2 [1 + \alpha (1 - {\gamma })]}}.$ Next, suppose on the center manifold

where $\rho_3 = \sqrt {X_n^2+\omega_n^2}$. According to

one has

By comparing the corresponding coefficients of terms in the above equation, we get

That's to say the system (3.4) which is restricted to the center manifold can be written as

and

Therefore, we have

i.e.,

According to (21.1.42)–(21.1.46) in ^[24], all of the conditions for the occurrence of a flip bifurcation are satisfied. The proof is over.

3.2.   Neimark-Sacker bifurcation

When $\beta = \beta_1$, a pair of imaginary roots with $|{\lambda _1}| = |{\lambda _2}| = 1$ occur, which is the necessary condition for a Neimark-Sacker bifurcation to occur. That is to say, there may be an occurrence of a Neimark-Sacker bifurcation at the fixed point $E_2$. In fact, one has the following results.

Theorem 3.2. Suppose the parameters $(\alpha, \beta, \delta) \in {S_{{E_{_ + }}}}$ and $\alpha > \frac{1}{1-\gamma^2}$. Let ${\beta _1} = \frac{{\alpha (1 - {\gamma ^2}) - 1}}{{\alpha \gamma {{(1 - \gamma)}^2}}}$. Then, the system (1.7) undergoes a Neimark-Sacker bifurcation at the fixed point $E_2$ if the parament $\beta$ varies in a small neighborhood of the critical value ${\beta _1}$. Moreover, if $L < (>) 0$ in (3.9), then an attracting (repelling) invariant closed curve bifurcates from the fixed point $E_2$ for $\beta > (<) \beta_{1}$.

Proof. Firstly, set a perturbation $\beta^{**}$ of the parameter $\beta$ around $\beta_1$ in the system (3.1), i.e., $\beta^{**} = \beta-\beta_1$ with $0 < |\beta^{**}|\ll1$, the perturbation of the system (3.1) reads

The characteristic equation of the linearized eqution of the system (3.5) at the fixed point (0, 0) is

where

$\alpha > \frac{1}{1-\gamma^2}$ implies $p^2(0)-4q(0) < 0$, so, the two roots of $F(\lambda) = 0$ are

The Neimark-Sacker bifurcation needs to satisfy the following two conditions to occur:

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.$

Due to

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

Secondly, one expands (3.5) into power series up to the third-order term around the origin to get the normal form of the system (3.5) as follows:

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

Then, we can obtain the two roots (eigenvalues) of Eq (3.6) are a pair of conjugate complex as follows:

where $\theta = \sqrt {4{\alpha ^2}{{(1 - \gamma)}^2} - {{[1- \alpha \gamma (1 - \gamma) + {\alpha ^2}(1 - {\gamma ^2})(1 - \gamma)]}^2}}.$

Their corresponding eigenvectors are

where $R = -1 + \alpha \gamma (1 - \gamma) + {\alpha ^2}(1 - {\gamma ^2})(1 - \gamma).$

Let

Transform the variables

then, the system (3.7) is changed to the form as follows

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

where $C = -\frac{R}{{2\theta\alpha^2\gamma^2 (1 - \gamma)}}$, $D = \frac{1}{{2\alpha^2\gamma^2(1 - \gamma)}}$.

One can easily calculate that

To determine the local stability of the closed orbit bifurcated from the Neimark-Sacker bifurcation of the system (3.8), the discriminating quantity $L$ should be calculated and not to be zero ^[25,26,27], where

By calculation we get

Based on the above analysis, it is clear that a Neimark-Sacker bifurcation of the system (1.7) occurs at the fixed point $E_2$ and that the stability of the invariant closed curve bifurcated from the fixed point $E_2$ is determined by the value of $L$. Up to here, the proof for Theorem 3.2 is complete.

4.   Numerical simulation

In this section, numerical simulations are performed to validate above theoretical analysis, including bifurcation diagrams, phase portraits, maximum Lyapunov exponents and fractal dimensions of system (1.7) at the fixed point $E_2$.

Firstly, vary $\beta$ in the range $(9, 10)$, and fix $\alpha = 1.5, \gamma = 0.5$ with the initial value $(x_0, y_0) = (0.3, 0.3)$. Figure 1(a) shows that the existence of a flip bifurcation at the fixed point $E_{2} = (0.25, 0.25)$ when $\beta = \beta_0 = 9.7$, which is in accordance with the result in Theorem 3.1. Figure 1(b) means the spectrum of the maximum Lyapunov exponent. Flip bifurcation may lead to chaos, which makes the system more complex.

Then, let $\beta \in(0.5, 1.5)$. The bifurcation diagram is depicted in Figure 1(c), which illustrates that the fixed point $E_{2}$ is stable for $\beta > \beta_{1} = 0.67$, and unstable when $\beta < \beta_{1}$. Hence, a Neimark-Sacker bifurcation occurs at $E_{2}$ when $\beta = \beta_{1}$, whose multipliers are $\lambda_{1, 2} = \frac{47\pm\sqrt{95}i}{48}$ with $|\lambda_{1, 2}| = 1$. The maximum Lyapunov exponents related to Figure 1(c) are disposed of in Figure 1(d), which exhibits the existence of periodic orbits and chaos as the parameter $\beta$ decreases.

Take the initial values $(x_0, y_0) = (0.3, 0.3)$ in Figure 2. These phase portraits illustrate that the dynamical properties of the fixed point $E_2$ have big changes with the parameter $\beta$ increasing. They change from unstable to stable. What's more, an invariant closed curve around $E_2$ occurs when $\beta = \beta_{1}$. These phenomenons verify the result of Theorem 3.2.

The occurrence of a Neimark-Sacker bifurcation causes the system to jump from stable window to chaotic states through periodic and quasi-periodic states, and trigger a route to chaos.

5.   Discussion and conclusions

This work is concerned with a Michaelis-Menten predator-prey model. Using the semi-discretization method, the system (1.4) is transformed into the system (1.7). Comparing the corresponding continuous system in ^[9], the discrete model has more rich dynamical behaviors. With the given parametric conditions, one demonstrates the existence and local stability of two nonnegative fixed points $E_{1} = (1, 0)$ and $E_{2} = (1-\alpha(1-\gamma), \frac{(1-\gamma)[1-\alpha(1-\gamma)]}{{\gamma}})$. Utilizing the center manifold theorem, one determines the existence conditions of the flip bifurcation and Neimark-Sacker bifurcation of the system (1.7) around the fixed point $E_2$. Especially, $E_2$ is asymptotically stable when $\beta > {\beta _1} = \frac{{\alpha (1 - {\gamma ^2}) - 1}}{{\alpha \gamma {{(1 - \gamma)}^2}}}$ and unstable when $\beta < \beta_1$. So, it is clear that the system (1.7) undergoes a Neimark-Sacker bifurcation when the parameter $\beta$ goes through the critical value $\beta_1$. This phenomenon indicates that the coexistence of prey and predator when the parameter $\beta = \beta_1$. Our results clearly display that the system (1.7) is very sensitive to its parameters: different parameter perturbations will lead to different bifurcations.

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).

Conflict of interest

The authors declare that they have no competing interests.