Period-doubling bifurcation and Neimark-Sacker bifurcation of a discrete predator-prey model with Allee effect and cannibalism

1.   Introduction and preliminaries

Predator-prey relationship is a kind of basic relationships of population interaction in natural ecosystem. Predator-prey model has always been the research focus of biological mathematics, and also one of the most concerned problems by biologists and mathematicians. Recently the most interesting research topic in biological mathematics is the involvement of Allee effect and cannibalism in the population of the prey and predator.

The raise of Allee effect is the most significant phenomenon in the entire biosphere, and it has been considered to be the extremely meaningful factor in both ecology and population dynamics ^[1]. Originally in 1930, the famous ecologist Allee proposed the Allee effect at low-density population, which is considered to be a significant factor in the positive density dependence of low-density population ^[2]. The existence of Allee effect shows that a population must maintain at least one minimum of population size in nature. There are varirties of populations in nature, including insects, birds, mammals and plants ^[3,4,5], in which Allee effect has been widely involved. Recently, the rapid development in the economy and society has proved that adding Allee effect factor to a predator-prey model will lead to an influence on the dynamics of the system, which may result in instability, but it depends on the additional position of Allee effect ^[6].

Cannibalism, as a noteworthy and interesting topic in the predator-prey interaction, plays a crucial role in the dynamics of this interaction. The emphasis on behaviors of cannibalism has been marked in a large number of animals such as non-carnivorous insects, flour beetles, spiders, fish and locusts ^[6,7,8,9,10]. By and large, cannibals and their sufferers are on different stages of life and in various types of categories. This situation regresses the predator-prey interaction in the same species, and the equivalent mathematical model is different from the predator-prey model in structure only for different species ^[11,12]. Polis ^[12] made a great contribution to cannibalism and mentioned about 1300 different species involved in this factor. In the interaction between predator and prey, cannibalism is considered to be a mechanism of natural selection which is actually a common phenomenon ^[13]. In many species, cannibalism occurs when the least resources are provided to the corresponding high population density of the species ^[14]. Although there exists ecological support in both experimental and field work, cannibalism can often be observed in prey populations ^[15,16,17]. However, the experimental work greatly promotes researchers to put forward some ground-breaking ideas in the current situation ^{[18,19,20,21]}.

Deng et al. ^[22] studied a famous Lotka-Volterra model involving cannibalism in predator population, and dealt with the stability of the model, which revealed the effect of cannibalism on stability. Zhang et al. ^[23] proposed a new approach based on non-dimensionalization and applied it to a stage-structured model including predator cannibalism, and also analyzed the dynamics of the model including global stability, supercritical and subcritical Hopf bifurcation and the biological meaning of system parameters.

In order to study a population model correlated with non-overlapping generations, Danca et al. ^[24] investigated the following systems:

where $x_n$ and $y_n$ represent the densities of prey and predator in the $n$th generation respectively, $r$ denotes the intrinsic growth rate of prey, $\alpha$ indicates per capita searching efficiency, and $d$ is the conversion rate of predator. Taking the natural death rate $c$ of predator into consideration, the system (1.1) is changed to the form as follows ^[25]:

Recently, Işık ^[26] further modified the system (1.2) with the addition of Allee effect in prey equation into the following form:

where the constant $m$ is known as Allee constant.

Not long ago, Shabbir et al. ^[27] discussed the complex dynamics of the system (1.2) with the prey cannibalism, i.e., the following system

In fact, the system (1.4) is generated from Euler forward difference method by the following oringinal system:

which is further considered in this paper. Here prey cannibalism has Holling-I type functional response, and we suppose that the prey population $x(t)$ depredates on its own species. The meanings of all parameters in the system (1.5) are shown in Table 1. In the biological sense one assumes that all the parameters are nonnegetive constants.

Without loss of generality, we take $r = 1$ and $\alpha = 1$ in the system (1.5). Actually the transformation $(\frac{\alpha}{r}y, rt, \frac{\beta}{r}, \frac{c}{r}, \frac{d}{r})\rightarrow (y, t, \beta, c, d)$ is sufficient. Denote $a = 1-\beta$. That is to say, in the sequel we consider the dynamical properties for the following system:

Generally speaking, it is very difficult to solve a complicate differential equation (system) without computer. So, one tries to use discretization method to derive and study the discrete model of a complicate differential equation (system) so that one can understand the properties of corresponding continuous system. As for how to derive a discrete version of a complicate differential equation (system), many discretization methods can be utilized, including Euler forward difference method, Euler backward difference method and semidiscretization method, etc.. In this paper, we adopt the semidiscretization method to educe the discrete model of the system (1.6) instead of Euler forward difference method. For the semidiscretization method, refer to ^{[28,29,30,31,32,33,37,38,39,40,41]}.

For this purpose, firstly imagine that $[t]$ represents the maximal integer not exceeding $t$. Then take into account the average rate of change of the system (1.6) at integer points in the following form:

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

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

2) $\frac{\text dx(t)}{\text dt}$ and $\frac{\text dy(t)}{\text dt}$ exist anywhere when $t\in[0, +\infty)$ but for $t = 0, 1, 2, 3, \cdots$;

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

The following system can be derived by incorporating the system (1.7) in 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.8) leads to

where all the parameters

In this paper, we consider the dynamic properties of the system (1.9), primarily for its stability and bifurcation. We always suppose the parameters $(a, c, d, m)\in\Omega$.

In most of the papers that have analyzed prey–predator models, the systems with Allee effect or cannibalism have been widely studied, while those with both Allee effect and cannibalism have not been considered much. Thus the study of the system (1.9) is of significant necessity and novelty to some extent. For related work, see also ^[42,43,44].

The rest of this paper is as follows: In Section 2, we discuss the existence and stability of the fixed points of the system (1.9). In Section 3, we derive the sufficient conditions for the occurrences of period-doubling bifurcation and Neimark-Sacker bifurcation of the system (1.9) at the fixed point $E^*$. In Section 4, we propose numerical simulations to illustrate the occurrence of Neimark-Sacker bifurcation. In Section 5, we make some conclusions and discussions about the system (1.9).

2.   Existence and stability of fixed points

In this section, we analyze the existence of fixed points, and then dissect the local stability of every fixed point of the system (1.9).

The fixed points of the system (1.9) satisfy

Given the biological meanings of the system (1.9), we only take into consideration its nonnegative fixed points. Thus, the system (1.9) has and only has three nonnegative fixed points $E_0(0, 0)$, $E_1(a, 0)$ for $a > 0$, and $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ for $ad-c > 0$.

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

The characteristic polynomial of $J(E)$ reads

where

Prior to studying the stability of fixed points of the system (1.9), we recollect the following lemma ^[32,33].

Lemma 2.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\neq-1$ if and only if $F(-1) = 0$ and $B\neq2$;

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

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

For the stability of the fixed points $E_0(0, 0)$, $E_1(a, 0)$ and $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$, we can get the following Theorems 2.2 - 2.4, respectively.

Theorem 2.2. The following statements about the fixed point $E_0(0, 0)$ of the system (1.9) are true.

1) If $a > 0$, then $E_0$ is a saddle.

2) If $a = 0$, then $E_0$ is non-hyperbolic.

3) If $a < 0$, then $E_0$ is a stable node, i.e., a sink.

Theorem 2.3. The following statements about the fixed point $E_1(a, 0)$ $(a > 0)$ of the system (1.9) are true.

a) If $a > \frac{c}{d}$, then,

1) for $a > 2$, $E_1$ is an unstable node, i.e., a source;

2) for $a = 2$, $E_1$ is non-hyperbolic;

3) for $0 < a < 2$, $E_1$ is a saddle.

b) If $a = \frac{c}{d}$, then $E_1$ is non-hyperbolic.

c) If $a < \frac{c}{d}$, then,

(a) for $a > 2$, $E_1$ is a saddle;

(b) for $a = 2$, $E_1$ is non-hyperbolic;

(c) for $0 < a < 2$, $E_1$ is a stable node, i.e., a sink.

The proofs for Theorems 2.2 and 2.3 are simple and omitted here.

Theorem 2.4. When $ad-c > 0$, $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ is a positive fixed point of the system (1.9). Let $a_1\triangleq\frac{c^2(c+dm+2)-4d(c+dm)}{d[c^2+(c-2)dm]}$ for $d \ne d_1$, $a_2\triangleq\frac{c^2(c+dm+1)}{d[c^2+(c-1)dm]}$ for $d \ne d_2$, $d_1\triangleq\frac{c^2}{m(2-c)}$ for $0 < c < 2$ and $d_2\triangleq\frac{c^2}{m(1-c)}$ for $0 < c < 1$, then the results about the fixed point $E^*$ of the system (1.9) are summarized in Tables 2–3.

Proof. The Jacobian matrix of the system (1.9) at the fixed point $E^*$ can be simplified as follows:

Hereout we obtain its characteristic polynomial

where

Therefore,

Accordingly, we can derive

The scenario $c > 0$ can be divided into three cases as follows.

1) $c\geq2$.

2) $1\leq c < 2$. Then $Q < ( =, > )1\Leftrightarrow a < ( =, > )a_2$.

i) For $c^2+(c-2)dm > 0$, i.e., $d < \frac{c^2}{m(2-c)}\triangleq d_1$,

ii) for $c^2+(c-2)dm = 0$, i.e., $d = d_1$,

iii) for $c^2+(c-2)dm < 0$, i.e., $d > d_1$,

3) $0 < c < 1$. The relationship $F(-1) < ( =, > )0$ is the same as the case $1\leq c < 2$.

i) When $c^2+(c-1)dm > 0$, i.e., $d < \frac{c^2}{m(1-c)}\triangleq d_2$,

ii) when $c^2+(c-1)dm = 0$, i.e., $d = d_2$, $Q < 1$ always holds;

iii) when $c^2+(c-1)dm < 0$, i.e., $d > d_2$,

Thus, all the results discussed above can be summarized in Tables 2 and 3.

The proof is totally finished.

3.   Bifurcation analysis

In this section, we apply the center manifold theorem and local bifurcation theorem ^[34,35,36] to mainly investigate the local bifurcation problems of the system (1.9) at the positive fixed point $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ ($ad-c > 0$), considering the practical biological meaning. For related work, refer to ^{[37,38,39,40,41]}.

3.1.   Period-doubling bifurcation at the fixed point $E^*$

Theorem 2.4 indicates that a bifurcation of the system (1.9) may occur at the fixed point $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ in the space of parameters $\Omega_i, i = 1, 2, ..., 15$ in Table 4 for $a = a_1$. We have the following result.

Theorem 3.1. Consider the system (1.9). Suppose the parameters $(a, c, d, m) \in \Omega_i, i = 1, 2, ..., 15$. Take $a_0 = a_1$. Assume (3.5) and (3.7) hold, then the system (1.9) undergoes a period-doubling bifurcation at the fixed point $E^*$ when the parameter $a$ varies in a small neighborhood of the critical value $a_0$.

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

Step one. Transform the fixed point $E^*$ to the origin $O(0, 0)$ by taking the changes of variables $u_n = x_n-\frac{c}{d}$, $v_n = y_n-\frac{(ad-c)(c+dm)}{cd}$. Then the system (1.9) is changed to

Step two. Giving a small perturbation $a^* = a-a_0 = a-a_1$ with $0 < |a^*|\ll1$, of the parameter $a$ around the critical value $a_0$, the system (3.1) is perturbed into

Letting $a_{n+1}^* = a_n^* = a^*$, the system (3.2) can be written as

Step three. Taylor expansion of the system (3.3) at $(u_n, v_n, a_n^*) = (0, 0, 0)$ to the third order obtains

where $\rho_1 = \sqrt{u_n^2+v_n^2+(a_n^*)^2}$,

Step four. Let

then we can solve the three eigenvalues of $M_1(E^*)$ to get

(in fact, $\lambda_1 = -1$ and $\lambda_2 = 3-E$) and the corresponding eigenvectors

respectively, where

and

Here it must be assumed that $E \ne 2, 4$, namely,

to ensure $|\lambda_2| \ne 1$.

Using the transformation

where

the system (3.4) is turned into the following form

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

Step five. Assume on the center manifold

where $\rho_3 = \sqrt{X_n^2+\delta_n^2}$, then,

and

Comparing the corresponding coefficients of terms with the identical orders in the above center manifold equation produces

Therefore, the system (3.6) restricted to the center manifold takes the following form:

where

Then,

Hence, we have

According to ^[35,36], the occurrence of a period-doubling bifurcation of the system (1.9) at the fixed point $E^*$ needs to satisfy the following conditions

which is equivalent to

The proof is over.

3.2.   Neimark-Sacker bifurcation at the fixed point $E^*$

When $a = a_2 = \frac{c^2(c+dm+1)}{d[c^2+(c-1)dm]}$ for $c \ge 1$ or $0 < c < 1$ and $d \ne d_2$, Theorem 2.4 shows that the system (1.9) has a pair of conjugate complex roots $\lambda_1$ and $\lambda_2$ with $|\lambda_1| = |\lambda_2| = 1$ at the fixed point $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ in the space of parameters $\Omega_i, i = 16, 17, ..., 23$ in Table 5. Here, based on the following analysis, we derive the occurrence of Neimark-Sacker bifurcation according to ^[35,36] in each parameter space in Table 5.

Similiar to the first step in Subsection 3.1, transform the fixed point $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ to the origin $O(0, 0)$, and the system (1.9) to

Give a small change $a^*$ to the parameter $a$, i.e., $a^* = a-a_0 = a-a_2$, then the perturbation of the system (3.8) can be regarded as

The characteristic equation of the linearized equation of the system (3.9) at the coordinate origin point (0, 0) can be read as

where

and

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

Moreover,

which means

Since $p(a^*)\big|_{a^* = 0} = 2-\frac{a_2d^2m+c^2}{d(c+dm)}$ and $q(a^*)\big|_{a^* = 0} = 1$, we have

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

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

In order to discover the normal form of the system (3.9), we reshape the system (3.9) in Taylor expansion to the third-order form around the origin $(0, 0)$ as follows:

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

Let

It is not difficult to get the two eigenvalues of $M_2(E^*)$

where

with the corresponding eigenvectors

Let

Make a transformation of variables

then the system (3.10) is changed into the form as follows:

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

Furthermore,

In order to determine the stability of an invariant closed orbit bifurcated from Neimark-Sacker bifurcation of the system (1.9), we need to calculate the discriminatory quantity $L$ and require it is not equal to zero, where

By calculation we get

Based on the above analysis, we obtain the following result.

Theorem 3.2. Suppose the parameters $a, c, d, m$ in the space $\Omega_i,$ $i = 16, 17,$ $..., 23$. Take $a_2 = \frac{c^2(c+dm+1)}{d[c^2+(c-1)dm)]}$. Let $L$ be defined in (3.12). If the parameter $a$ varies in the small neighborhood of the critical value $a_2$, then a Neimark-Sacker bifurcation occurs at the fixed point $E^*$ of the system (1.9). Furthermore, if $L < (>)0$, then an attracting (repelling) invariant closed curve bifurcates from the fixed point $E^*$ for $a > (<)a_2$.

4.   Numerical simulation

In this section, we employ the bifurcation diagrams, Lyapunov exponents and phase portraits of the system (1.9) at the fixed point $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ for $ad-c > 0$ to corroborate the theoretical results previously derived with some specific parameter values by Matlab software.

Vary $a$ in the interval $(4.2, 4.4)$, and fix the parameter values $c = 3$, $d = 1$ and $m = 8$ with the initial value $(x_0, y_0) = (0.45, 0.35)$. Figure 1(a) illustrates the bifurcation diagram of $(a, x)$-plane. It is quite easy to derive the critical values $a_1 = 4.2941$ and $a_2 = 4.32$, the positive fixed point $E^*(3, 4.84)$, and the eigenvalues of $J(E^*)$ which are $\lambda_{1, 2} = -0.98\pm0.199i$ with $|\lambda_{1, 2}| = 1$. We can also clearly see that the fixed point $E^*$ is unstable when $a < a_1$, stable when $a_1 < a < a_2$ and unstable when $a > a_2$, and that a period-doubling bifurction occurs when $a = a_1 = 4.2941$ and a Neimark-Sacker bifurcation occurs when $a = a_2 = 4.32$. Figure 1(b) depicts the corresponding maximal Lyapunov exponents, which are positive for the parameter $a\in(4.2, 4.4)$, which means the occurrence of chaos in the system (1.9).

Now choose different initial values $(x_0, y_0) = (3, 4.84)$. Then the corresponding phase portraits of the system (1.9) are drawn in Figures 2 and 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, there occurs stable invariant closed curve around the fixed point $E^*$ when $a = a_2 = 4.32$. This agrees to the content of Theorem 3.2.

Biologically, an invariant curve in a prey-predator system is bifurcated from a fixed point, meaning that the prey and the predator can live side by side and produce their own densities. The dynamics of the system on the invariant curve may be periodic or quasi-periodic. A Neimark-Sacker bifurcation in the system indicates that both prey and predator populations can coexist and fluctuate around some mean values of prey survival thresholds, and these fluctuations are stable.

5.   Conclusions and discussion

In this paper, we analyze the dynamical behaviors of a discrete predator-prey system with Allen effect and cannibalism which is derived by using the semidiscretization method and nondmensionalization. Given the parameter conditions, we completely formulate the existence and stability of three nonnegative equilibria $E_0(0, 0)$, $E_1(a, 0)$ for $a > 0$, and $E^*(\frac{c}{d}, \frac{(ad-c)(c+dm)}{cd})$ for $ad-c > 0$. We also derive the sufficient conditions for the period-doubling bifurcation and Neimark-Sacker bifurcation to occur at the fixed points $E^*$ in certain parameter space. Finally some interesting dynamical properties for Neimark-Sacker bifurcation are obtained via numerical simulations.

Neimark-Sacker bifurcation is an important mechanics for one system to produce complicate dynamical behaviors. The occurrence of a Neimark-Sacker bifurcation often causes the system to jump from stable window to chaotic states through periodic and quasi-periodic states, and triggers a route to chaos. Our numerical results just show this point.

Acknowledgements

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 there is no conflicts of interest.