Change language
English Deutsch
Change currency
€ EUR £ GBP $ USD
BY 4.0 license Open Access Published by De Gruyter Open Access August 18, 2022

Bifurcation and chaos in a discrete predator-prey system of Leslie type with Michaelis-Menten prey harvesting

  • Jialin Chen , Zhenliang Zhu , Xiaqing He and Fengde Chen EMAIL logo
From the journal Open Mathematics
https://doi.org/10.1515/math-2022-0054
Download article (PDF)

Abstract

In this paper, a discrete Leslie-Gower predator-prey system with Michaelis-Menten type harvesting is studied. Conditions on the existence and stability of fixed points are obtained. It is shown that the system can undergo fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation by using the center manifold theorem and bifurcation theory. Numerical simulations are presented to illustrate the main theoretical results. Compared to the continuous analog, the discrete system here possesses much richer dynamical behaviors including orbits of period-16, 21, 35, 49, 54, invariant cycles, cascades of period-doubling bifurcation in orbits of period-2, 4, 8, and chaotic sets.

Keywords: discrete predator-prey system; Leslie-Gower; Michaelis-Menten type harvesting; fold bifurcation; flip bifurcation; Neimark-Sacker bifurcation
MSC 2010: 34C25; 92D25; 34D20; 34D40

1 Introduction

Due to its universal existence and importance, the predation relationship between predator and prey is one of the dominant themes in ecology. Since the pioneering works of Lotka and Volterra, their classical models (called Lotka-Volterra predator-prey models) have been reasonably modified to incorporate real biological backgrounds. Among these modifications are Leslie-Gower models, which have been extensively studied (see, for example, [1,2, 3,4] and the references cited therein). In such models, the environmental carrying capacity of the predator species is determined by the abundance of the prey species. A general Leslie-Gower predator-prey model can be written as follows:

d x d t = r x 1 − x K − ϕ ( x ) y , d y d t = s y 1 − h y x ,

where x and y are the densities of the prey and predator, respectively; ϕ is a functional response, and all the parameters are positive constants.

As we know, humans need to exploit natural resources to satisfy their own needs. To ensure the sustainable development of the ecosystem and to maximize economic benefits, harvesting models have therefore attracted the attention of many scholars [1,5,6, 7,8,9, 10,11,12, 13,14,15, 16,17,18, 19,20]. The study of harvesting models can provide suggestions to relevant government agencies on the development of policies to regulate and discipline. Among different forms of harvesting are the constant harvesting and linear harvesting proposed by May et al. [5]. In constant harvesting, the harvest rate is independent of the number of the harvested population, which is not practical. In linear harvesting [7,9,10, 11,12,13, 14,15,19], the harvest rate is proportional to the size of the harvested species, and it has the expression h = q E x . Note that when E or x tends to infinity and the other one is fixed, h ( x ) would tend to infinity. Obviously, this is against the fact that the harvesting capacity and the size of species are limited in reality. To overcome the drawbacks of these two harvesting forms, Clark and Mangel [6] proposed nonlinear harvesting, i.e., Michaelis-Menten type harvesting h ( x ) = q E x m E + n x . If E → ∞ , then h → q x m , the linear harvesting and while if x → ∞ , then h → q E n , the constant harvesting. Therefore, nonlinear harvesting is much more realistic, and it has attracted the attention of many researchers (to name a few, see [1,8,11,12,14,16, 17,18,20]).

The aforementioned works mainly focus on continuous predator-prey systems and obtained results include stability, bifurcation, limit cycles, and so on [21,22,23,24,25]. However, when species have nonoverlapping generations or their sizes are too small, discrete models described by difference equations are more appropriate than continuous-time ones. Over the past few decades, dynamic behaviors of discrete predator-prey systems have been widely studied, some of which can be found in [26,27,28, 29,30,31, 32,33,34, 35,36,37, 38,39,40] and the references cited therein. All these investigations have demonstrated that discrete systems tend to have more complex dynamic behaviors than continuous ones. In particular, He and Lai [28] studied a discrete Lotka-Volterra predator-prey system of Leslie-Gower with Holling type II functional response, whereas Ren et al. [41] proposed a system with Crowley-Martin functional response. In [36], Zhao and Yan investigated the aforementioned system with a modified Holling-Tanner functional response. In addition, many scholars have considered adding harvesting terms to these systems [42,43, 44,45]. For example, Hu and Cao [42] studied a discrete system of Holling type II functional response and Leslie type with constant-yield prey harvesting; Zhu et al. [44] proposed a discrete May type cooperative model incorporating Michaelis-Menten type harvesting. These papers mainly investigated the dynamical behaviors of bifurcation and chaos of the systems.

The aforementioned discussion motivates us to study a discrete Leslie-Gower predator-prey model with nonlinear harvesting. Precisely, the model is based on the following continuous Leslie-Gower predator-prey model with Michaelis-Menten type prey-harvesting and linear functional response proposed by Gupta et al. [46],

(1.1) d x d t = r x 1 − x K − α x y − q E x m E + n x , d y d t = s y 1 − h y x ,

where x and y denote the densities of prey and predator at time t , respectively, and r , K , α , q , E , m , n , s , and h are all positive constants. Here, r is the intrinsic growth rate of the prey and K represents the carrying capacity of the prey in the absence of predators; α is the maximum predation rate of the predator and s is the intrinsic grow rate of the predator; q stands for the harvesting coefficient and E is the harvesting effort of the prey species; h represents the number of prey required to provide one predator at equilibrium when y equals to x h .

For the sake of simplicity, we first rescale system (1.1) by introducing

x ¯ = x K , y ¯ = h K y , t ¯ = r t .

Dropping the bars, system (1.1) becomes

(1.2) d x d t = x 1 − x − b y − e c + x , d y d t = p y 1 − y x ,

where b = α K h r , c = m E n K , e = q E r n K , and p = s r are positive constants. Then applying the forward Euler scheme to system (1.2) yields the model to be studied in this paper,

(1.3) x → x + δ x 1 − x − b y − e c + x , y → y + δ p y 1 − y x ,

where δ is the step size. The aim of this paper is to study the dynamical behaviors of system (1.3) including the stability of fixed points and bifurcation phenomena.

The rest of the paper is arranged as follows. We first study the existence and stability of fixed points in Section 2. Then, in Section 3, we show that system (1.3) can undergo fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation under appropriate conditions on parameter values. The feasibility of the main theoretic results is illustrated with numerical simulations in Section 4. The paper ends with a brief conclusion.

2 The existence and stability of fixed points

2.1 The existence of fixed points

We start with the existence of fixed points of system (1.3). Note that fixed points of (1.3) are solutions of

(2.1) x = x + δ x 1 − x − b y − e c + x , y = y + δ p y 1 − y x .

First, we consider boundary fixed points, that is, x ≠ 0 and y = 0 . It follows easily from (2.1) that x satisfies

(2.2) x 2 + ( c − 1 ) x + e − c = 0 .

Let Δ 1 denotes the discriminant, namely,

(2.3) Δ 1 = ( c + 1 ) 2 − 4 e .

A necessary condition on the existence of boundary fixed points is Δ 1 ≥ 0 , i.e., ( c + 1 ) 2 ≥ 4 e . In this case, the roots of (2.2) can be written as follows:

(2.4) x 1 = 1 − c − Δ 1 2 , x 2 = 1 − c + Δ 1 2 .

Specifically, when Δ 1 = 0 , x 1 , and x 2 collide, and we denote it as

(2.5) x 3 = 1 − c 2 .

Taking into consideration of seeking positive solutions of (2.2), we easily obtain the following result on the existence of boundary fixed points of (1.3).

Theorem 2.1

The following statements on the boundary fixed points of (1.3) hold.

  1. If e > c , c < 1 , and ( c + 1 ) 2 > 4 e , then there are two boundary fixed points E 1 ( x 1 , 0 ) and E 2 ( x 2 , 0 ) ;

  2. There is only the boundary fixed point E 3 ( x 3 , 0 ) if c < 1 and ( c + 1 ) 2 = 4 e ;

  3. There is only the boundary fixed point E 2 ( x 2 , 0 ) if either e < c or e = c < 1 .

Here, x 1 and x 2 are given in (2.4), and x 3 is given in (2.5).

Now, we consider positive fixed points, which satisfy

(2.6) 1 − x − b y − e c + x = 0 , y = x

by (2.1). Substituting y = x into the first equation of system (2.6) produces

(2.7) ( b + 1 ) x 2 + ( b c + c − 1 ) x + e − c = 0 .

Similarly as for the discussion of boundary fixed points, let

Δ 2 = [ c ( b + 1 ) + 1 ] 2 − 4 e ( b + 1 ) .

Then we require Δ 2 ≥ 0 or [ c ( b + 1 ) + 1 ] 2 ≥ 4 e ( b + 1 ) . Under this condition, when Δ 2 > 0 , equation (2.7) has two solutions,

(2.8) x 1 ∗ = 1 2 1 b + 1 − c − Δ 2 b + 1 , x 2 ∗ = 1 2 1 b + 1 − c + Δ 2 b + 1 ,

and when Δ 2 = 0 , it has a unique root,

(2.9) x 3 ∗ = 1 2 1 b + 1 − c .

Simple calculations give the existence of positive fixed points of (1.3), which is summarized below.

Theorem 2.2

The following statements are valid for the positive fixed points of (1.3).

  1. If e > c , 1 b + 1 > c , and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) , there are two positive fixed points E 1 ∗ ( x 1 ∗ , y 1 ∗ ) and E 2 ∗ ( x 2 ∗ , y 2 ∗ ) ;

  2. If 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 = 4 e ( b + 1 ) , there is only one positive fixed point E 3 ∗ ( x 3 ∗ , y 3 ∗ ) ;

  3. If either ( e = c and 1 b + 1 > c ) or e < c , there is only one positive fixed point E 2 ∗ ( x 2 ∗ , y 2 ∗ ) .

Here, x i ∗ are given in (2.8) or (2.9), and y i ∗ = x i ∗ , i = 1 , 2 , 3 .

2.2 The stability of fixed points

In this section, we study the (local) stability of the fixed points obtained earlier. Denote the Jacobian matrix of system (1.3) evaluated at a fixed point E ( x , y ) by J ( E ) , where

J ( E ) = J 11 J 12 J 21 J 22

with

J 11 = 1 + δ 1 − x − b y − e c + x + δ x − 1 + e ( c + x ) 2 , J 12 = − δ b x , J 21 = δ p y 2 x 2 , J 22 = 1 + δ p 1 − y x − δ p y x .

We write the characteristic equation of J ( E ) as F ( λ ) = λ 2 + B λ + C = 0 . E can be classified according to the two roots λ 1 and λ 2 of F ( λ ) as follows (see, for example, [47, pp. 27]).

Definition 2.1

The fixed point E is called

  1. a sink if ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ < 1 , and it is locally asymptotically stable;

  2. a source if ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ > 1 , and it is unstable;

  3. a saddle if either ( ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ < 1 ) or ( ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ > 1 );

  4. nonhyperbolic if either ∣ λ 1 ∣ = 1 or ∣ λ 2 ∣ = 1 .

The following result gives sufficient and necessary conditions on the moduli of λ 1 and λ 2 with respect to 1 in the case of F ( 1 ) > 0 .

Lemma 2.1

[47, Lemma 2] Assume that F ( 1 ) > 0 . Then

  1. ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ < 1 if and only if F ( − 1 ) > 0 and C < 1 ;

  2. ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ > 1 if and only if F ( − 1 ) > 0 and C > 1 ;

  3. ( ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ < 1 ) or ( ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ > 1 ) if and only if F ( − 1 ) < 0 ;

  4. λ 1 = − 1 and ∣ λ 2 ∣ ≠ 1 if and only if F ( − 1 ) = 0 and B ≠ 0 , 2;

  5. λ 1 and λ 2 are conjugate complex roots and ∣ λ 1 ∣ = ∣ λ 2 ∣ = 1 if and only if B 2 − 4 C < 0 and C = 1 .

We mention that F ( 1 ) > 0 and F ( − 1 ) = 0 imply that B ≠ 0 . Therefore, B ≠ 0 is redundant in statement 4 of the aforementioned lemma, which will be dropped in the coming discussion.

Similarly, we have the following result for the case of F ( 1 ) < 0 .

Lemma 2.2

Assume that F ( 1 ) < 0 . Then

  1. ∣ λ 1 ∣ > 1 and  ∣ λ 2 ∣ > 1 if and only if F ( − 1 ) < 0 ;

  2. ( ∣ λ 1 ∣ > 1 and  ∣ λ 2 ∣ < 1 ) or ( ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ > 1 ) if and only if  F ( − 1 ) > 0 ;

  3. λ 1 = − 1 and  ∣ λ 2 ∣ ≠ 1 if and only if F ( − 1 ) = 0 .

Proof

Note that F ( λ ) as a quadratic with the leading coefficient being positive. When F ( 1 ) < 0 , we know that F ( λ ) = 0 has two real distinct roots λ 1 and λ 2 on the two sides of 1. Without loss of generality, say λ 1 < 1 < λ 2 . Then the results follow immediately from the intermediate value theorem.□

Theorem 2.3

Under the conditions on the existence of boundary fixed points of (1.3) in Theorem 2.1,

  1. E 1 ( x 1 , 0 ) is always a source and it is unstable;

  2. E 2 ( x 2 , 0 ) is

    1. a saddle if δ < 2 ( c + x 2 ) Δ 1 x 2 ;

    2. a source if δ > 2 ( c + x 2 ) Δ 1 x 2 ;

    3. nonhyperbolic if δ = 2 ( c + x 2 ) Δ 1 x 2 ;

  3. E 3 ( x 3 , 0 ) is always nonhyperbolic.

Proof

At a boundary fixed point E i ( i = 1 , 2 , 3 ), the Jacobian matrix reduces to

J ( E i ) = 1 + δ x i − 1 + e ( c + x i ) 2 − δ b x i 0 1 + δ p

as 1 − x i = e c + x i . Thus, the eigenvalues of J ( E i ) are λ 1 = 1 + δ x i − 1 + e ( c + x i ) 2 and λ 2 = 1 + δ p > 1 .

  1. Since x 1 = 1 2 ( 1 − c − Δ 1 ) , we obtain

    ( c + x 1 ) 2 − e = 1 4 ( 1 + c − Δ 1 ) 2 − e = 1 2 ( Δ 1 − Δ 1 − c Δ 1 ) = Δ 1 2 [ Δ 1 − ( 1 + c ) ] < 0

    or ( c + x 1 ) 2 < e . It follows that λ 1 > 1 . According to Definition 2.1, E 1 is a source and it is unstable when it exists.

  2. It follows from x 2 = 1 − c + Δ 1 2 , we can obtain ( c + x 2 ) 2 − e = Δ 1 2 ( Δ 1 + 1 + c ) > 0 , and hence, λ 1 < 1 . Furthermore, since x 2 satisfies 1 − x 2 − e c + x 2 = 0 , we have

    δ x 2 − 1 + e ( c + x 2 ) 2 = δ x 2 − 1 + 1 − x 2 c + x 2 = − δ Δ 1 x 2 c + x 2 .

    Therefore, if δ < 2 ( c + x 2 ) Δ 1 x 2 holds, then ∣ λ 1 ∣ < 1 , and hence, E 2 is a saddle. This proves (ii)(a). In a similar way, we can prove (ii)(b) and (ii)(c).

  3. It is easy to obtain that ( c + x 3 ) 2 = e since Δ 1 = 0 . Hence, the eigenvalues of J ( E 3 ) are λ 1 = 1 and λ 2 = 1 + δ p > 1 . According to Definition 2.1, E 3 ( x 3 , 0 ) is always nonhyperbolic.

This completes the proof.□

Now, the characteristic equation of the Jacobian matrix J of system (1.3) at each positive fixed point E i ∗ ( i = 1 , 2 , 3 ) is

F ( λ ) ≜ λ 2 − ( 2 + G δ ) λ + ( 1 + G δ + H p δ 2 ) = 0 ,

where

G = x i ∗ e ( c + x i ∗ ) 2 − 1 − p , H ( x i ∗ ) = x i ∗ ( b + 1 ) − e ( c + x i ∗ ) 2 .

Thus,

F ( 1 ) = H p δ 2 , F ( − 1 ) = 4 + 2 G δ + H p δ 2 .

We use Lemmas 2.1 and 2.2 to investigate the stability of the positive fixed points one by one.

Theorem 2.4

Assume that e > c , 1 b + 1 > c , and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) . Then the positive fixed point E 1 ∗ in Theorem 2.2 is

  1. a saddle if 0 < δ < − G + G 2 − 4 H p H p ;

  2. a source if δ > − G + G 2 − 4 H p H p ;

  3. nonhyperbolic if δ = − G + G 2 − 4 H p H p .

Proof

We first determine the sign of F ( 1 ) , or equivalently the sign of H . In fact, since E 1 ∗ ( x 1 ∗ , y 1 ∗ ) satisfies (2.7), we obtain

H ( x 1 ∗ ) = x 1 ∗ ( b + 1 ) − e ( c + x 1 ∗ ) 2 = x 1 ∗ ( b + 1 ) − 1 − ( b + 1 ) x 1 ∗ c + x 1 ∗ = x 1 ∗ c + x 1 ∗ [ 2 ( b + 1 ) x 1 ∗ + c ( b + 1 ) − 1 ] .

Noting x 1 ∗ = y 1 ∗ = 1 2 1 b + 1 − c − Δ 2 b + 1 , we have H ( x 1 ∗ ) = − Δ 2 x 1 ∗ c + x 1 ∗ < 0 , and hence, F ( 1 ) < 0 .

As H ( x 1 ∗ ) < 0 , F ( − 1 ) = 0 as a quadratic equation of δ has a negative solution − G − G 2 − 4 H p H p and a positive solution − G + G 2 − 4 H p H p . If 0 < δ < − G + G 2 − 4 H p H p then F ( − 1 ) > 0 . By Lemma 2.2(ii), E 1 ∗ is a saddle. This proves (i). (ii) and (iii) can be proved similarly by using Lemma 2.2(i) and (iii), respectively.□

For the positive fixed point E 2 ∗ ( x 2 ∗ , y 2 ∗ ) , a similar calculation as for H ( x 1 ∗ ) will produce H ( x 2 ∗ ) = Δ 2 x 2 ∗ c + x 2 ∗ > 0 , and hence, F ( 1 ) > 0 . With the assistance of Lemma 2.1, we easily obtain

Theorem 2.5

Suppose one of the following three conditions holds,

  1. e > c , 1 b + 1 > c , and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) ;

  2. e = c and 1 b + 1 > c ;

  3. e < c .

Then the positive fixed point E 2 ∗ of system (1.3) is

  1. a sink if one of the following conditions holds,

    1. − 2 H p ≤ G < 0 and 0 < δ < − G H p ;

    2. G < − 2 H p and 0 < δ < − G − G 2 − 4 H p H p ;

  2. a source if one of the following conditions holds,

    1. − 2 H p ≤ G < 0 and δ > − G H p ;

    2. G < − 2 H p and δ > − G + G 2 − 4 H p H p ;

    3. G ≥ 0 ;

  3. a saddle if G < − 2 H p and − G − G 2 − 4 H p H p < δ < − G + G 2 − 4 H p H p ;

  4. nonhyperbolic if one of the following conditions holds:

    1. G < − 2 H p , δ = − G ± G 2 − 4 H p H p , and δ ≠ − 4 G ;

    2. − 2 H p < G < 0 and δ = − G H p .

Finally, from Theorem 2.2, the positive fixed point E 3 ∗ exists when [ c ( b + 1 ) + 1 ] 2 = 4 e ( b + 1 ) and 1 b + 1 > c . It follows that F ( 1 ) = 0 , and hence,

Theorem 2.6

When [ c ( b + 1 ) + 1 ] 2 = 4 e ( b + 1 ) and 1 b + 1 > c , the positive fixed point E 3 ∗ of (1.3) is always nonhyperbolic.

3 Bifurcation analysis

In this section, we analyze different bifurcation types at fixed points of system (1.3) by using the center manifold theorem [48] and bifurcation theory [49,50]. We begin with the fold bifurcation.

3.1 Fold bifurcation

Recall from Theorem 2.2.2 that if

(3.1) e ∗ = [ c ( b + 1 ) + 1 ] 2 4 ( b + 1 ) , 1 b + 1 > c ,

then system (1.3) only has one positive fixed point E 3 ∗ ( x 3 ∗ , y 3 ∗ ) and the eigenvalues of the Jacobian matrix J ( E 3 ∗ ) are λ 1 = 1 and λ 2 = 1 + G δ . Suppose that

(3.2) e ∗ ≠ ( x 3 ∗ + p ) ( x 3 ∗ + c ) 2 x 3 ∗ , e ∗ ≠ ( δ p + δ x 3 ∗ − 2 ) ( x 3 ∗ + c ) 2 δ x 3 ∗ .

Then, ∣ λ 2 ∣ ≠ 1 .

Let u = x − x 3 ∗ , v = y − y 3 ∗ , and η = e − e ∗ . Then system (1.3) can be rewritten as follows:

(3.3) u η v → a 1 u + a 2 η + a 3 v + a 11 u 2 + a 12 u η + a 13 u v + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ η ∣ ) 3 ) η b 1 u + b 3 v + b 11 u 2 + b 13 u v + b 33 v 2 + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ η ∣ ) 3 ) ,

where

a 1 = 1 + δ x 3 ∗ e ∗ ( c + x 3 ∗ ) 2 − 1 , a 2 = − δ x 3 ∗ c + x 3 ∗ , a 3 = − δ b x 3 ∗ , a 11 = δ e ∗ c ( c + x 3 ∗ ) 3 − 1 , a 12 = − δ c ( c + x 3 ∗ ) 2 , a 13 = − δ b , b 1 = δ p , b 3 = 1 − δ p , b 11 = − δ p x 3 ∗ , b 13 = 2 δ p x 3 ∗ , b 33 = − δ p x 3 ∗ .

We choose

T 1 = a 2 ( a 1 − λ 2 ) a 2 λ 2 − b 3 0 1 − λ 2 0 a 2 b 1 0 b 1 ,

which is invertible. Then with the transformation

u η v = T 1 x ˜ η 1 y ˜ ,

we transform (3.3) into

(3.4) x ˜ η 1 y ˜ → 1 1 0 0 1 0 0 0 λ 2 x ˜ η 1 y ˜ + ϕ ( x ˜ , y ˜ , η 1 ) 0 ψ ( x ˜ , y ˜ , η 1 ) ,

where

ϕ ( x ˜ , y ˜ , η 1 ) = a 11 b 1 + b 11 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) b 1 u 2 + a 12 a 2 ( 1 − λ 2 ) u η + a 13 b 1 + b 13 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) b 1 u v + b 33 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) b 1 v 2 + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ η ∣ ) 3 ) , ψ ( x ˜ , y ˜ , η 1 ) = b 11 ( a 1 − λ 2 ) − a 11 b 1 ( 1 − λ 2 ) b 1 u 2 − a 12 1 − λ 2 u η + b 13 ( a 1 − λ 2 ) − a 13 b 1 ( 1 − λ 2 ) b 1 u v + b 33 ( a 1 − λ 2 ) ( 1 − λ 2 ) b 1 v 2 + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ η ∣ ) 3 ) , u = a 2 ( a 1 − λ 2 ) x ˜ + a 2 η 1 + ( λ 2 − b 3 ) y ˜ , η = ( 1 − λ 2 ) η 1 , v = a 2 b 1 x ˜ + b 1 y ˜ .

By the center manifold theory, in a small neighborhood of η 1 = 0 , there exists a center manifold W c ( 0 ) of (3.4) at the fixed point ( x ˜ , y ˜ ) = ( 0 , 0 ) , which can be represented as follows:

W c ( 0 ) = { ( x ˜ , y ˜ , η 1 ) ∈ R 3 ∣ y ˜ = h ( x ˜ , η 1 ) , h ( 0 , 0 ) = 0 , D h ( 0 , 0 ) = 0 } ,

where x ˜ and η 1 are sufficiently small. Assume that the expression of h is

(3.5) h ( x ˜ , η 1 ) = m 1 x ˜ 2 + m 2 x ˜ η 1 + m 3 η 1 2 + O ( ( ∣ x ˜ ∣ + ∣ η 1 ∣ ) 3 ) ,

which must satisfy

(3.6) h ( x ˜ + η 1 + ϕ ( x ˜ , h ( x ˜ , η 1 ) , η 1 ) , η 1 ) = λ 2 h ( x ˜ , η 1 ) + ψ ( x ˜ , h ( x ˜ , η 1 ) , η 1 ) .

Substituting (3.5) into (3.6) and comparing the coefficients of terms like x ˜ k η 1 l in the resultant, we obtain

m 1 = a 2 2 ( a 1 − λ 2 ) b 1 ( 1 − λ 2 ) 2 { ( a 1 − λ 2 ) [ ( a 1 − λ 2 ) b 11 − a 11 b 1 + b 13 b 1 ] − a 13 b 1 2 + b 33 b 1 2 } , m 2 = 2 a 2 2 ( a 1 − λ 2 ) ( a 1 − λ 2 ) b 11 − a 11 b 1 b 1 ( 1 − λ 2 ) 2 − a 12 a 2 ( a 1 − λ 2 ) 1 − λ 2 + a 2 2 [ ( a 1 − λ 2 ) b 13 − a 13 b 1 ] ( 1 − λ 2 ) 2 − 2 m 1 1 − λ 2 , m 3 = a 2 2 [ ( a 1 − λ 2 ) b 11 − a 11 b 1 ] b 1 ( 1 − λ 2 ) 2 − a 12 a 2 1 − λ 2 − m 1 + m 2 1 − λ 2 .

Therefore, the map (3.4) restricted to the center manifold can be written as follows:

F 1 : x ˜ → x ˜ + η 1 + k 1 x ˜ 2 + k 2 x ˜ η 1 + k 3 η 1 2 + O ( ( ∣ x ˜ ∣ + ∣ η 1 ∣ ) 3 ) ,

where

(3.7) k 1 = a 2 2 ( a 1 − λ 2 ) ( a 1 − λ 2 ) a 11 b 1 + b 11 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) b 1 + a 13 b 1 + b 13 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) + a 2 b 1 b 33 ( b 3 − λ 2 ) 1 − λ 2 , k 2 = ( a 1 − λ 2 ) 2 a 2 2 a 11 b 1 + 2 a 2 2 b 11 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) b 1 + a 12 + a 2 2 a 13 b 1 + a 2 2 b 13 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) , k 3 = a 2 2 a 11 b 1 + a 2 2 b 11 ( b 3 − λ 2 ) a 2 ( 1 − λ 2 ) b 1 + a 12 .

Since F 1 ( 0 , 0 ) = 0 , ∂ F 1 ∂ x ˜ ( 0 , 0 ) = 1 , ∂ F 1 ∂ η 1 ( 0 , 0 ) = 1 , ∂ 2 F 1 ∂ x ˜ 2 ( 0 , 0 ) = 2 k 1 , and ∂ 2 F 1 ∂ x ˜ η 1 ( 0 , 0 ) = k 2 , we obtain the following result.

Theorem 3.1

In addition to (3.1) and (3.2), suppose that k 1 ≠ 0 or k 2 ≠ 0 . Then system (1.3) undergoes a fold bifurcation at E 3 ∗ , where k 1 and k 2 are given by (3.7). Moreover, the fixed points E 1 ∗ and E 2 ∗ bifurcate from E 3 ∗ for e < e ∗ , coalesce at E 3 ∗ for e = e ∗ , and disappear for e > e ∗ .

3.2 Flip bifurcation

Next we discuss the flip bifurcations of system (1.3).

Let

F A = ( c , e , δ ) ∣ δ = 2 ( c + x 2 ) Δ 1 x 2 , c , e > 0 ,

where x 2 and Δ 1 are given by (2.3) and (2.4), respectively. System (1.3) can undergo flip bifurcation at the boundary fixed point E 2 ( x 2 , 0 ) when parameters belong to F A . Since a center manifold of system (1.3) at E 2 is y = 0 and system (1.3) restricted to it is the logistic model:

x → f ( x ) = x + δ x 1 − x − e c + x .

Its nontrivial fixed point is x 2 = 1 − c + Δ 1 2 . It is easy to see that

f ′ ∣ x 2 = 1 + δ x 2 − 1 + e ( c + x 2 ) 2 < 0

when parameters vary in a small neighborhood of F A . Thus, flip bifurcation can occur (see Figure 2). In this case, the predator species becomes extinct, and the prey species undergoes the flip bifurcation to chaos by choosing the appropriate value of the bifurcation parameter δ .

Since E 1 ∗ is always unstable, we now focus on flip bifurcation at E 2 ∗ due to the biological significance. Here, we choose δ as the bifurcation parameter.

From Theorem 2.5(iv(a)), we can easily obtain that one of the eigenvalues of J ( E 2 ∗ ) is − 1 and the other one is neither 1 nor − 1 . We rewrite the conditions in Theorem 2.5(iv(a)) as the following two sets:

F B 1 = ( b , c , e , p , δ ) (i) b , c , e , p > 0 , δ = − G − G 2 − 4 H p H p , G < − 2 H p , δ ≠ − 4 G , (ii) e > c > 0 , 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) , or e = c and 1 b + 1 > c , or e < c

and

F B 2 = ( b , c , e , p , δ ) (i) b , c , e , p > 0 , δ = − G + G 2 − 4 H p H p , G < − 2 H p , δ ≠ − 4 G , (ii) e > c > 0 , 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) , or e = c and 1 b + 1 > c , or e < c .

We shall see that flip bifurcation may undergo when parameters belong to F B 1 or F B 2 . We only provide the detail on discussing parameters belonging to F B 1 and omit that for F B 2 as it is similar.

Take parameter values ( b , c , e , p , δ 1 ) arbitrarily from F B 1 . Then system (1.3) with these parameter values becomes

(3.8) x → x + δ 1 x 1 − x − b y − e c + x , y → y + δ 1 p y 1 − y x .

The map (3.8) has a positive fixed point E 2 ∗ , whose eigenvalues are λ 1 = − 1 and λ 2 = 3 + G δ 1 with ∣ λ 2 ∣ ≠ 1 by Theorem 2.5. We consider a perturbation of (3.8) as follows:

(3.9) x → x + ( δ 1 + δ ∗ ) x 1 − x − b y − e c + x , y → y + ( δ 1 + δ ∗ ) p y 1 − y x ,

where ∣ δ ∗ ∣ ≪ 1 , which is a small perturbation parameter.

Let u = x − x 2 ∗ and v = y − y 2 ∗ , which transform the fixed point E 2 ∗ of map (3.9) into the origin. Moreover, (3.9) is transformed into

(3.10) u v → c 1 u + c 2 v + c 11 u 2 + c 12 u v + c 111 u 3 + c 13 u δ ∗ + c 23 v δ ∗ + c 113 u 2 δ ∗ + c 123 u v δ ∗ + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ δ ∗ ∣ ) 4 ) d 1 u + d 2 v + d 11 u 2 + d 12 u v + d 22 v 2 + d 111 u 3 + d 112 u 2 v + d 122 u v 2 + d 13 u δ ∗ + d 23 v δ ∗ + d 113 u 2 δ ∗ + d 123 u v δ ∗ + d 223 v 2 δ ∗ + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ δ ∗ ∣ ) 4 ) ,

where

(3.11) c 1 = 1 + δ 1 x 2 ∗ e ( c + x 2 ∗ ) 2 − 1 , c 2 = − δ 1 b x 2 ∗ , c 11 = δ 1 e c ( c + x 2 ∗ ) 3 − 1 , c 12 = − δ 1 b , c 111 = − δ 1 e c ( c + x 2 ∗ ) 4 , c 13 = x 2 ∗ e ( c + x 2 ∗ ) 2 − 1 , c 23 = − b x 2 ∗ , c 113 = e c ( c + x 2 ∗ ) 3 − 1 , c 123 = − b , d 1 = δ 1 p , d 2 = 1 − δ 1 p , d 11 = d 22 = − δ 1 p x 2 ∗ , d 12 = 2 δ 1 p x 2 ∗ , d 111 = d 122 = δ 1 p ( x 2 ∗ ) 2 , d 112 = − 2 δ 1 p ( x 2 ∗ ) 2 , d 13 = − d 23 = p , d 113 = d 223 = − p x 2 ∗ , d 123 = 2 p x 2 ∗ .

With the transformation

u v = T 2 x ˜ y ˜ , where T 2 = c 2 c 2 − 1 − c 1 λ 2 − c 1 ,

the map (3.10) becomes

(3.12) x ˜ y ˜ → − 1 0 0 λ 2 x ˜ y ˜ + f ( x ˜ , y ˜ , δ ∗ ) g ( x ˜ , y ˜ , δ ∗ ) ,

where

f ( x ˜ , y ˜ , δ ∗ ) = [ c 11 ( λ 2 − c 1 ) − c 2 d 11 ] c 2 ( λ 2 + 1 ) u 2 + [ c 12 ( λ 2 − c 1 ) − c 2 d 12 ] c 2 ( λ 2 + 1 ) u v − d 22 λ 2 + 1 v 2 + [ c 13 ( λ 2 − c 1 ) − c 2 d 13 ] c 2 ( λ 2 + 1 ) u δ ∗ + [ c 23 ( λ 2 − c 1 ) − c 2 d 23 ] c 2 ( λ 2 + 1 ) v δ ∗ + [ c 111 ( λ 2 − c 1 ) − c 2 d 111 ] c 2 ( λ 2 + 1 ) u 3 − d 112 λ 2 + 1 u 2 v − d 122 λ 2 + 1 u v 2 + [ c 113 ( λ 2 − c 1 ) − c 2 d 113 ] c 2 ( λ 2 + 1 ) u 2 δ ∗ + [ c 123 ( λ 2 − c 1 ) − c 2 d 123 ] c 2 ( λ 2 + 1 ) u v δ ∗ − d 223 λ 2 + 1 v 2 δ ∗ + O ( ( ∣ u ∣ + ∣ v ∣ + ∣ δ ∗ ∣ ) 4 ) , g ( x ˜ , y ˜ , δ ∗ ) = [ c 11 ( 1 + c 1 ) + c 2 d 11 ] c 2 ( λ 2 + 1 ) u 2 + [ c 12 ( 1 + c 1 ) + c 2 d 12 ] c 2 ( λ 2 + 1 ) u v + d 22 λ 2 + 1 v 2 + [ c 13 ( 1 + c 1 ) + c 2 d 13 ] c 2 ( λ 2 + 1 ) u δ ∗ + [ c 23 ( 1 + c 1 ) + c 2 d 23 ] c 2 ( λ 2 + 1 ) v δ ∗ + [ c 111 ( 1 + c 1