DYNAMICS OF A FOOD CHAIN MODEL WITH RATIO-DEPENDENT AND MODIFIED LESLIE-GOWER FUNCTIONAL RESPONSES

. The paper is concerned with a diﬀusive food chain model subject to homogeneous Robin boundary conditions, which models the trophic interac- tions of three levels. Using the ﬁxed point index theory, we obtain the existence and uniqueness for coexistence states. Moreover, the existence of the global attractor and the extinction for the time-dependent model are established under certain assumptions. Some numerical simulations are done to complement the analytical results.


1.
Introduction. Understanding the spatial and temporal behaviors of interaction species has become a central issue in population ecology, where one aspect of great interest for a multi-species model is whether the involved species can persist or have a coexistence steady state. At the beginning of 20th century, it is indicated by a constant positive solution of a ordinary differential equation system in the case where the species are homogeneously distributed; in recent years, stationary pattern induced by diffusion has been studied extensively, and many important phenomena have been observed.
In population dynamics, the traditional Lotka-Voltera type predator-prey systems are very important in the models of multi-species population dynamics. Yet the traditional prey-dependent predator-prey models have been challenged by several biological and physiological evidence that in many situations, especially, predators have to search, share or compete for food. A more suitable general predatorprey model appeared based on the so-called ratio-dependent theory, which can be roughly stated as that the per capita predator growth rate should be a function of the ratio of prey to predator abundance, and so there appears the so-called ratio-dependent function response. These hypotheses are strongly supported by numerous fields and laboratory experiments and observations. A general form of a ratio-dependent model is: . There are many authors who have studied the above system based on various response functions. In [2], the permanence and global stability of positive equilibria were studied for a delayed ratio-dependent predator-prey system; in [16,17], the long time behavior and stability of equilibria and stationary pattern formation were investigated to a diffusive one-prey and two-competing-predator system; in [19], the authors showed the global asymptotic stability of possible three steady states and the nonexistence of nontrivial positive periodic solutions. More details, see [13,6,10] and references therein.
More realistic models are derived in view of laboratory experiments and observations, the "carrying capacity" of the predator's environment is proportional to the population size of the prey, Leslie [20,21] introduced the following predator-prey model (also see [22,18,4,5]): . Furthermore, many kinds of modified Leslie-Gower models have been discussed, such as [28,14,12,15,11,1] and references therein. For example, Ji and Jiang proved the existence of the stationary distribution for a predator-prey model with modified Leslie-Gower and Holling-type II response functions in [14]; Aziz-Alaoui and Daher Okiye proved the boundedness of solutions and global asymptotic stability of the interior equilibrium for a predator-prey model with modified Leslie-Gower and Holling-type II response function in [1].
Recently, Upadhyay et al [29] have proposed the following model: . For the detailed background of the above model, one can see [29], where the boundedness of the system, existence of an attracting set, local and global stability of non-negative equilibrium points are established there. Particularly, we point out that the square term cZ 2 signifies the fact that mating frequency is directly proportional to the number of males as well as that of females present at any instant of time t [30].
In this paper, taking account of the ratio-dependent theory, we consider the following diffusive food chain system: subject to the boundary conditions and the initial conditions In this model, Ω denotes a bounded domain in R N (N ≥ 1) with smooth boundary ∂Ω, n is the outward unit normal over ∂Ω, k i > 0 (i = 1, 2, 3) are constants and for simplicity, we may take k i = 1 (i = 1, 2, 3) throughout the whole paper. A prey population of size u is the only food for the predator population of size v, which, in turn, is the food for the predator population of size w. a 1 is the intrinsic growth rate of the prey population u; a 2 is the intrinsic death rate of the predator v in the absence of the only food u; c measures the rate of self-reproduction of predator w. k is the carrying capacity, D and D 1 quantify the extent to which environment provides protection to the prey u and can be regarded as a measure of the effectiveness of the prey in evading a predator's attack. D 2 is the value at which per capita removal rate of v becomes β 1 /2, α 1 , α 2 , β 1 , β 2 are the maximum values which per capita growth rate can attain. And, we always assume α 2 > a 2 throughout the whole paper, so that the predator population of size v can persist. The interaction between the prey and the predator v is the ratio-dependent Holling type II functional response q(u, v) = uv u+Dv . The interaction between the predator v and the predator w is the Beddington-DeAngelis type functional response q(v, w) = vw D2+v+bw and the modified Leslie-Gower type functional response q(v, w) = w 2 v+D3 , where D 3 is added in the denominator to normalize the residual reduction in predator w due to severe scarcity of the predator v.
The present paper is built up as follows. In section 2, we give some known results about the eigenvalue problem and the fixed point index on positive cone; in section 3, we get the sufficient and necessary condition for the existence of positive steady states by using monotone method and the fixed point index theory on a positive cone, furthermore, the uniqueness is shown. The main result in this section is given in Theorems 3.2, 3.12 and 3.14; section 4 is devoted to consider global attractor and the extinction for the model; we illustrate our results with numerical simulations in section 5.
Consider the following equation and its steady state equation, the elliptic equation where Ω is a bounded domain in R N with smooth boundary ∂Ω. Assume that function f (x, u) : Ω × [0, ∞) → R satisfies the following hypotheses: (iii) If λ 1 (−f (x, 0)) < 0, then (2.2) has a unique positive solution which is globally asymptotically stable. In this case, the trivial solution u(x) = 0 is unstable.
Let E be a real Banach space and P ⊂ E be the natural positive cone of E. For y ∈ P , define P y = {x ∈ E : y + rx ∈ P f or some r > 0} and S y = {x ∈ P y : −x ∈ P y }. Then P y is a wedge containing P, y and −y, while S y is a closed subset of E containing y. Let T be a compact linear operator on E, which satisfies T (P y ) ⊂ P y . We say that T has property α on P y if there is a t ∈ (0, 1) and an element y 0 ∈ P y \S y such that (I − tT )y 0 ∈ S y . Assume A : P → P is a compact operator with a fixed point y ∈ P. Let L = A (y) be the Fréchet derivative of A at y. Then L maps P y into itself. Dancer's theorem can be stated as follows.  Suppose E 1 and E 2 are ordered Banach spaces with positive cones C 1 and C 2 , respectively. Let E = E 1 × E 2 and C = C 1 × C 2 . Then clearly E is an ordered Banach space with positive cone C. Let Γ be an open set in C containing 0 and A i : Γ → C i be completely continuous operators, i = 1, 2. Denote by (u, v) a general element in C with u ∈ C 1 and v ∈ C 2 . Let A : Γ → C be defined by We also define Suppose A 2 : Γ → C 2 extends to a continuously differentiable mapping of a neighborhood of Γ into E 2 , C 2 − C 2 is dense in E 2 and T = {u ∈ U : u = A 1 (u, 0)}.Then the following conclusions are true: (i) deg C (I − A, U × C 2 (δ), 0) = 0 for δ > 0 small if for any u ∈ T , the spectral radius r(A 2 (u, 0)| C2 ) > 1 and 1 is not an eigenvalue of A 2 (u, 0)| C2 corresponding to a positive eigenvector; 3. Existence and uniqueness. Thanks to Lemma 2.2, we give the following notations. Suppose a, m, d are positive constants, A(x) ∈ C α (Ω) (0 < α < 1) and A(x) > 0 on Ω. If a > λ 1 , we denote by θ[a, A]the unique positive solution of (3.1). If m − d > λ 1 , we denote by ζ[m, d, A] the unique positive solution of (3.2).
Theorem 3.1. Assume α 2 > a 2 holds. Then any nonnegative solution (u, v, w) of (1.4)-(1.5) has a priori bounds: The proof is standard. Since u satisfies it follows from the maximum principle and the Hopf lemma that u ≤ k .
= M 1 . The other two inequalities can be proved similarly. Now we give the sufficient and necessary conditions for the existence of positive solution to system (1.4)-(1.5) by means of the upper-lower solution method.
By Lemma 2.2 and the upper-lower solution method for elliptic equation, we con- ).
v+D3 . Then (f, g, h) has a mixed quasimonotone form. It is easy to check that (u * , v * , w * ), (u * , v * , w * ) are the couple upper and lower solutions for (1.4)-(1.5), i.e., u * ≥ u * , v * ≥ v * , w * ≥ w * on Ω, and they satisfy the following inequalities: Evidently Therefore, the existence of positive solution follows from the existence-comparison argument [24] for the elliptic system with a mixed quasimonotone form. This finishes the proof of the theorem. The proof of Theorem 3.3 is fundamental, we omit it here.
In the rest of the section, we give the existence and uniqueness of positive solution by using the fixed point index theory on a positive cone. It is obvious that (0, 0, 0) is the singular point of (1.4)-(1.5), we follow the idea of Kuang and Berreta in [19] and assume f (0, 0, ·) = 0 and g(0, (ii) Nonnegative solutions with exactly one component identically zero It is easy to see that all the nonnegative solutions of (1.4)-(1.5) are in Ξ by and also an operator A by By the definition of M , A is a positive operator, which maps Ξ to P . It follows from the standard elliptic regularity theory that A is a completely continuous operator.

Proof. (i) Define an operator
where θ ∈ [0, 1], A θ is also a positive, completely continuous operator. If (u, v, w) is a nonnegative fixed point of A θ , then one can show that (u, v, w) must satisfy (3.3). Hence, A θ has no fixed point on ∂Ξ. It then follows from the homotopy invariance of the fixed point index that index P (A θ , Ξ) ≡ constant. We see that A 0 has a unique fixed point (0, 0, 0) in Ξ and index P (A 0 , (0, 0, 0)) = 1 by Lemma 2.3.
, and also an operator A by Since the second part of the proof is long, we divide it into two propositions. Proposition 1. If a 1 > λ 1 and c = λ 1 or c > λ 1 and a 1 = λ 1 , then for all > 0, index P (A , (0, 0, 0)) = 0.
Proposition 2. If a 1 − α1 D > λ 1 and c = λ 1 , then index P (A, (0, 0, 0)) = 0. Proof. From Proposition 1, we know (0, 0, 0) is an isolated fixed point of A in P for all > 0 if a 1 − α1 D > λ 1 . We claim that (0, 0, 0) is an isolated fixed point of A in P if a 1 − α1 D > λ 1 . Suppose that the statement was not true. Then there nonnegative nontrivial solution U n = (u n , v n , w n ) such that U n → (0, 0, 0) as n → ∞ and (I − A)U n = 0. Then U n must be one of the following three cases for large n: U n with w n ≡ 0; U n with u n ≡ 0; U n with v n ≡ 0.
For the case U n with w n ≡ 0, by the maximum principle, w n > 0. Let w n = w n /||w n || ∞ . Then w n satisfies Passing to a subsequence if necessary, by L p estimates and Sobolev embedding theorems, there exists w, such that w n → w in C 1 as n → ∞. Taking the limit of the equation of w n , we obtain −∆w = cw, then c = λ 1 , a contradiction. Hence, w n ≡ 0 and then u n , v n satisfy un+D1vn . It is obvious that v n ≡ 0 if u n ≡ 0; u n = θ[a 1 , a1 k ] if v n ≡ 0. So, we only need to consider the case U n = (u n , v n , 0) with u n ≥, ≡ 0, v n ≥, ≡ 0. Due to the maximum principle, it is sufficient to consider the case U n = (u n , v n , 0) with u n > 0, v n > 0.
k ] is defined in Theorem 3.2. So, u n → 0 can not hold forever, a contradiction with U n → (0, 0, 0) as n → ∞. The claim has been proved.
It is obvious that t 0 > 1. On the contrary, we suppose that L has property α on P (0,0,θ[c, β 2 Now we assume that µ is an eigenvalue of L which are greater than 1. Then for all (ϕ, ψ, η) ∈ E, we have L(ϕ, ψ, η) = µ(ϕ, ψ, η), i.e., in Ω,  We now turn to consider the following two sub-systems, which are very important for the existence of semi-trivial solutions.
Now, we consider the existence of positive solution of (3.4).
are defined in Theorem 3.1. Define an operator A : E → E by By the denotation of M , A is a positive operator, which maps Ξ to P . Similarly to the proof before, we only give Lemmas 3.8 and 3.9 and omit the proof here. The proof of Lemma 3.10, which uses the same arguments as those of Theorem 3.12, will be proved later.
According to Theorem 3.1, it follows that u, v, w are uniformly bounded on Ω independent of α 1 and β 1 . By the standard regularity theory, one can get the result of Theorem 3.13.
5. Numerical simulation. The goal of this section is to present the results of numerical simulations which complement the analytic results in section 3. We simulate the corresponding system of (1.1)-(1.3) with k i = 1 (i = 1, 2, 3) in the one-dimensional space domain. Without loss of generality, we take Ω = (0, l), l > 0. Therefore, (1.1)-(1.3) become We perform the initial-boundary-value problem numerically based on the Crank-Nicholson scheme and let l be 10. In each simulation, the figures are plotted at sufficiently final time (here, we take t = 100), which allow us to regard the solutions as steady states. In the finite difference scheme, we take the temporal axis with grid spacing t = 0.01 and the spatial axis with grid spacing x = 0.1. In addition, √ λ 1 satisfies cot(10s) = 1 2 s − 1 2s , 0 < s < π/10. Through MATLAB software, we can easily get the approximate number λ 1 ≈ 0.0691. Our numerical simulations indicate the following outcomes.
(i) In Figures 1-2 Here, we consider the weaker condition a 1 > λ 1 replaced the one in the previous hypotheses, under which the positive steady states maybe exist numerically, and proceed with the similar argument in outcome (ii). We know that if α 2 − a 2 > λ 1 + (β 1 cD 3 )/(D 2 β 2 + bcD 3 ), which is denoted by (H), then α 2 − a 2 > λ 1 ( β1θ[c, (ii) In Figure 3, we are concerned with the effects of the parameter α 1 , β 1 on the density of positive steady state solutions. Here, we fix k = 1, α 2 = 1.0, a 1 = 1.0, a 2 = 0.1, β 2 = 1.0, D = 1.0, D 1 = 1.0, D 2 = 2.0, D 3 = 3.0, b = 1.0, c = 0.15, and fix α 1 , let β 1 change or fix β 1 , let α 1 change. Figure 3 (a)-(d) present a phenomenon that the population density of u is increasing as β 1 increases, but the population density of v and w is in the opposite direction; Figure 3 (e)-(h) present a phenomenon that the population density of u, v and w decreases with increasing of the parameter α 1 . In addition, for some parameter region we have chosen, there appears a small bump in the distribution of predator population of size v near the boundary, but the prey population of size u and the predation population of size w distribute by a tendency of "a high and flat central portion, but two low end portions". Remark 1. In the above simulations, we consider the weaker condition a 1 > λ 1 replaced the one in the hypotheses of Theorem 3.12, under which the positive steady