1 Introduction

Many ecological phenomena among different populations can be characterized or simulated by various mathematical models. By analyzing different kinds of mathematical models people may give scientific predictions and explanations on the dynamics of these models, and further, put forward reasonable projects corresponding to some ecological problems. In many kinds of ecological-mathematical models, the predator–prey model is a very important branch. Since the classical Lotka–Volterra ecological model [13] was introduced into investigations, the predator–prey-type models received extensively attentions. In recent decades, different kinds of predator–prey models have been established, and many great progress and valuable achievements have been made; see [4,6,9,20,22,23] for examples.

In predator–prey models, the functional responses are very important terms in describing the relations between predator and prey, they determine many dynamical properties of the systems in some extent. Specifically, the functional response describes the change of predation rate of predator capturing prey, that is, how the predator consuming prey is dependent on the response terms. In large number of predator–prey models, the functional responses are of prey-dependent type, such as the classical Lotka–Volterra-type functional response, the Holling-type II (sometimes referred to as the Michaelis–Menten functional response), III and IV functional responses, the Ivlev and inverse Ivlev functional responses, etc. Taking these functions as the functional responses, a great many of predator–prey models have been studied thoroughly, and quite a lot of important research achievements have been given; see [1, 10, 18, 21, 24], etc.

The introduction of prey-dependent-type functional response is based on the assumption that the predation rate of predator is only effected by prey, and the predators are mutual noninterference within themselves. In fact, generally speaking, the phenomenon in sharing or competing their foods always occurs during the process of predation. When the population density of predator is high, the predation rate will decrease. Moreover, as the well-known environmental paradox being proposed by Hairston, et al. [8] and Rosenzweig [19], the type of functional responses gave rise to many strong controversies among ecologists. Just in this situation, based on the Holling-type II functional response, by a great deal of experimental and numerical analysis, Beddington [3] and DeAngelis et al. [5] introduced a new type functional response, the predator–prey dependent functional response f (x) = bx/(a + x + cy), which is called the Beddington–DeAngelis-type functional response, where x and y represent the predator and prey, respectively, and a, b, c are parameters. This type of functional response is very like with the Hollingtype II functional response, the only difference is that the Beddington–DeAngelis-type functional response involves both prey and predator. Compared with the prey-dependenttype functional responses, the predator–prey-dependent functional response not only reflects the mutual interference with predators, but also remains the merits of the ratiodependent functional responses. Moreover, it avoids some controversies induced by low population density and better describes the predation effect of the predator versus prey. It is just because of these reasons, different kinds of models with predator–prey-dependent functional response have become a center topic in population dynamics.

Except the Beddington–DeAngelis-type functional response, another attentive predator–prey-dependent functional response hx/y is called the Leslie–Gower-type functional response, where h is a constant. It was introduced by Leslie [11] and discussed jointly by Leslie and Gower [12]. This functional response measures the self-consumption of predator when the foods are scarce, where the growth of predator population is of logistic form. Since in actual ecological systems, the carrying capacity set by environmental resources is proportional to prey abundance and the predators always try to survive by catching other preys when their conventional food is short seriously, in such situations, the growth rate of the predator would be affected. For this reason, Aziz-Alaoui [2] modified the logistic form and proposed the Leslie–Gower-type functional response hx/(a + y), where a is a positive constant and measures the environmental protection for predator. With such a functional response, Aziz-Alaoui et al. studied some predator–prey models in spatially homogeneous or inhomogeneous cases; see [7, 15,16,17, 25]. Based on and motivated by all the above-mentioned, in this paper, we introduce the following nondimensional reaction– diffusion system

with two modified predator–prey-dependent nonlinear functional responses αv/(c + u + mv), βv/(r + u). Here Ω is a bounded domain in Rn with smooth boundary ∂Ω, u and v are the densities of prey and predator populations. d1, d2 and a, b are the diffusion rates and the growth rates of u and v, respectively; α is the predation rate of predator, and c is the self-saturation density of prey; m represents the semi-saturation term of the system; β measures the conversion rate of v predating u, and r reflects the protection of the surroundings on predators; ν is the outward unit normal vector on ∂Ω. All parameters are positive constants. u0(x), v0(x) are booth smooth functions in Ω.

For system (1), we consider the large-time behavior of positive solutions. In addition, corresponding to (1), we mainly analyze the steady-state system

An outline of this paper is as follows: In Section 2, we consider the large-time behavior of positive solutions to system (1) mainly by the comparison principle of parabolic equations. In Section 3, we give a priori estimate for the positive solutions of (2) by the maximum principle. Then, in Section 4, we investigate the coexistence of predator and prey to system (2) by the energy integral techniques and topological degree computation. In Section 5, by taking the diffusion rate of predator as a parameter, we study the bifurcating phenomenon, which emits from the unique positive constant equilibrium also by using the topological degree techniques. Since pattern is a very interesting nonlinear phenomena and pattern dynamics is a primary branch of nonlinear science, finally, we give some numerical examples and depict the corresponding trajectory graphs or spatiotemporal patterns to simulate the related theoretical results in Section 6. It is worth mentioning that the numerical simulated graphs imply that the available food resource for predator may not single if the system parameters are controlled properly.

2 Large-time behavior

In this section, we consider the large-time behaviors of positive solutions to system (1) including the global attractor and the persistence.

Lemma 1. Let k > 0 be a constant. If the function f (x, t) satisfies

with Ω being the same as that of in (1) and f0 being a smooth function in Ω, then limt→∞ f (x, t) = k holds uniformly in Ω.

Proof. For any χ0 > 0, it is well known that the problem

has a unique solution χ = χ(t; χ0) satisfying limt→∞ χ(t; χ0) = k. Here we use the mark χ(t; χ0) to represent that the solution χ is related to χ0.

Let M = maxΩ f0(x). Then M > 0 and f (M, t) is an upper-solution of (3). Obviously, 0 is a lower-solution of (3). These lead to (3) permitting a unique nonnegative solution f (x, t), and f (x, t) satisfies 0 ≤ f (x, t) ≤ f (M, t). Then the maximum principle induces that f (x, t) > 0, x ∈ Ω, t > 0. Take some τ > 0. Then f (x, τ ) > 0 for x ∈ Ω. Denote z(x, t) = f (x, t + τ ). Then z(x, t) satisfies

By the comparison principle of parabolic equations, for any t ≥ 0, we have

with k = minΩ f (x, τ ), k = maxΩ f (x, τ ). Since limt→∞ f (k, t) = limt→∞ f (k, t) = k, then limt→∞ f (x, t + τ ) = k. Therefore, limt→∞ f (x, t) = k holds uniformly in Ω.

Theorem 1. If (u(x, t), v(x, t)) is a positive solution of (1), then

Hence, for any ϵ > 0, the rectangle [0, a + ϵ] × [0, b + ϵ] is a global attractor of (1).

Proof. Clearly, u(x, t) satisfies

By the comparison principle of parabolic equations and Lemma 1 we know that for any 0 < ε « 1, there exists t1 > 0 such that u(x, t) ≤ a + ε, x ϵ Ω when t > t1. Hence, we get lim supt →maxΩ u(x, t) ≤ a + ε.

Likewise, using the comparison principle of parabolic equations and Lemma 1 again, we also know that there is t2 with t2 > t1 such that v(x, t) ≤ b + ε, x Ω for 0 < ε 1 when t > t2. So, lim supt→∞maxΩ v(x, t) ≤ b + ε. Then the arbitrariness of ε induces that the conclusion holds.

The result shows that any positive solution (u(x, t), v(x, t)) of (1) lies in a bounded region when t → ∞, that is, (u(x, t), v(x, t)) exists globally, and the rectangle [0, a + ϵ] × [0, b + ϵ] is a global attractor of (1) for any ϵ > 0.

Theorem 2. If ac > bα, then (1) is persistent. Specifically, any positive solution (u(x, t), v(x, t)) of (1) satisfies

Proof. For any ε > 0, the proof of Theorem 1 implies that u(x, t) satisfies

where t2 is defined in the proof of Theorem 1. Similar as that we prove Theorem 1, there exists t3, t3 > t2 such that u(x, t) ≥ a −(α(b +ε))/c− ε, x ∈ Ω when t > t3. Therefore,

Since u(x, t) ≥ 0, x ∈ Ω, then v(x, t) satisfies

Hence, there is t4 with t4 > t3 such that v(x, t) ≥ br/(r + β) ε > 0, x ϵ Ω when

Then the arbitrariness of ε implies that

This shows that system (1) is persistent.

Remark 1. Theorem 2 implies that the predator and prey always coexist in their living surroundings despite the location and time as long as ac > bα, moreover, this coexistence is independent of their diffusion situation.

3 A priori estimate

In this section, we seek a priori estimate for the positive solutions of (2), which will be used frequently in the latter contents. As preparation, we introduce the following lemma due to [14].

Lemma 2. Suppose the functions ω(x) C2(Ω) C1(Ω) and g(x, ω) C(Ω R. Then the followings hold.

(i) If ω(x) satisfies ∆ω(x) + g(x, ω(x)) ≥ 0, x ϵ Ω, ∂ω/∂ν ≤ 0, x ϵ ∂Ω, and ω(x0) = maxΩ ω, then g(x0, ω(x0)) ≥ 0.

(ii) If ω(x) satisfies ∆ω(x) + g(x, ω(x)) ≤ 0, x ϵ Ω, ∂ω/∂ν ≥ 0, x ϵ ∂Ω, and ω(x0) = minΩ ω, then g(x0, ω(x0)) ≤ 0.

Theorem 3. Let (u(x), v(x)) be any positive solution of (2). Then

Proof. Since (u(x), v(x)) satisfies (2), then regularity theory of elliptic equations implies that u(x) and v(x) must attain their maximum and minimum values in Ω. By the first equation of (2) we have

Then Lemma 2(i) induces u(a u) ≥ 0 at the maximum of u, so maxΩ¯ u(x) ≤ a.

Then by maxΩ¯ u(x) ≤ a and the second equation of (2) we get

Using Lemma 2(i) again, there holds v(b – v βv/(a + r)) ≥ 0 at the maximum of v, and further, maxΩ¯ v(x) ≤ (a + r)/(a + r + β) < b. The proof is finished.

In the following, for convenience in use, we denote ∧ = ∧(a, b, c, m, r, α, β).

Theorem 4. Let d > 0 be a given number. If d1 ≥ d, then there is a constant C = C(d, n, Ω, ∧) such that any positive solution (u(x), v(x)) of (2) satisfies

Proof. Denote

Similar to the proof of Theorem 3, apply Lemma 2(ii) directly to the first equation of (2) to yield

which implies a ≤ u(x1) + αv(x1)/(c + u(x1) + mv(x1)). Then use Lemma 2(i) and (ii) to the second equation of (2) continuously to get

respectively. Combining with Theorem 3, we get

Thus,

where P > 0 is a constant satisfying P > bα(r + u(x2))/(cu(x2)(r +β +u(x2))). Now, let

Then h(x) and u satisfy

Further, by the Harnack inequality when d1 > d, there is C∗ = C∗(d, n, Ω, h ∞) such that

Substitute (6) into (5) to get

Combining (4) with (7), we have

Take

The result follows.

In a same way, if we take d2 as the parameter, then we can get a very similar estimate on the lower bound.

4 Coexistence

4.1 Nonexistence

This subsection devotes to the nonexistence of nonconstant positive solutions to system (2) by energy integral procedure. The following two nonexistent results are based on the priori estimates obtained in the previous section.

Theorem 5. If ( √ (am a c α)2 4am(a + c) + am a c α)/2m < C or b(a + r)/(a + r + β) < C holds, then (2) has no nonconstant positive solution, where C is given in Theorem 4.

Proof. We only prove the result holds under the first condition since the second case can be proved similarly. Suppose that (u, v) is a nonconstant positive solution to system (2). Integrating the first equation of (2) in Ω, we get

Using the positivity of u , we obtain further mC2 − (am − a − c − α)C − a(a + c) ≤ 0, which induces C ≤ (√ (am – a – c α)2 – 4am (a + c) + am – a – c – α)/2m, a contradiction occurs.

Theorem 6. Let d2 > bµ1−1 with µ1 being the second eigenvalue of −∆ with homogenous Neumann boundary condition. Then there exists a positive constant D˜ = D˜(µ1, α) or D = D(µ1, ∧) such that (2) does not permit any nonconstant positive solution if d1 > D.

Proof. Let (u, v) be a nonconstant positive solution. Denote u = (1/|Ω|) ∫Ω u(x) dx, v = (1/| Ω|) ∫ Ωv (x) dx. Multiply the first equality in (2) by u 􀀀 u and then integrate in Ω to get

Similarly, we can get

Add these two inequalities and use the Cauchy inequality to yield

where ε > 0 is arbitrary and L = max {α, βv2/((r + C)(r + u))} . Using the Poincaré equality, we have

Since d2 > bµ1−1, take ε small enough such that µ1d2 > b + εL. Then there holds

Denote D = µ−1 1(a + ε−1L). Then this inequality holds if and only if u ≡ u is a constant in view of d1 > D. This is a contradiction, and the proof is complete.

4.2 Existence

This subsection deals with the existence of nonconstant positive solutions to system (2) by choosing d2 as the parameter and by using the Leray–Schauder degree theory.

Firstly, we give some preliminaries. Denote

Denote by 0 = µ0 < µ1 < µ2 < …· the eigenvalues of −∆ under homogenous Neumann boundary condition in Ω and E(µi) the eigenspace corresponding to each µi in C2(Ω) ∩ C1(Ω). Let {ϕil, l = 1, 2, . . . , m(µi)} be an orthogonal basis in E(µi) and Xil = {Cϕil: C ∈ R2} with m(µi) being the multiplicity of µi, i = 0, 1, 2, . . . . Then Xi = ⊗l=im(µi) Xil and X = ⊗ i=0∞ Xi.

with

So, U is a solution of (2) if and only if U is a solution of (8), this is equivalent to

namely, U is a zero point of G(d1, d2; ), where (I ∆)−1 is the inverse of I ∆.

In the following, we need to use the positive constant solution of (2). By a direct analysis on the derivative of a cubic function it can be easily shown that (2) admits a unique positive constant solution denoted by (u∗, v∗) =: U ∗ if one of the following two cases holds.

Clearly, U ∗ is a zero point of G(d1, d2; . ), and the Fréchet derivative of G(d1, d2; . ) on U at U ∗ is

DU G(d1, d2; U ∗) = I − (I − ∆)−1 (D−1B + I ),

where with

It easy to check that λ is an eigenvalue of DU G(d1, d2; U ∗) on Xi if and only if

λ(1 + µi) is an eigenvalue of the matrix

Thus, DU G(d1, d2; U ∗) is invertible if and only if Mi is nonsingular. By decomposition Xi = ⊗ l = 1 m(µi) Xil it is seen that λ is an eigenvalue of DU G(d1, d2; U ∗) on each Xil, also, the multiplicity of λ is equivalent to the multiplicity of the eigenvalues λ(1 + µi) of Mi. Therefore, if λ is an q-multiplicity eigenvalue of DU G(d1, d2; U ∗) on each Xi and λ(1 + µi) is an n-multiplicity eigenvalue of Mi, then q = nm(µi).

Direct computation yields

Denote

In the following, if it involves H(d1d2; θ) but does not emphasize the effects of d1 and d2 on H, we use h(θ) to denote H(d1, d2; θ) for simplicity. Then h(µi) = d1d2 det Mi. If

then h(θ) = 0 has two real roots denoted by θ1 = θ1(d1, d2), θ2 = θ2(d1, d2), respectively, and

Denote by the set of all eigenvalues of -∆ under homogenous Neumann boundary condition in Ω. Let T = {(d1, d2) = θ: θ ≥ 0, θ2(d1, d2) < θ < θ1(d1, d2)} with θ1, θ2 being given by (11). Then we have the following lemma.

Lemma 3. If h(µi) /= 0 for µi ∈ S, then the index equality index(G(d1, d2; •), U ∗) = (−1)σ holds, where

Proof. Since h(µi) =/ 0, we know det Mi ≠ 0. So the matrix Mi is nonsingular, and DUG(d1; d2;U*) is invertible. Then index(G(d1; d2; . );U* ) = (-1) holds with y = ∑i≥0 ∑Reλi<0m(λi) and λi being an eigenvalue of DUG(d1; d2;U*) on Xi.

Now, we need to show γ = σ. Denote by τi the eigenvalue of Mi. Then τi = λi(1 + µi), m(λi) = m(τi)m(µi). The sum of algebraic multiplicity of the eigenvalues with negative real parts of Mi modulo 2 can be expressed as

where

For each µi ϵ , if µi ϵ , then det Mi < 0. By (12) we have ρi = 1. If µi ϵ T , then

det Mi > 0. By (12), again, we get ρi = 0. Therefore,

Further, γ = Σi≥0 ΣRe τi<0 m(τi)m(µi) = Σi≥0 ρim(µi) = σ.

Theorem 7. Suppose a11 > 0. If there is n > 1 such that a11d1-1 ϵ (µn, µn+1) and δn = ∑i=1 n m (µi) is odd, then there exists d* > 0 such that (2) has at least a nonconstant positive solution when d2 > d∗.

Proof. Since a11 > 0, it is easy to see that (10) holds if d2 large enough, moreover, θ1 > θ2 > 0 and limd2 →+∞ θ1 = a11d−1, limd2 →+∞ θ2 = 0. Furthermore, since a11d1−1 ∈ (µn, µn+1), there is large d0 such that

holds if d2 ≥ d0. By Theorem 5 there exists d ≥ d0 such that (2) has no nonconstant positive solution for d1 = d, d2 ≥ d. Take d large enough such that 0 < a11d−1 1 < µ1. Then there is d∗ > d such that the following holds for d2 ≥ d∗:

Now, we show that (2) has at least a nonconstant positive solution when d2 ≥ d∗. We prove this by contradiction. Suppose that there exists d∗2 ≥ d∗ such that (2) has no nonconstant positive solution. Fix d2 = d2∗ and define the homotopy operator W(t) as

Denote U = (u(x), v(x)). Consider the problem

Obviously, U is a solution of (2) if and only if U is a solution of (15) for t = 1. Especially, for any t ∈ [0, 1], U ∗ is the unique positive constant solution of (15).

For t ∈ [0, 1], U is a positive solution of (15) if and only if the equality

holds. By the discussion above we know that (16) has no nonconstant positive solution for t = 0. By our assumption (16) has no nonconstant positive solution for d2 = d∗ and t = 1, furthermore, the followings equalities hold:

where W˜ = diag(d, d∗). Then (13) and (14) induce that

Since σn is odd, Lemma 3 implies that

Thanks to Theorems 3 and 4, we know that there exist positive constants P1, P2 with P1 < C and P2 > max{a, b(a + r)/(a + r + β)} such that any positive solution (u(x), v(x)) of (15) satisfies P1 < u(x), v(x) < P2, x ∈ Ω for t ∈ [0, 1]. Define Θ = {U ∈ X: P1 < u(x), v(x) < P2, U = (u(x), v(x))T}. Then, clearly, K(U ; t) /= 0 for t ∈ [0, 1] if U ∈ ∂Θ. By the homotopy invariance of degree we get

However, the equations (U ; 0) = 0 and (U ; 1) = 0 both have a unique positive solution U ∗ in Θ, therefore,

These contradict to (17), and the proof is completed.

5 Bifurcation

In this section, taking d2 as a parameter, we investigate the existence of bifurcation solutions emitting from U ∗ by using the topological degree techniques.

For some d˜2 ϵ (0, + ∞ ), if (d˜2; U ∗) satisfies: for any δ ϵ (0, d˜2), there is d2 ϵ [d˜2 δ, d˜2 + δ] such that system (2) admits a nonconstant positive solution, then (d˜2; U ∗) ϵ (0, + ∞ ) x X is called a bifurcation point of (2). Otherwise, (d˜2; U ∗) is called a regular point of (2).

Define N (d2) = {θ > 0: h(θ) = 0} , where h(θ) is H(d1, d2; θ) defined by (9).

Obviously, N (d2) has at most two elements for any given d1, d2 > 0.

For some given d˜2 > 0, the following holds.

Theorem 8.

(i) If S ∩ N (d˜2) = ∅, then (d˜2; U ∗) is a regular point of (2).

(ii) Suppose S ∩ N (d˜2) ≠ ∅ and (bd1 + a11d˜2)2 ≠ 4d1d˜2(a11b − a12a21). If Σµi ∈N (d˜2 ) m(µi) is odd, then (d˜2; U ∗) is a bifurcation point of (2).

Proof. Let W (x) = U (x) − U ∗. Then (8) is equivalent to

and then it is equivalent to W satisfying

The Fréchet derivative of R(d2; •) on W at W = 0 is

DW R(d2; 0) = I − (I − ∆)−1 D−1B + I = DU G(d1, d2; U ∗),

where B is given in Section 4.2. As what we do in Section 4, λ is an eigenvalue of DW R(d2; 0) on Xi if and only if λ(1 + µi) is an eigenvalue of Mi.

(i) If S ∩ N (d˜2) = ∅, that is, µi ∈ S, H(d1, d˜2; µi) ≠ 0, then det Mi ≠ 0, and Mi is nonsingular, so DW (d2; 0) is regular. By the implicit function theorem we know that there exists ε, 0 < ε 1, such that W = 0 is the unique solution of (d2; ) = 0 in Ξ = {U ϵ X: ││U││ X < ε} for any d2 (d2 is in a neighborhood d˜2), this means that (d˜2; U ∗) is a regular point of (2).

(ii) By assumptions we take µi ϵ (d˜2), then H (d1, d˜2; µi) = 0 and µi ≠ (a11d1−1 + bd2−1)/2. Further,

Thus, the two eigenvalues of Mi are 0 and tr Mi, respectively, and 0 is simple.

If the result is false, then there is d˜2 > 0 satisfying the following.

(a) S ∩ N (d˜2) ≠ ∅, (bd1 + a11d˜2)2 ≠ 4d1d˜2(a11b − a12a21) and Σµ ∈N (d˜ ) m(µi) is odd

(b) There exists δ ∈ (0, d˜2) such that W = 0 is the unique solution of R(d2; .•) = 0 in Ξ for d2 ∈ [d˜2 − δ, d˜2 + δ].

Now, it needs to deduce a contradiction. Since R(d2; •) is a compact perturbation of I, then (b) implies that the Leray–Schauder degree deg(R(d2; •), Ξ, 0) is well defined, furthermore, deg(R(d2; •), Ξ, 0) is independent of d2 ∈ [d˜2 − δ, d˜2 + δ]. If DW R(d2; 0) is invertible for d2 ∈ [d˜2 − δ, d˜2 + δ], then

where γ(d2) is the sum of algebraic multiplicity of the eigenvalues with negative real parts of DW R(d2; 0). For any µi ∈ S ∩ N (d˜2), we have

The Fréchet derivative of H(d1, d2; µi) on d2 at d˜2 is Hd (d1, d˜2; µi) = d1µ2 − a11µi.

It is easy to see that H d2 (d1, d˜2; µi) ≠ 0. (Otherwise, if H d2 (d1, d˜2; µi) = 0, then µi = a11d1−1. Substituting µi = a11d1−1 into (19), we get a12a21 = 0, which contradicts to a12a21 < 0.) Therefore, for all d2 ∈ [d˜2 − δ, d˜2 + δ], we have

Combine (20) with (19) to get

Clearly, S has no accumulation point. So, S ∩ N (d2) = Φ for d2 ϵ [d˜2 δ, d˜2) (d˜2, d˜2 + δ] and small δ. Then by the proof of (i) we know that DW R (d2; 0) is invertible when d2 ϵ [d˜2 δ, d˜2) U (d˜2, d˜2 + δ].

The proof of Lemma 3 tells us that the sum of algebraic multiplicity of eigenvalues with negative real parts of DW R(d2; 0) is Re λi<0 m(λi) =: γi(d2), then

where ρi is defined by (12).

If µi ∈/ N (d˜2), then γi(d2) is independent of d2 ∈ [d˜2 − δ, d˜2 + δ] since deg( R (d2; . ), Ξ, 0) is independent of d2. If µi ϵ (d˜2), then (21) and (22) imply that the sum of algebraic multiplicity of the eigenvalues with negative real parts of DW R(d˜2 − δ; 0) and DW R(d˜2 + δ; 0) on Xi satisfies

Hence,

Since ∑µiϵN(d2) m(µi) is odd, by (18) we get

However, deg(R(d2; . ); Ξ; 0) is independent of d2 ϵ [d2 δ; d2+ δ], so (23) cannot hold.

Therefore, ( d2;U*) is a bifurcation point of (2). The proof is accomplished.

6 Numerical examples, discussions and conclusions

In the previous sections, we investigate the coexistence of prey and predator and find out their coexistent conditions. Nevertheless, from the realistic point of view, we do hope the prey and predator can coexist since our aim is to maintain the balance of ecosystems. The conditions in Theorems 7 and 8 ensure that the prey and predator can coexist and, in this case, the ecological balance may be sustainable. To clarify this point, in this section, we give some numerical examples to support our theoretical analysis. We fix a, b, c, m, r and take , as parameters. For the sake of simple calculation, the parameters a, b, c, m and r are referred to [18] in part and taken as follows: a = 1:2, b = 0:08, c = 2:4, m = 50, r = 12. The parameters , are chosen, and the values (u*; v*), ac, b, br= (cr + c + bmr) =: s1, c + r + + bm =: s2, a11 are calculated and listed in Table 1.

Summarizing these numerical simulations, the following results follow: In Examples 1–2, we see that ac > b holds, then Theorem 2 shows that system (1) has persistent property and any positive solution (u(x; t); v(x; t)) of (1) satisfies

respectively. However, ac < b holds in Examples 3–5, so, we do not know whether system (1) has persistence or not in these three cases. Nevertheless, we also see that the existent condition (H1) of (u*; v*) in Examples 1–5 holds and (u*; v*) is worked out in each example, and then (2) may have a nonconstant positive solution according to Theorem 7.

Based on the parameter values taken above, we depict some trajectory graphs or spatiotemporal pattern formation to simulate our theoretical results.

The followings Figs. 1–5 and Figs. 6–10 are, respectively, the trajectory graphs and the spatiotemporal patterns corresponding to Examples 1–5 with d1 = 0:1, d2 = 1 and different initial values.

All the trajectory graphs show that the quantity of prey and predator decrease and increase gradually as the predation rate of predator, the conversion rate of predator due to capturing prey increasing and time going on, and tend to stabilize eventually. Specifically, the quantity of prey and predator are as followings in turn: 0:08, 0:07, 0:06, 0:05, 0:04 and 0:22, 0:41, 0:65, 0:84, 0:98, respectively, from Figs. 1,2,3,4,5. We also see that the declining and ascending scales between the prey and predator are inconformity, in contrast, the increase of the predator is much stronger than the decrease of the prey. Though the quantity of the prey is decreasing, it is clear that such changes are very small, that is, the quantity of the prey remains roughly stable regardless some fluctuations of the predation rate and the conversion rate of predator. This situation probably implies that the available food resource for predator is not single, that is, what we want to see in population dynamics since which sustains the persistence of the species when the food resource for predator is scarce. This helps to maintain the ecological balance. Figures 6,7,8,9,10 show that the spatiotemporal patterns always occur for our given parameter values, which imply that the system exhibits rich dynamical behavior.

In fact, we made a lot of numerical examples and trajectory graphs, we found that the differences among these simulation and graphic results are very small. Here, we only

present the above five numerical examples and their corresponding trajectory graphs, the spatiotemporal pattern formation to state our findings.

The numerical examples tell us that the ecosystem reflected by model (1) or (2) is easy to maintain stability when the predation rate and the conversion rate of predator do not change dramatically. This might imply that neither the prey nor the predator will disappear within a certain period of time, which could also mean that other species is difficult to invade such system, this is exactly what we hope to happen biologically.

References

1 N. Apreutesei, G. Dimitriu, On a prey–predator reaction–diffusion system with Holling type III functional response, J. Comput. Appl. Math., 235(2):366–379, 2010, https://doi.org/10.1016/j.cam.2010.05.040

2 M.A. Aziz-Alaoui, Study of a Leslie–Gower-type tritrophic population model, Chaos Solitons Fractals, 14(8):1275–1293, 2002, https://doi.org/10.1016/S0960-0779(02)00079-6

3 J.R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, J. Anim. Ecol., 44(1):331–340, 1975, https://doi.org/10.2307/3866

4 R.S. Cantrell, C. Cosner, Models for predator–prey systems at multiple scales, SIAM Rev.,38(2):256–286, 1996, https://doi.org/10.1137/1038041

5 D.L. DeAngelis, R.A. Goldstein, R.V. O’Neill, A model for trophic interaction, Ecology,56(4):881–892, 1975, https://doi.org/10.2307/1936298

6 Y. Du, Y. Lou, Qualitative behaviour of positive solutions of a predator–prey model: Effects of saturation, Proc. R. Soc. Edinb., Sect. A, Math., 131(2):321–349, 2001, https://doi.org/10.1017/S0308210500000895

7 R.P. Gupta, M. Banerjee, P. Chandra, Bifurcation analysis and control of Leslie–Gower predator–prey model with Michaelis–Menten type prey-harvesting, Differ. Equ. Dyn. Syst., 20(3):339–366, 2012, https://doi.org/10.1007/s12591-012-0142-6

8 N.G. Hairston, F.E. Smith, L.B. Slobodkin, Community structure, population control and competition, Am. Nat., 94(879):421–425, 1960, https://doi.org/10.1086/282415

9 S.-B. Hsu, T.-W. Huang, Global stability for a class of predator–prey systems, SIAM J. Appl. Math., 55(3):763–783, 1995, https://doi.org/10.1137/S0036139993253201

10 R.E. Kooij, A. Zegeling, A predator-prey model with Ivlev’s functional response, J. Math. Anal. Appl., 198(2):473–489, 1996, https://doi.org/10.1006/jmaa.1996.0093

11 P.H. Leslie, Some further notes on the use of matrices in population mathematics, Biometrika, 35(3-4):213–245, 1948, https://doi.org/10.1093/biomet/35.3-4.213

12 P.H. Leslie, J.C. Gower, The properties of a stochastic model for the predator-prey type of interaction between two species, Biometrika, 47(3–4):219–234, 1960, https://doi.org/10.1093/biomet/47.3-4.219

13 A.J. Lotka, Elements of physical biology, J. Am. Stat. Assoc., 20(151):452–456, 1925, https://doi.org/10.2307/2965538

14 Y. Lou, W.-M. Ni, Diffusion, self-diffusion and cross-diffusion, J. Differ. Equations, 131(1): 79–131, 1996, https://doi.org/10.1006/jdeq.1996.0157

15 A. Moussaoui, M.A. Azi-Alaoui, R. Yafia, Permanence and periodic solution for a modified Leslie–Gower type predator–prey model with diffusion and non constant coefficients, Biomath., 6(1):1707107, 2017, https://doi.org/10.11145/j.biomath.2017. 07.107

16 A.F. Nindjin, M.A. Aziz-Alaoui, M. Cadivel, Analysis of a predator-prey model with modified Leslie–Gower and Holling-type II schemes with delay, Nonlinear Anal., Real World Appl., 7(5):1104–1118, 2006, https://doi.org/10.1016/j.nonrwa.2005.10.003

17 A.F. Nindjin, K.T. Tia, H. Okou, A. Tetchi, Stability of a diffusive predator-prey model with modified Leslie-Gower and Holling-type II schemes and time-delay in two dimensions, Adv. Difference Equ., 2018:177, 2018, https://doi.org/10.1186/s13662-018-1621-z

18 F.A. Rihan, S. Lakshmanan, A.H. Hashish, R. Rakkiyappan, E. Ahmed, Fractional-order delayed predator–prey systems with Holling type-II functional response, Nonlinear Dyn., 80(1–2):777–789, 2015, https://doi.org/10.1007/s11071-015-1905-8

19 M.L. Rosenzweig, Paradox of enrichment: Destabilization of exploitation ecosystems in ecological time, Science, 171(3969):385–387, 1971, https://doi.org/10.1126/science.171.3969.38520

20 E. Sáez, E. González-Olivares, Dynamics of a predator–prey model, SIAM J. Appl. Math.,59(5):1867–1878, 1999, https://doi.org/10.1137/S0036139997318457

21 A. Sharma, A.K. Sharma, K. Agnihotri, Analysis of a toxin producing phytoplankton–zooplankton interaction with Holling IV type scheme and time delay, Nonlinear Dyn., 81(1–2): 13–25, 2015, https://doi.org/10.1007/s11071-015-1969-5

22 A.R.E. Sinclair, S. Mduma, J.S. Brashares, Patterns of predation in a diverse predator– prey system, Nature, 425(6955):288–290, 2003, https://doi.org/10.1038/nature01934

23 J. Wang, J. Shi, J. Wei, Predator–prey system with strong Allee effect in prey, J. Math. Biol.,62(3):291–331, 2011, https://doi.org/10.1007/s00285-010-0332-1

24 W. Wang, L. Zhang, H. Wang, Z. Li, Pattern formation of a predator–prey system with Ivlevtype functional response, Ecol. Model., 221(2):131–140, 2010, https://doi.org/10.1016/j.ecolmodel.2009.09.011

25 H. Yin, X. Xiao, X. Wen, Analysis of a Lévy-diffusion Leslie–Gower predator–prey model with nonmonotonic functional response, Discrete Contin. Dyn. Syst., Ser. B, 23(6):2121–2151, 2018, https://doi.org/10.3934/dcdsb.2018228