Full Text Article

Introduction

Predator-prey interaction is one of the most significant phenomena among species [1]. In recent years, many mathematicians and biologists have studied the dynamical behaviors between predator and prey [1], especially using the traditional Lotka-Volterra predator-prey model, which takes the form as follows:

Thanks to the pioneering work of Lotka and Volterra, the study of practical mathematical models in ecology has become a hot topic that has attracted a large number of mathematicians and biologists to join in. In the past few years, many mathematicians have been involved in the dynamical behaviors of predator-prey systems with the theory of dynamical system so that abundent significant results have been yielded [1-6].

In the course of studying the dynamics of predator-prey models, many scholars take into account the effect of the functional response for predator to prey. Yu [7] researched the global asymptotic stability of a predator-prey model with modified Leslie-Gower and Holling-II scheme:

in which x(t)and y(t) denote the densities of the prey and the predator at time t respectively, and are the growth rate of the prey and the predator respectively, r₁ and r₂measures the strength of competition among the individuals of the prey, a₁and₂ are the maximum values that per capita reduction rate of the prey and the predator can attain respectively, k₁andk₂ measure the extent to which the environment provides protections to the prey and the predator respectively [8]. In addition, Yu [7] offered two sufficient conditions on the global asymptotic stability of a positive equilibrium of the system (1.2). After that, Yue [9] considered the dynamics of the following modified Leslie-Gower predator-prey model with Holling-II scheme and a prey refuge:

where mx denotes the part of the refuge protection of the prey, and m∈[0,1] . Yue [9] also found that an increase of the amount of refuge may guarantee the coexistence and attractivity of the two species with no difficulty.

Nowadays, cannibalism, a special phenomenon in nature [10-12], has attracted scholars’ attetion, which means the behavior of consuming the same species. Many species in biology have the phenomenon of cannibalism. For example, some mature organisms eat young individuals, and the stronger ones prey on the weaker ones, etc.. Cannibalism occurs in fish, bird and insect, such as Atlantic salmon, red backed spiders and some copepods [10-12]. Due to the need of energy acquisition and others, the behavior of cannibalism will be widely followed in the whole population [10-12]. This strategy helps adult individuals to preserve energy. Cannibalism of biological species will be helpful to the sustainable survival of organisms to a certain extent. It is universally acknowledged that cannibalism has a quite important effect on the dynamical behaviors of the species.

Scholars once used the bilinear function βxy to describe the cannibalism (refer to [10-16] and the references therein). Till recently the thought of the functional response of predator-prey models was adopted [17-18], and the nonlinear cannibalism model was then proposed.

In 2016 Basheer et al. [17] proposed the predator-prey model with nonlinear prey cannibalism in the following form:

where and represent the densities of prey and predator at time respectively, and the parameters c, c₁,α,β,δ and are nonnegetive constants. Unlike the previous works [13-16], Basheer et al. [17] described the cannibalism in the Holling-II type functional response. The general cannibalism term is C(x)=c.x.x/x+a in the prey equation, in which is the cannibalism rate. This term C(x) is manifestly more appropriate to the reality of ecology and has obvious addition of energy to the cannibalistic prey. The addition leads to an increase in reproduction in the prey by adding a term c₁ to the prey equation. Apparently c₁ < c, as the cannibals need to ingest a lot of prey to produce a new offspring [19]. Scholars concluded that prey cannibalism changes the dynamics of the predator–prey model. The system (1.4) is stable without cannibalism, while it’s unstable with prey cannibalism in the same conditions [17]. After that, Basheer et al. [18] studied the prey-predator model with cannibalism in both predator and prey population and got more comprehensive results.

Model Development

In this paper, we further consider the following predator-prey model with predator cannibalism:

where c₁< c. The meanings of parameters in (1.5) are shown in Table 1, and (cy ²/(y+d) denotes the cannibalism of the predator. In biological sense one assumes that all the parameters are nonnegetive constants.

Without loss of generality, we can assume m=n=1 in the system (1.5). In fact, to do this, the transformation (nx,my,md,α/n)→(x,y,d,α) is sufficient. That is to say, in the sequel, we consider the dynamical properties for the following system:

In general, it is of little possiblility to obtain an exact solution for a complex differential equation or system, so we usually derive its appropriate solution by computer. Thus we should study the corresponding discrete model. For a given differential system, many discretization methods can be utilized, including Euler backward difference method, Euler forward difference method and semidiscretization method, etc.. In this paper, we use the semidiscretization method to derive the discrete model of the system (1.6). For the semidiscretization method, refer also to [20-22, 26-29].

For this purpose, firstly suppose that [t] represents the greatest integer not exceeding t. Then consider the average rate of change of the system (1.6) at integer number points in the following form:

It is quite 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∈[0,+∞) is in possession of the following three characteristics:

1.x(t) and y(t) are continuous on the interval [0,+∞);
2.(dx(t))/dt and (dy(t))/dt exist anywhere when t∈[0,+∞) except for the points t∈{0,1,2,3,⋯};
3.the system (1.6) is true in each interval [n,n+1) with n=0,1,2,3,....

The following system can be obtained by integrating the system (1.7) over the interval [n,t] for any t∈[n,n+1) and n=0,1,2,...:

where x_n=x(n) and y_n=y(n).Letting t→(n+1)^- in the system (1.8) leads to

where all the parameters b,c,c_1,d,α,β>0.

In this paper, our main aim is to consider the dynamics properties of the system (1.9), primarily for its stability and bifurcation. We always assume the space of parameters Ω={(b,c,c_1,d,α,β)∈R₊⁶|c₁< c}.

The rest of this paper is organized 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 occurence of the transcritical bifurcation of the system (1.9). In Section 4, we make some conclusions and discussions about the system (1.9).

Existence and Stability of Fixed Points

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

The fixed points of the system (1.9) satisfy

Considering 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 four nonnegative fixed points E₀ (0,0), E₁ (b/α,0), E₂ (0,(d(c_1-β))/(β+c-c_1 )) for β< c ₁< β+c, and E^* (x^*,y^*) where

for B²≥4αC and x^*< b/α. For the existence of E^* (x^*,y^*), one can also refer to the discussions in its corresponding continuous system in [19].

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

The characteristic polynomial of Jacobian matrix J(E) reads

F(λ)=(λ)²+Pλ+Q,

where

Before we analyze the fixed points of the system (1.9), we recall the following lemma [20-22, 26-29].[lem:201]Let F(λ)=λ²+Pλ+Q, where P and Q are two real constants. Suppose λ₁ and λ₂ are two roots of F(λ)=0. Then the following statements hold.

(i) If F(1)>0, then
(i.1) |λ₁ |<1 and |λ₂ |<1 if and only if F(-1)>0 and Q<1;
(i.2) λ₁=-1 and λ_2≠-1 if and only if F(-1)=0 and P≠2;
(i.3) |λ₁ |<1 and |λ₂ |>1 if and only if F(-1)<0;
(i.4) |λ_1 |>1 and |λ_2 |>1 if and only if F(-1)>0 and Q>1;
(i.5) λ₁ and λ_2 are a pair of conjugate complex roots with |λ₁ |=|λ₂ |=1
if and only if -2< P<2 and Q=1;
(i.6) λ₁=λ₂=-1 if and only if F(-1)=0 and P=2.
(ii) If F(1)=0, namely, 1 is one root of F(λ)=0, then the other root λ
satisfies |λ|=(<,>)1 if and only if |Q|=(<,>)1.
(iii) If F(1)<0, then F(λ)=0 has one root lying in (1,∞). Moreover,
(iii.1) the other root λ satisfies λ<(=)-1 if and only if F(-1)<(=)0;
(iii.2) the other root -1< λ<1 if and only if F(-1)>0.

For the stability of fixed points E₀ (0,0), E ₁ (b/α,0) and E₂ (0,(d(c_1-β))/(β+c-c_1 )), we can get the following Theorems 2.2-2.4, respectively.

[theorem:202] The following statements about the fixed point E_0 (0,0) of the system (1.9) are true.

1. If β>c₁, then E_0 is a saddle.
2. If β=c₁, then E_0 is non-hyperbolic.
3. If β< c₁, then E_0 is a source.

Proof. The Jacobian matrix of the system (1.9) at the fixed point E₀ (0,0) is given by

Obviously,λ₁=e^b > 1 and λ₂= e^-β+c₁ . If β>c₁,then , so is a saddle; if β=c₁,then |λ₂| =1 ,so E_o is non-hyperbolic; if β< c₁,then |λ₂|>1,so E_o is a source. The proof is over.

[theorem:203] The following statements about the fixed point E₁(a/b, 0) of the system (1.9) are true.

1. for b >2, E₁ is an unstable node, i.e., a source;
2. for b=2, E₁ is non-hyperbolic;
3. for 0 < b <2, E₁ is a saddle.

Proof. The Jacobian matrix of the system (1.9) at E₁(a/b, 0) is

We know λ₁ =1-b and λ₂=e ^-β+c1+b/a explicitly.

If β> c₁+b/a,then |λ₂|<1, so we can get the following results:When b>2,|λ₁|>1,so E₁ is a saddle;when b=2,|λ₁|=1,which says E₁ is non-hyperbolic;when 0< b<=2,|λ₁|<1,reading E₁ is a sink.

If β> c₁+b/a,then|λ₂|=1,so is non-hyperbolic.

If β> c₁+b/a,|λ₂|>1 then.Hence,when b>2,|λ₁|>1,so E₁is a source;when b=2,|λ₂|=1, therefore,E₁is non-hyperbolic;when 0< b<2,|λ₁|<1 ,implying E₁ is a saddle.

The proof is complete.

[theorem:204] When β< c₁< β+c, E₂ (0,(d(c₁-β))/(β+c-c₁ )) is a nonnegative fixed point of the system (1.9). Let c₂ ≜(b+d)/b(c_1-β), c₃≜((c₁-β)²)/(c₁-β-2), and d₀≜(2(b+d))/d, then the results about the fixed point E₂ of the system (1.9) are summarized in Table 2.

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

Hereout we obtain λ₁=e^{(b-(d(c_1-β))/(β+c-c_1 ))} and λ₂=1-((β+c-c₁)(c₁-β))/c. Therefore,

Since β< c₁< β+c, here come two cases as follows.

If c≤2, then β< c₁< β+c≤β+2, so 0< c₁-β<2, which says that c(c₁-β-2)<0<(c₁-β)², and hence |λ₂ |<1 holds under all circumstances. Therefore, c<(=,>)(b+d)/b(c₁-β)≜c_2⇔E_2 is a stable node (non-hyperbolic, a saddle).

If c>2, then, for β< c₁≤β+2, c(c₁-β-2)≤0<(c₁-β)². So, |λ² |<1. The type of E₂ is determined by the relation |λ₁ |<(=,>)1 (hence c<(=,>)c_2). For β+2< c₁< β+c,|λ₂ |<(=,>)1⇔c<(=,>)((c₁-β)²)/(c₁-β-2)≜c₃.

Furthermore, one can see

Thereout, we can summarize all the results discussed above in Table 2.

The proof is totally finished.

Bifurcation Analysis

In this section, we use the center manifold theorem and bifurcation theorem to analyze the local bifurcation problems of the fixed points E₀ (0,0), E₁ (b/α,0) and E₂ (0,(d(c_1-β))/(β+c-c_1 )) of the system (1.9). For related bifurcation analysis work for biological systems, refer to the references [26-29] and the references cited therein.

for the Fixed Point E₀(0,0)

Theorem 2.2 shows that a bifurcation of the system (1.9) may occur at the fixed point E₀ in the space of parameters Ω={(b,c,c₁,d,α,β)∈R₊⁶|c₁< c} for β=c₁. We have the following result.

[theorem:301] Consider the parameters (b,c,c₁,d,α,β) in the space Ω. Let β₀=c₁, then the system (1.9) may undergo a 1:1 strong renonance bifurcation at the fixed point E₀ when the parameter β varies in a small neighborhood of the critical value β₀.

Proof. We draw the conclusion through the following analysis.

Giving a small perturbation β^*=β-β₀ with 0<|β^*|<<1 of the parameter β, the system (1.9) is perturbed into

Taylor expanding of the system (3.2) at (x_n,x_n,β_n^*)=(0,0,0) takes the form

Let

then we derive the three eigenvalues of M(E₀) to be λ_1,2,3=1, which displays that a 1:1 strong renonance bifurcation may occur at E₀(0,0). This will be reserved in the futhur discussion.

For the Fixed Point E₁(b/α,0)

Theorem 2.3 shows that a bifurcation of the system (1.9) may occur at the fixed point E₁ (b/α,0) in the space of parameters

Ω₁={(b,c,c₁,d,α,β)∈Ω|β>c₁+bn/α}

or

Ω₂={(b,c,c₁,d,α,β)∈Ω|β< c₁+bn/α}

for b=2. Next one will consider these two cases.

In this case, one has the following result.

[theorem:302] Assume the parameters (b,c,c₁,d,α,β)∈Ω_1. Let b₀=2, then the system (1.9) may undergo a fold-flip bifurcation at the fixed point E₁ when the parameter b varies in a small neighborhood of the critical value b₀.

Proof. Transform the fixed point E¹ (b/α,0) to the origin O(0,0). Giving a small perturbation b^*=b-b₀ with 0<|b^* |≪1 of the parameter b, the system (1.9) is perturbed into

Letting b_(n+1)^*=b_n^*=b^*, the system (3.4) can be written as

Taylor expanding of the system (3.5) at (u_n,v_n,β_n^*)=(0,0,0) obtains

then we derive the three eigenvalues of M(E₁) to be λ₁=-1 and λ₂,3=1. A fold-flip bifurcation may occur at E₁ (b/α,0) when b varies in the neighborhood of the critical value b₀. This will be focused on our further study.

Case II: b₀=2,β< c₁+b/α

In this case, one has the following result, whose proof is similiar to the above Subsection 3.2.1 Case I and omitted here.

[theorem:303] Suppose the parameters in Ω₂. Put b₀=2. Then the system (1.9) may undergo a fold-flip bifurcation at the fixed point E₁ when the parameter b varies in a small neighborhood of the critical value b₀.

For the fixed point E₂ (0,(d(c₁-β))/(β+c-c₁))

Theorem 2.4 shows that a bifurcation of the system (1.9) may occur at the fixed point E₂ (0,(d(c₁-β))/(β+c-c₁ )) when the parameters occur in the spaces Ω_i,i=3,4,...,7 in Table 3, which can be classified into two cases: one is that the parameter c is perturbed around the critical value c₂, and the other is that the parameter c is perturbed around the critical value c₃.

Case I: c₀=c₂

One has the following result.

[theorem:304] Assume the parameters (b,c,c₁,d,α,β)∈Ω_i,i=3,4,5,6,7 in Table 3. Let c₀=c₂=(b+d)/b(c₁-β). If (d(c₁-β))/(b(b+d))≠2, then the system (1.9) undergoes a transcritical bifurcation at the fixed point E2 when the parameter c varies in a small neighborhood of the critical value c₀.

Proof. For convenience of statement, here we take the space of parameters Ω_3={(b,c,c₁,d,α,β)∈Ω|c≤2,β< c₁< β+c} in Table 3 as an example of verification. The proof for parameters in the other spaces is similar and omitted. The proof of this proposition is based on the following analysis.

Take the changes u_n=x_n-0 and v_n=y_n-(d(c₁-β))/(β+c-c₁ ) to translate the fixed point E₂ (0,(d(c₁-β))/(β+c-c₁ )) into the coordinate origin, and the system (1.9) to

Giving a small perturbation c^* of the parameter c, i.e., c^*=c-c₀=c-c₂, with 0<|c_* |<<1, the system (3.7) is perturbed into

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

Thereout we have

Then, according to [23, (21.1.42-46), p507] or [24,25], it is valid for the occurrence of a transcritical bifurcation in the system (1.9) at the fixed point E₂ (0,(d(c₁-β))/(β+c-c₁)) when the parameter c is perturbed around the critical value c₀=c₂=(b+d)/b(c₁-β) and (d(c₁-β))/(b(b+d))≠2 for its parameters in the space Ω₃. Then, a transcritical bifurcation of the system (1.9) at the fixed point E₂ takes place. The proof is over.

Case II: c₀=c₃

Similar to the previous Case I, one can get the following result in this case, whose proof is similar, and hence omitted here.

[theorem:305] Suppose the parameters in the space Ω_5 (or Ω_7) stated in Table 3. Let c₀=((c₁-β)²)/(c₁-β-2),If (d(c₁-β)²)/((c₁-β-2)(b+d)²)≠2, then the system (1.9) undergoes a transcritical bifurcation at the fixed point E₂ when the parameter c varies in a small neighborhood of the critical value c₀.

Conclusion

In this paper, we analyse the dynamical behaviors of a discrete Lokta-Volterra predator-prey system with the predator cannibalism. Given the parameter conditions, we completely formulate the existence and stability of three nonnegative equilibria E_0 (0,0), E_1 (b/α,0) and E_2 (0,(d(c₁-β))/(β+c-c₁)). We also derive the sufficient conditions for the transcritical bifurcation of the system to occur at the fixed points E_2 under variaous different conditions. However, it should be pointed out that the positive equilibria E^* (x^*,y^*) and its complex bifurcations have not been discussed yet, which are worthy being considered in our future study.

Acknowledgements

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

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

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