Stability and Coexistence of a Diffusive Predator-Prey System with Nonmonotonic Functional Response and Fear Effect

This paper investigates the diﬀusive predator-prey system with nonmonotonic functional response and fear eﬀect. Firstly, we discussed the stability of the equilibrium solution for a corresponding ODE system. Secondly, we established a priori positive upper and lower bounds for the positive solutions of the PDE system. Thirdly, suﬃcient conditions for the local asymptotical stability of two positive equilibrium solutions of the system are given by using the method of eigenvalue spectrum analysis of linearization operator. Finally, the existence and nonexistence of nonconstant positive steady states of this reaction-diﬀusion system are established by the Leray–Schauder degree theory and Poincar´e inequality.


Introduction
In order to describe the evolution of biological populations in the ecosystem, some mathematical theories and methods have been used to establish the corresponding biological mathematical model, which has become a research hotspot. In recent years, the research on biological models such as the predator-prey model has aroused the attention of many scientists and biologists. e predator-prey model of PDE forms is an important branch of reaction-diffusion equations. e dynamic relationship between predator and their prey is one of the dominant themes in ecology and mathematical ecology. During these thirty years, the investigation on the prey-predator models has been developed, and more realistic models are derived in view of laboratory experiments. Moreover, the research on the prey-predator models has been studied from various views and obtained many good results (see  and the references therein).
However, many studies have shown that only the presence of predators in front of the prey can affect the size of the prey population, and the effect is even greater than the effect of direct predation. Although some biologists have realized that the relationship between the prey and the predator cannot be simply described as direct killing, we should take the fear of the prey population into account. At present, there are few research studies on establishing corresponding mathematical models to explain this phenomenon.
For every specific prey-predator system, we know that the functional response of the predator to the prey density is very important, which represents the specific transformation rule of the two organisms. In [8], Pang and Wang considered a predator-prey model incorporating a nonmonotonic functional response which is called the Monod-Haldane or Holling type IV function: where z n is the outward directional derivative normal to zΩ. Model (1) describes a prey population u which serves as food for a predator with population v. e parameters r, d, m, p, c, and k are assumed to be only positive values: the positive constant k is the carrying capacity of the prey and the positive constant m is the death rate of the predator; r is the growth rate of prey u; and the positive constants d 1 and d 2 are the diffusion coefficients.
In this paper, based on the above model, in order to describe the evolution law of the population in the ecosystem more specifically, we will consider the natural mortality and fear effect of the prey population and establish the corresponding PDE model within a fixed bounded domain Ω ⊂ R N with smooth boundary at any given time and the natural tendency of each species to diffuse to areas of smaller population concentration [7][8][9][10]. Hence, we will investigate the following reaction-diffusion system under the homogeneous Neumann boundary conditions as follows: where u 0 and v 0 are continuous functions of x. u and v stand for the densities of prey and predators, respectively. e parameters a, b, c, d, r, k, p, and m are assumed to be only positive constants. a and m denote the intrinsic death rate of prey u and predator v, respectively. k stands for the fear factor of prey to predator. e remaining parameters refer to (1). Here, f(u, v) � (uv/d + u 2 ) stands for Monod-Haldane functional response. e main aim of this paper is to study the nonconstant positive steady states of (2), that is, the existence and nonexistence of nonconstant positive classical solutions of the following elliptic system: e rest of this paper is arranged as follows. In Section 2, we discuss the stability of the equilibrium of the ODE system which corresponds to system (2). In Section 3, we establish a priori positive upper and lower bounds for the positive solutions of the PDE system. In Section 4, sufficient conditions for the local asymptotical stability of two positive equilibrium solutions of the system are established by using the method of eigenvalue spectrum analysis of linearization operator. In Section 5, the existence and nonexistence of nonconstant positive steady states of this reaction-diffusion system are established by using the Leray-Schauder degree theory, which demonstrates the effect of large diffusivity.

Stability of the ODE Model
e goal of this section is to discuss the stability of the ODE model; we give the ordinary differential equation of system (3) as follows: By the similar method to [7], for (4), we can get the following result.
Next, we will calculate the equilibrium point of system (4), and the result is given as follows.
Theorem 1. System (4) always has an extinction equilibrium point E 0 � (0, 0). If r > a, then system (4) has only the equilibrium point Proof. It is easy to see that all equilibrium points of system (4) satisfy the following equations: It follows that system (4) obviously has equilibrium points E 0 � (0, 0) and E 1 � (r − a/b, 0) with r > a. Next, we consider the existence of positive constant equilibrium point E 2,i . By calculating the second equation of (6), we directly get 2 Complexity where (cp − 2m (6) and combining the two equations of system (6), we can obtain the following equation: rough the same solution deformation calculation, we can get where Obviously, if a 0 ≤ 0, then f(v) � 0 has no positive constant solution; if a 0 > 0, then f(v) � 0 has only one positive constant solution. anks to the same sign a 0 and r − a − bu i , a 0 > 0 implies r − a − bu i > 0, which ensures that f(v) � 0 has only one positive constant solution denoted by v i . us, system (4) has two positive constant equilibrium Proof. e proof of eorem 2 is similar to that of eorem 2 of [9]; hence, we omit it.
rough mathematical calculation, we obtain the Jacobian matrix of system (5) at the equilibrium point E 1 � (r − a/b, 0) as follows: Obviously, when (r − a)(cpb + ma − mr) − mb 2 d < 0 and both eigenvalues of J E 1 have negative real parts, then E 1 � (r − a/b, 0) is locally asymptotically stable; when , the corresponding Jacobian matrix is given by By simplifying, we can get where It is easy to get that det J E 2,i > 0 and trac en, two eigenvalues of the matrix J E 2,i have negative real parts. erefore, the equilibrium E 2,i � (u i , v i ) is locally asymptotically stable. If u i < (2m 2 v i /c 2 pb) and the matrix J E 2,i has one positive eigenvalue, then E 2,i is unstable.

A Priori Estimates on Equation (3)
e main purpose of this section is to give a priori upper and lower bounds for the positive solutions. To this aim, we first recall the following maximum principle due to [23,24].

Theorem 5.
If r > a and ack is a positive solution of (3). en, the solution (u(x), v(x)) of (3) yields Proof. By Lemma 2, if u reaches its maximum at x ∈ Ω, it follows from the first equation of (3) that that is, If ϖ reaches its maximum at x 0 ∈ Ω, then which results in anks to ϖ � cd 1 u + d 2 v, we know that that is, If φ reaches its maximum at x 1 ∈ Ω, then where N � (2cp(r − a)/b d)ℵ, which means that Letting h(u) � cru/acu + bcu 2 + N, it is easy to get the maximum of h(u), that is, h(u) < r/a. us, v(x 1 ) > (r − a/ak): By (23) and (28), we have proved eorem 5. According to eorems 5 and 1, we can easily get the following conclusion.

Theorem 7. Suppose that (u, v) is a nonnegative classical solution of (3). If r ≤ a, then (u, v) is always zero solution.
Proof. Integrating the equation for u in (3) over Ω by parts, we get us, 4 Complexity Hence, u ≡ 0. Substituting u � 0 into the second equation of (3), we get and then, v ≡ 0. e proof is complete.

Stability of the Equilibrium of Equation (3)
e goal of this section is to investigate the local and global stability of the positive constant steady state (u i , v i ) � U. We first discuss the local stability of U i . To this end, we need to introduce some notations for developing our result. Let (32) erefore, system (3) becomes the following forms: It follows that two positive solutions (u i , v i ) (i � 1 and 2) satisfy where v i (i � 1, 2) satisfies In order to get the linearization operator of (3) at the positive constant solution (u i , v i ), for (33), we calculate the partial derivatives with respect to u and v, respectively, at the equilibrium point (u i , v i ), as follows: (37) Next, give some results as follows: (i) 0 � μ 0 < μ 1 < μ 2 < μ 3 < · · · < μ i · · · < ∞ are the eigenvalues of − Δ on Ω under homogeneous Neumann boundary condition, and m i is the algebraic multiplicity of eigenvalue μ i . (ii) ϕ ij , 1 ≤ j ≤ m i , are the normalized eigenfunctions corresponding to μ i , and then ϕ ij (i ≥ 0, 1 ≤ j ≤ m i ) are the orthonormal basis of L 2 (Ω). (39) (2) If f 1 < 0, then the positive constant steady state (u 2 , v 2 ) of (3) is locally asymptotically stable; if f 1 > 0, then the positive constant steady state (u 2 , v 2 ) of (3) is unstable.
e linearization operator of (3) at the positive constant solution (u i , v i ) (i � 1 and 2) can be written as where f 1 , f 2 , g 1 , and g 2 are defined in (36)-(37). According to the linear stability theory, if the real parts of all eigenvalues of L σ are negative, then (u i , v i ) is locally asymptotically stable; if there exists the positive real part of the eigenvalue of L σ , then (u i , v i ) is unstable. Let (ϕ(x) and ψ(x)) be the eigenfunctions corresponding to the eigenvalue λ. en, that is, us, the eigenvalue equation of system (3) is equivalent to λ is an eigenvalue of L σ if and only if there exists i ≥ 0 such as det (B i ) � 0, which is equivalent to where Next, we check the stability of (u 1 , v 1 ) and (u 2 , v 2 ), respectively.

Nonconstant Positive Steady States of Equation (3)
e main purpose of this section is to provide some sufficient conditions for the existence and nonexistence of a nonconstant positive solution of (3) by using the Leray-Schauder degree theory [12,24,25]. Next, we will establish these results by dividing into two sections.

Nonexistence.
e goal of this part is to establish some sufficient conditions for the nonexistence of nonconstant positive solutions of (3) by the energy norm method. Some related research studies can refer to [8][9][10]. For the ease of notation, we set where (u, v) is a positive solution of (3).
Multiplying the second equation of v by χ and integrating over Ω by parts, we obtain anks to the boundary of u and v (see in eorem 5), we get 6 Complexity where C 1 � cp dℵ > 0 and C 2 � m − cp/2 . Applying Cauchy inequality, we obtain Substituting (50) into (49) and using Poincaré inequality, we get Because u and v are nonnegative, we obtain Multiplying the above equation of u by ω and integrating over Ω by parts, using Poincaré inequality again, we obtain that is, which implies that ω and χ are always constant. e proof is complete. For simplicity, we only consider the existence of nonconstant positive classical solutions near (u 2 , v 2 ) which are denoted by (u, tv). (3) can be written as follows: where f 3 ≐ o(|u|, |v|) and g 3 ≐ o(|u|, |v|). Define the space S and E as follows: Set U � (u, tv), and then, (56) becomes where (59) If the principal eigenvalue μ 1 has an odd multiple eigenfunction and d 2 > d (1) 2 , then system (3) has at least one nonconstant positive solution.
Proof. It is easy to see that system (3) has no solution on the boundary of the space S. According to Homotopy invariance of degree theory, for all d 1 > 0, deg(I − K − H, E ∩ S, 0) is well defined and constant. Next, we will prove Assume that (0, 0) is an isolated fixed point of I − K − H, then where τ is the sum of algebraic multiplicity of all eigenvalues greater than 0. Assume that λ is the eigenvalue of K − I and the corresponding eigenfunction is denoted by (ξ, η), then − d 1 (λ + 1)Δξ � (1 − λ)f 1 ξ + f 2 η, − d 2 (λ + 1)Δη � g 1 ξ + λg 1 η. (62) us, the eigenvalue equation of system (3) is equivalent to where us, all eigenvalues of K − I satisfy Complexity 7 Set P i � 2d 1 d 2 μ 2 i + d 1 μ i g 1 − f 1 g 1 and Notice that f 2 < 0, g 1 > 0, and f 1 > 0 are defined in (36)-(37), and it follows that L σ and K(d 1 ) − I have the same number of eigenvalues. anks to eorem 9, let d 1 � max r − a/μ 1 , f 1 /μ 1 + 1, then Next, we will calculate the sum of algebraic multiplicity of all eigenvalues of K − I greater than 0.