Qualitative analysis of a modified Leslie-Gower predator-prey model with Crowley-Martin functional responses

In this paper, we study a modified Leslie-Gower predator-prey model with Crowley-Martin functional response. We show the existence of a bounded positive invariant attracting set and establish the permanence conditions. The parameter regions for the stability and instability of the unique constant steady state solution are derived, and the existence of time-periodic orbits and non-constant steady state solutions are proved by bifurcation method.

1. Introduction. In this paper, we consider the following reaction-diffusion equations: where Ω ⊂ R N , N ≥ 1, is a bounded domain with smooth boundary ∂Ω, u and v represent the population densities at time t and location x; a 1 , a 2 , b, c, d, e, f and K are model parameters assuming only positive values and defined as follows: • a 1 and K are the intrinsic growth rate and the carrying capacity of prey population u respectively. • bv/((1+cu)(1+dv)) is the Crowley-Martin functional responses term, in which c measures the magnitude of interference among prey (see [6]). • a 2 describes the growth rate of predator v.
• ev/(u + f ) is the modified Leslie-Gower term, in which e is the maximum value which per capita reduction rate of v can attain, f measure the extent to which environment provides protection to predator v (see [2]).
• ν is the outward unit normal vector of the boundary ∂Ω and the homogeneous Neumann boundary condition "∂ ν u = ∂ ν v = 0 on ∂Ω"indicates that the system is self-contained with zero population flux across the boundary (see [5]). • The positive constants d 1 and d 2 are the diffusion coefficients, and the initial data u 0 (x) and v 0 (x) are nonnegative nontrivial continuous functions due to its biological sense. Next we introduce some background of the model (1) and we assume all the parameters in the models listed below are positive. In [11], Leslie and Gower proposed the following classical predator-prey model Several ecological treatises regarded (2) as a prototypical predator-prey system (see [13,16]). The term ev/u is called the Leslie-Gower term. It measures the loss in the predator population due to rarity of its favorite food. In the case of severe scarcity, v can switch over to other population but its growth will be limited by the fact that its most favorite food (u) is not available in abundance (see [2]). By considering the above reasons, and using Holling type-II functional response (see [10]) in both prey and predator interaction terms, a modified Leslie-Gower predator-prey system was proposed in [2]: The boundedness of solutions, existence of an attracting set and global stability of the coexisting interior equilibrium of (3) are studied in [2]. When the spatial dispersal is also considered, system (3) becomes Hopf and Turing bifurcations of (4) with no-flux boundary condition were analyzed in [3,4,19] and numercial simulations showed rich pattern formation dynamics and self-organization of various patterns. Positive steady state solutions of (4) and its related problems with d 1 = d 2 = 1 and zero Dirichlet boundary condition were considered in [15,26,27,28,29], and the existence, multiplicity and uniqueness of positive solutions were given in these papers. The Holling type II functional response bu/(1 + cu) in (3) is classified as one of prey-dependent functional response, it assumes that predators do not interfere with one another's activities; thus competition among predators for food occurs only via the depletion of prey. However, Crowley and Martin proposed a functional that can accommodate interference among predators (see [6]). This type of functional response is classified as one of predator-dependent functional response, i.e. that are functions of both prey and predator abundance because of predator interference. It is assumed that predator-feeding rate decreases by higher predator density even when prey density is high, and therefore the effects of predator interference on feeding rate remain important all the time whether an individual predator is handing or searching for a prey at a given instant of time. The per capita feeding rate in this formulation is given by mu/(1 + Au + Bv + ABuv), where m, A and B describe the effects of capture rate handling time and the magnitude of interference among predators, respectively, on the feeding rate (see [1,7,12,18,21,22] for related works).
By taking above considerations, Ali and Jazar [1] consider a predator-prey model which incorporates a modified version of the Leslie-Gower with Crowley-Martin functional response: which is the ODE system corresponding to (1). In [1], the asymptotic behavior of the positive equilibrium and the existence of Hopf-bifurcation of nonconstant periodic solutions surrounding the interior equilibrium are considered.
On the other hand, the spatial component of ecological interactions has been identified as an important factor in how ecological communities are shaped, and understanding the role of space in challenging both theoretically and empirically [14]. Empirical evidence suggests that the spatial scale and structure of environment can influence population interactions [5]. By considering the above reasons, we consider the problem (1). Problem (1) is rendered dimensionless using the following change of variables and parameters: With the scaling, system (1) takes the form where In this paper we focus on the question of existence and stability of steady state solutions and periodic orbits of (6). The steady state equations associated with (6) are It is obvious that (7) always have three boundary equilibria equilibria E 0 (0, 0), E 1 (1, 0) and E 2 (0, E/D). By [1], (6) has a unique positive equilibrium and any of the following three conditions is satisfied Furthermore, (u * , v * ) is given by where In the rest of this article, we always assume the following condition hold, which implies (6) has a unique positive equilibrium E * = (u * , v * ).
In this paper, we mainly discuss for what parameter ranges, time-periodic patterns and non-constant stationary patterns are possible. These patterns have been predicted by Turing [20] for a wide class of reaction-diffusion models. The organization of the remaining part of the paper is as follows. In Section 2 we study the permanence conditions. In Section 3 we study the stability and Hopf bifurcation results to the ODE model corresponding to (6). In Section 4 we analyze the local and global stability of the uniform steady state (u * , v * ) of the PDE model (6). In Sections 5 and 6 we use bifurcation theory to prove the existence of periodic and non-constant solutions. Some numerical simulations of periodic orbits and non-constant steady state solutions are also shown at the end of Section 5 and Section 6, respectively. Throughout this paper, N is the set of natural numbers and N 0 = N ∪ {0}. The eigenvalues of operator −∆ with homogeneous Neumann boundary condition in Ω are denoted by 0 = µ 0 < µ 1 ≤ µ 2 ≤ · · · ≤ µ n ≤ · · · , and the eigenfunction corresponding to µ n is denoted by φ n (x).

2.
Permanence. In this section, we present the permanence of the system (6), that is, uniform persistence plus dissipation. Firstly, we give the definition of permanence (see for example [5]). Definition 2.1. System (6) is said to be permanent if there are positive constants M 1 and M 2 such that the solution (u(x, t), v(x, t)) of (6) with nonnegative and nontrivial initial datas u 0 and v 0 satisfies In order to get the permanence, we first give a dissipative result to the solution of (6).
The proof of Theorem 2.2 depends on the asymptotic behavior of the following logistic equation: where d, r and K are positive constants. The following result was obtained in [24]: Assume w 0 (x) be a nonnegative, nontrivial and continuous function. The solution w(x, t) of (13) satisfies Proof of Theorem 2.2. We only give the proof of Part 2 since Part 1 follows easily from maximum principle of parabolic equations (see [8]). From the first equation of (6), we have u t ≤ D 1 ∆u + u(1 − u), from the comparison principle of parabolic equations and Lemma 2.3, for any ε > 0, there exists a positive constant t 1 such that which implies lim sup It follows from the second equation of (6) and (14) that Hence there exists a constant and the arbitrariness of ε implies the second inequality in Part 2 holds. Now, we show that system (6) is persistent. From the viewpoint of biology, this implies that the two species of prey and predator will always coexist at any time and any location of the inhabited domain, no matter what their diffusion coefficients are (see [5]).
Theorem 2.4. Assume u 0 and v 0 are nonnegative and nontrivial continuous functions. Let (u(x, t), v(x, t)) be the solution of (6). If B ≥ 1 or it is easy to see (8) holds under the assumptions in Theorem 2.4.
Proof of Theorem 2.4. From the first equation of (6) and (15), we have where v := 1+ε+E D + ε. Since (16) holds, we can choose ε small enough such that Hence there exists a constant t 3 ≥ t 2 such that It follows from the second equation of (6) and (18) that for ε > 0 small enough. From (18) and (19), we obtain (17).
3. Stability and Hopf bifurcation to the ODE model. In this part, we collect some stability and bifurcation results to the ODE model corresponding to (6), the ODE model is given by Obviously, (u * , v * ) giving in (12) is the unique positive equilibrium of (20) under the condition (H). In the following we fixed parameters A, B, D, E > 0 (which implies (u * , v * ) is fixed), and use C as the main bifurcation parameter. The Jacobian matrix of system (20) at (u * , v * ) is where The characteristic equation of L 0 (C) is where The following results hold (see [1]) with the assumption (H): 4. Stability with respect to the PDE model. In this part, we consider the local and global stability of the positive equilibrium E * with respect to the PDE model (6). Linearizing the system (6) about the equilibrium E * , we obtain an eigenvalue problem x ∈ Ω, where C * is defined in (22). Denote For each n ∈ N 0 , we define a 2 × 2 matrix The following statements hold true by using Fourier decomposition: 1. If µ is an eigenvalue of (25), then there exists n ∈ N 0 such that µ is an eigenvalue of L n (C). 2. The constant equilibrium (u * , v * ) is locally asymptotically stable with respect to (6) if and only if for every n ∈ N 0 , all eigenvalues of L n (C) have negative real part. 3. The constant equilibrium (u * , v * ) is unstable with respect to (6) if there exists an n ∈ N 0 such that L n (C) has at least one eigenvalue with nonnegative real part. The characteristic equation of L n (C) is where and In the rest of this article, we always assume that R(u * , v * ) > 0 to study the stability and bifurcation results.
Next we consider the global stability of E * by using two methods. Firstly, we construct a Lyapunov function and establish the global stability. Secondly, we get the global stability by using iterative method.
By using Lyapunov function, we obtain the following result: Theorem 4.3. Assume (H) holds and E * is locally asymptotically stable. E * is globally asymptotically stable if and any one of the following conditions holds: and It can be easily verified that the function V (u, v) is zero at the equilibrium E * and it is positive for all other positive values of u and v, and thus, E * is the global minimum of V . Using the relation between u * and v * , a direct calculation implies (see [1, Theorem 6.3]) Then E * is globally asymptotically stable.
Next we study the global stability of E * by using iterative method. To this end we define a function F(α, β) by We claim {u n } ∞ n=1 , {u n } ∞ n=1 , {v n } ∞ n=1 and {v n } ∞ n=1 satisfy (a): u 1 < u 2 < · · · u n < u n < · · · < u 2 < u 1 ; Proof of the claim. By the monotone property of Φ and Ψ, (a) and (b) hold if and only if u 1 < u 2 < u 2 < u 1 .
After a series of calculations, one can get So u 1 < u 2 and v 1 < v 1 . Furthermore, by the definition of Φ, we have u 2 < 1 = u 1 . Finally, u 2 < u 2 follows from u 1 < u 1 , v 1 < v 1 and the monotone property of Φ. So, (45) holds, and (a), (b) hold. Next we will prove (c) and (d) by induction. Let (u(x, t), v(x, t)) be the solution of (6). It follows from Theorems 2.2 and 2.4 that So (c) and (d) hold for n = 1. Next, we assume (c) and (d) hold for n = k, i.e., and prove (c) and (d) hold for n = k + 1. For a fixed positive constant ε small enough, it follows from (46) that there exists a positive constant T 1 1 such that By (47) we have for x ∈ Ω and t ≥ T 1 It follows from Lemma 2.3 and comparison principle for parabolic equations that By letting ε → 0, we obtain (c) and (d) hold for n = k + 1.
The above analysis shows that the claim holds. From the the monotonicity of the sequence, we may assume where u, u, udv and v are four positive constants satisfy
Next we give a corollary from Theorem 4.4.

Hopf bifurcation.
In this part, we study the existence of periodic solutions of (6) by analyze the Hopf bifurcations from the positive equilibrium E * , and we will show the existence of spatially homogeneous and spatially non-homogeneous periodic orbits. In this section and also in section 6, we assume that all eigenvalues µ i are simple, and denote the corresponding eigenfunction by φ i (x), where i ∈ N 0 . Note that this assumption always holds when N = 1 for Ω = (0, π), as for i ∈ N 0 , µ i = i 2 / 2 and φ i (x) = cos(ix/ ), where is a positive constant; and it also holds for generic class of domain in higher dimensions. To identify possible Hopf bifurcation value C H , we recall the following necessary and sufficient condition from [9,25]. For i ∈ N 0 , we define where C H (µ) is defined in (32). So T i (C H i ) = 0 and T j (C H i ) = 0 for all j ∈ N 0 \ {i}. By Lemma 4.1, it is easy to see that C H i is strictly decreasing in i and max i∈N0 C H i = C * and lim where C * is defined in (22). Then there exists a unique n 0 ∈ N 0 such that µ n0 < µ * # ≤ µ n0+1 and C H n0 > C * # ≥ C H n0+1 , where C * # = C H (µ * # ) and µ * # is defined in (39). We have n 0 + 1 possible Hopf bifurcation points at C = C H j (0 ≤ j ≤ n 0 ) defined by (55), and these points satisfy Furthermore, it follows from Lemma 4.1 (see also Fig. 1) that C H i > C S (µ i ), i.e. D i (C H i ) > 0 for all i = 0, 1, 2, · · · , n 0 . Finally let the eigenvalues close to the pure imaginary one near at C = C H i , i = 0, 1, 2, · · · , n 0 , be α(C) ± iω(C). So Now by using the Hopf bifurcation theorem in [25], we have Theorem 5.1. Suppose A, B, D, E > 0 are fixed such that C * > 0 and R(u * , v * ) > 0, D 1 and D 2 are fixed positive constants, where C * is defined in (22) and R(u * , v * ) is defined in (24). Let Ω be a bounded smooth domain so that the spectral set S = {µ i } i∈N0 satisfy that: (S 1 ): all eigenvalues µ i are simple for i ∈ N 0 . Then there exists unique n 0 ∈ N such that µ n0 < µ * # ≤ µ n0+1 and C H n0 > C * # ≥ C H n0+1 , where C * # = C H (µ * # ) and µ * # is defined in (39), and for (6) there exist n 0 +1 Hopf bifurcation points C H j , j = 0, 1, 2, · · · , n 0 , defined by (55), satisfying (56). At each C = C H j , the system (6) undergoes a Hopf bifurcation, and the bifurcation periodic orbit near (C, u, v) = (C H j , u * , v * ) can be parameterized as (C(s), u(s), v(s)), so that C(s) ∈ C ∞ in the form of C(s) = C H j + o(s) for s ∈ (0, δ) and δ is a small positive constant, and where ω(C H j ) = D j (C H j ) is the corresponding time frequency, φ j (x) is the corresponding spatial eigenfunction, and (l j , m j ) s the corresponding eigenvector, i.e. (26) and I is the identity map. Moreover, (1): the bifurcation periodic orbits from C = C H 0 = C * are spatially homogeneous, which coincide with the periodic orbits of the corresponding ODE system; (2): the bifurcation periodic orbits from C = C H j , 1 ≤ j ≤ n 0 , are spatially non-homogeneous.  Fig. 2 shows a numerical simulation with parameters chosen above, and the solution converges to a spatially homogeneous periodic orbit. 6. Steady state bifurcation. In this part, we analyze the properties of steady state solution bifurcation for (7). The following a priori estimate can be easily established via maximum principle, and we omit the proof. and Recall T n (C) and D n (C) defined in (29) and (30). Now we identify steady state bifurcation value C S of (7), which satisfies the following steady state bifurcation condition [9,25]: (A S ): there exists n ∈ N 0 such that D n (C S ) = 0, T n (C S ) = 0, D j (C S ) = 0 and T j (C S ) = 0 for j ∈ N 0 \ {n}, and D n (C S )) = 0.
Apparently, D 0 (C) = Q(u * , v * )C = 0 since we assumed C > 0 and Q(u * , v * ) > 0, where Q(u * , v * ) is defined in (31). Hence we only consider n ∈ N. In the following we fixed arbitrary A, B, D, D 1 , D 2 , E > 0, and determine C values satisfy condition (A S ). We notice that D n (λ) = 0 is equivalent to C = C S (µ), where C S (µ) is defined in (32). Here we make the following additional assumption on the spectral set S = {µ i } i∈N0 according to Lemma 4.1: (S 2 ): there exits p, q ∈ N such that µ p−1 < µ * # < µ p ≤ µ q < C * /D 1 ≤ µ q+1 , where µ * # and C * are defined in (39) and (22), respectively. In the following, we denote The points C S n defined above are potential steady state bifurcation points. It following from Lemma 4.1 that for each n ∈ p, q , there exists only one point C S n such that D n (C S n ) = 0 and T n (C S n ) = 0. One the other hand, it is possible for somẽ C ∈ (C * # , C # ) with C * # and C # defined in Lemma 4.1 and some i, j ∈ p, q with i < j such that C S i = C S j =C, where and Q(u * , v * ) is defined in (31). Then for C =C, (A S ) is not satisfied, and we shall not consider bifurcations at such point. On the other hand, it is possible such that However, from an argument in [25], for N = 1 and Ω = (0, π), there are only countably many , such that (61) or (62) occurs for some i = j. For general bounded domains in R N , one can also show above argument also hold (see [23]). According to above analysis, to satisfy the bifurcation condition (A S ) we only need to verify D n (C S n ) = 0. In fact, we have D n (C S n ) = D 1 µ n + Q(u * , v * ) > Q(u * , v * ) > 0. Summarizing the above discussion and using a general bifurcation theorem [25], we obtain the main result of this part on bifurcation of steady state solutions: Theorem 6.2. Suppose A, B, D, E > 0 are fixed such that C * > 0 and R(u * , v * ) > 0, D 1 and D 2 are fixed positive constants, where C * is defined in (22) and R(u * , v * ) is defined in (24). Let Ω be a bounded smooth domain so that the spectral set S = {µ i } i∈N0 satisfying (S 1 ) and (S 2 ). Then for any n ∈ p, q , which is defined in (60), there exists unique C S n ∈ (0, C # ], where C # is defined in (36) such that D n (C S n ) = 0 and T n (C S n ) = 0. If in addition, we assume (S 3 ): C S j = C S j for any i, j ∈ p, q and i = j, and C S k = C H l for any k, l ∈ p, q and k = l, where C H l is defined in (55). 1. There is a smooth curve Γ n of positive solutions of (7) bifurcating from (C, u, v) = (C S n , u * , v * ); near (C, u, v) = (C S n , u * , v * ), Γ n = {(C n (s), u n (s), v n (s)) : s ∈ (− , )}, where u n (s) = u * + sl n φ n (x) + sψ 1,n (s) v n (s) = v * + sm n φ n (s) + sψ 2,n (s), for some C ∞ functions C n , ψ 1,n , ψ 2,n such that C n (0) = C S n and ψ 1,n (0) = ψ 2,n (0) = 0. Here (l n , m n ) satisfies L(C S n )[(l n , m n ) T φ n (x)] = (0, 0) T , where L(C) is defined in (26).
2. In addition, if we assume (16) holds, then Γ n contained in a global branch Σ n of positive nontrivial solutions of (7); and either Σ n contains another (C, u, v) = (C S m , u * , v * ) for m = n, or the projection proj C Σ n of Σ n on the C-axis satisfies proj C Σ n ⊃ (0, C S n ) or proj C Σ n ⊃ (C S n , ∞). Proof. The condition (A S ) has been proved in the previous paragraphs, and the bifurcation of solutions to (7) occur at C = C S n . Note that we assume (S 3 ) holds, so C = C S n is always a bifurcation from simple eigenvalue point. From the global bifurcation theorem in [17] and Lemma 6.1, Γ n is contained in a global branch Σ n of nontrivial positive solutions. Furthermore, either Σ n contains another (C, u, v) = (C S m , u * , v * ) for m = n or Σ n intersects the boundary C = 0 or Σ n is unbounded, and proj C Σ n ⊃ (0, C S n ) if Σ n intersects the boundary C = 0, proj C Σ n ⊃ (C S n , ∞) if Σ n is unbounded.