Dynamics and pattern formation in homogeneous diffusive predator–prey systems with predator interference or foraging facilitation

Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa

Dynamics and pattern formation in homogeneous diffusive predator–prey systems with predator interference or foraging facilitation Xuebing Zhang a ,∗, Honglan Zhu b a

College of Mathematics and Statistics, Nanjing University of Information Science and Technology, Nanjing 210044, People’s Republic of China b Business School, Huaiyin Institute of Technology, Huaian 223003, People’s Republic of China

A R T I C L E

I N F O

Article history: Received 27 September 2018 Received in revised form 18 January 2019 Accepted 21 January 2019 Available online xxxx Keywords: Diffusion Persistence Stability Spatial pattern

A B S T R A C T

We study a general diffusive predator–prey system under Neumann boundary conditions. First, we examine the global attractor and persistence of the system, which characterize the long-time behavior of the time-dependent solution, and the stability of all nonnegative equilibria of the system. Then, by analyzing the associated characteristic equation, we derive explicit conditions for the existence of Hopf bifurcation and Turing instability. Furthermore, the existence and nonexistence of nonconstant positive steady states of this model are studied by considering the effect of large diffusivity. Finally, to verify our theoretical results, the selected numerical simulations are included. The results show that under the effect of diffusion, the system shows different spatial patterns that can be stationary, periodic and chaotic. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Let u and v denote the size of the prey and predator population, respectively. When the prey and predator densities are spatially inhomogeneous in a bounded domain Ω ⊂ RN with a smooth boundary, we consider the following modified Rosenzweig–MacArthur predator–prey model with the predator-dependent functional response: ⎧ u λ(v)u ∂u ⎪ ⎪ = d1 ∆u + ru(1 − ) − v (x, t) ∈ Ω × (0, +∞), ⎪ ⎪ ∂t K 1 + hλ(v)u ⎪ ⎪ ⎪ ⎪ λ(v)u ⎨ ∂v = d2 ∆v + e v − mv, (x, t) ∈ Ω × (0, +∞), (1.1) ∂t 1 + hλ(v)u ⎪ ⎪ ⎪ u(x, 0) = φ(x), v(x, 0) = ψ(x), x ∈ Ω , ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂u(x, t) = ∂v(x, t) = 0, x ∈ ∂Ω . ∂ν ∂ν In the above equation, ν is the normal vector that extends past the bounded domain Ω . The homogeneous Neumann boundary condition means that model (1.1) is self-contained and has no population flux across ∗ Corresponding author. E-mail addresses: [email protected] (X. Zhang), [email protected] (H. Zhu).

https://doi.org/10.1016/j.nonrwa.2019.01.016 1468-1218/© 2019 Elsevier Ltd. All rights reserved.

268

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

the boundary ∂Ω . All the parameters appearing in model (1.1) are assumed to be positive constants. The constants d1 and d2 are the diffusion coefficients corresponding to u and v, and the initial data φ(x) and ψ(x) are continuous functions. Since the variables u and v account for the densities of prey and predator, respectively, they are required to be nonnegative. The parameter h is the predator handling time of one prey. The function λ(v) is an encounter rate function that represents the predator interference (λ′ (v) < 0) or foraging facilitation (λ′ (v) > 0). The corresponding ODE version was studied in [1–3]. Berec [1] considered the specific encounter rate function λ0 , (1.2) λ(v) = (b + v)ω where b ≥ 0, λ0 > 0 and ω > 0 correspond to interference, ω = 0 recovers the type II functional response, λ(v)u and ω < 0 represents facilitation. Clearly, the term 1+hλ(v)u covers the Hassell–Varley functional response for b = 0, including the form used by Cosner [4] if ω = −1 and the ratio-dependent functional responses if ω = 1. In addition, this term covers the Beddington–DeAngelis functional response if b > 0 and ω = 1. Hence, all the common functional responses that model predator interference can be viewed as arising from encounter rates that decrease with increasing predator density. In addition, it is easily obtained that for a negative ω, λ(v) is an increasing unbounded function from above, concave for ω ∈ (−1, 0) and convex for ω < −1. Note that for ω > 0, the encounter rate (1.2) models a negative interaction. However, Berec [1] considered only the number and character (i.e., stability) of coexistence equilibria. Pˇribylov´ a [2] extended this study by looking more thoroughly at the system dynamics when some of the coexistence equilibria are unstable and characterized system bifurcations that occur to obtain an even more general predator-dependent functional response. Furthermore, Pˇribylov´ a [3] proposed an encounter rate function as follows: λ(v) = L −

λ0 , (b + v)ω

(1.3)

where b ≥ 0 λ0 > 0, ω > 0 and λ(0) ≥ 0. Encounter (1.3) is increasing, concave and bounded from above. λ(v)u For a large population density of predators, (1.3) saturates at L, and the functional response 1+hλ(v)u asymptotically approaches that of Holling type II. In both the specifically chosen encounter rate forms (1.2) and (1.3), one may interpret the parameter ω as a measure of the facilitation rate among the predator population. For ω = 0, there is no facilitation, and the encounter rate is not dependent on the predator population density. With increasing ω, the encounter rate for the same predator population density increases. In [3], using the tools of bifurcation analysis, the authors described all the nonlinear phenomena that occur in the system provoked by foraging facilitation, including the fold, Hopf, transcritical, homoclinic and Bogdanov–Takens bifurcation. Foraging facilitation was shown to stabilize the coexistence in the predator– prey system for specific rates, but in most cases, such facilitation can have fatal consequences for the predators themselves. On the other hand, the spatial component of ecological interactions has been identified as an important factor in determining the distribution of populations with the advancement of time and the structure of communities over their habitats. Empirical evidence suggests that movement rates of interacting individuals can dramatically affect population stability [5,6] and the composition of communities [7,8]. Furthermore, such rates can cause the system to generate a wide variety of spatial patterns such the Turing pattern and chaotic pattern, among others. Recently, nonconstant steady states have been suggested as a fundamental argument for the existence of spatial patterns and have received great attention from various researchers [9–18]. Therefore, in this paper, the main aim is to study the large-time behavior of solutions to (1.1) and the Hopf bifurcation, Turing instability, and existence and nonexistence of nonconstant positive steady states of (1.1). ¯ + ), λ(0) > 0, and In this paper, we assume that λ(v) satisfy the following conditions: (A) λ(v) ∈ C 1 (R λ(v) is monotonic. Clearly, the functions (1.2) and (1.3) satisfy condition (A).

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

269

For completeness, we recall well-known results such as the existence and local asymptotic stability conditions for coexisting steady states that are relevant to the forthcoming discussion of our present work. (1) The ODE system has two boundary equilibria: E0 = (0, 0) and E1 = (K, 0). The extinction equilibrium E0 is always a saddle point. The pre-only equilibrium E1 is locally stable for any K if e < mh and for K < m/[λ(0)(e − mh)] if e > mh; if K < m/[λ(0)(e − mh)] for e > mh, then E1 is a saddle point. (2) The coexistence equilibrium E ∗ = (u∗ , v ∗ ) satisfies ⎧ C2 ⎪ , ⎪ u∗ = K ⎨ λ(v ∗ ) ∗ ⎪ ⎪ v ∗ = C λ(v ) − C2 , ⎩ 1 λ2 (v ∗ ) er where C1 = e−hm and C2 = Therefore, we have

m K(e−hm) .

(1.4)

Hence, if e − hm ≤ 0, then no coexistence equilibrium exists. e − hm > 0, λ(v ∗ ) > C2 .

(1.5)

In [1], the authors obtain that if e < mh, then the prey-only equilibrium [K, 0] is locally asymptotically stable. However, for the PDE system (1.1), we can claim that [K, 0] is globally asymptotically stable under this condition; see Theorem 3.1. In addition, we obtain that the unique positive constant steady state E ∗ = (u∗ , v ∗ ) is globally asymptotically stable under certain conditions; see Theorem 3.3. Theorems 3.4 and 3.5 showed that the diffusion we incorporated into system (1.1) has a great effect on the dynamics on system (1.1). Such diffusion can cause not only Turing stability but also Hopf bifurcation of system (1.1). The organization of the remaining part of the paper is as follows. In Section 2, we investigate the long-time behavior of the time-dependent solution of (1.1), in particular, its dissipative and persistent properties. In Section 3, the stability of all nonnegative constant solutions of (1.1) is studied. In Section 4, by the energy method, we present some nonexistence results of nonconstant positive classical solutions of (4.1) for a certain range of the parameters; based on degree theory, we consider the existence of nonconstant positive classical solutions for (4.1). Our results demonstrate that spatial patterns can be found as a result of diffusion. 2. Global attractor and permanence First, we present the following lemma, which is used to prove the dissipation of system (1.1) based on Theorem 3.1 in [19] by using the L1 bound. Lemma 2.1. Consider the following system: ⎧ ∂w ⎪ = d∆w + wB(x, t) ⎪ ⎪ ⎪ ⎨ ∂t ∂w = 0, ⎪ ⎪ ⎪ ∂ν ⎪ ⎩ w(x, 0) = w0 (x),

(x, t) ∈ Ω × (0, +∞), x ∈ ∂Ω ,

(2.1)

¯. x∈Ω

Suppose that B(x, t) is uniformly bound and local Lipschitz compliant in (x, t). Assume that w(x, t) ≥ 0 and that supt≥0 ∥w(·, t)∥ < M , where M is a finite positive constant. Then, sup ∥w(·, t)∥L∞ (Ω) < M ∗ , t≥0

where M ∗ depends on M and on ∥w0 (·)∥L∞ (Ω) .

270

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Theorem 2.1. For system (1.1), we have (a) If u0 (x, t) ≥ 0 (̸≡ 0) and v0 (x, t) ≥ 0 (̸≡ 0), then system (1.1) has a unique solution (u(x, t), v(x, t)) such ¯. that u(x, t) > 0 and v(x, t) > 0 for t ∈ (0, ∞) and x ∈ Ω (b) Any solution (u, v) of (1.1) satisfies ∫ eK(m + r) |Ω |. lim sup u(t, x) ≤ K, v(x, t)dx ≤ m t→+∞ Ω Moreover, ∗ ∥u(·, t)∥C(Ω) ¯ ≤ K1 , ∥v(·, t)∥C(Ω) ¯ ≤C ,

where K1 = max{K, maxΩ¯ φ(x)} and C ∗ depends on r, K, b, λ0 , ω, φ(x), ψ(x) and Ω . ¯. on Ω (c) If d1 = d2 , then lim supt→+∞ maxΩ¯ v(x, t) ≤ lim supt→+∞ maxΩ¯ w(x, t) ≤ Ke(m+r) m Proof . Proof of (a): Define f (u, v) =

λ(v)u u , F (u, v) = ru(1 − ) − f (u, v)v, 1 + hλ(v)u K

G(u, v) = ef (u, v)v − mv. λ(v) ∂f ∂f ∂G Then, we can obtain that ∂F ∂v = − ∂v v − f and ∂u = ev ∂u = ev (1+hλ(v)u)2 > 0. It is easily obtained that the system ∂F is a quasimonotone nondecreasing system if ∂F ∂v < 0 and that the system is a mixed quasimonotone system if ∂v ≥ 0. Consider the following system: ⎧ du u λ(v)u ⎪ ⎪ = ru(1 − ) − v, ⎪ ⎪ dt K 1 + hλ(v)u ⎪ ⎨ λ(v)u dv (2.2) =e v − mv, ⎪ ⎪ ⎪ dt 1 + hλ(v)u ⎪ ⎪ ⎩ u(0) = u0 , v(0) = v0 .

Let (u(t; u0 , v0 ), v(t; u0 , v0 )) be the unique solution of system (2.2) and min φ(x) = ϕm , max φ(x) = ϕM , min ψ(x) = ϕm , max ψ(x) = ϕM . ¯ Ω

¯ Ω

¯ Ω

¯ Ω

(2.3)

Let φm > 0; if system (1.1) is a quasimonotone nondecreasing system, then (u(t; φm , ψm ), and v(t; φm , ψm )) and (u(t; φM , ψM ), v(t; φM , ψM )) are the lower solution and upper solution of system (1.1), respectively. Therefore, system (1.1) has a unique globally defined solution (u(x, t), v(x, t)) satisfying u(t; φm , ψm ) ≤ u(x, t) ≤ u(t; φM , ψM ), v(t; φm , ψm ) ≤ v(x, t) ≤ v(t; φM , ψM ). If system (1.1) is a mixed quasimonotone system, then (u(t; φm , ψM ), v(t; φm , ψm )) and (u(t; φM , ψm ), v(t; φM , ψM )) are the lower solution and upper solution of system (1.1), respectively. Therefore, system (1.1) has a unique globally defined solution (u(x, t), v(x, t)) satisfying u(t; φm , ψM ) ≤ u(x, t) ≤ u(t; φM , ψm ), v(t; φm , ψm ) ≤ v(x, t) ≤ v(t; φM , ψM ). If φm = 0(φ(x) ̸≡ 0), then u(x, t) > 0 for t > 0. Choose δ > 0; then, minΩ¯ u(x, δ) > 0. Letting minΩ¯ u(x, δ) as the initial value of system (2.2), we obtain that system (1.1) has a unique globally defined solution (u(x, t), v(x, t)). ¯. The strong maximum principle implies that u(x, t), v(x, t) > 0 when t > 0 for all x ∈ Ω Proof of (b): By the first equation of model (1.1), we have ∂u u − d1 ∆u = ru(1 − ) − f (u, v)v ∂t K u ≤ r(1 − ). K

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

271

According to the comparison principle of parabolic equations, we obtain lim supt→+∞ maxΩ¯ u(x, t) ≤ K. By the maximum principle, we obtain ∥u(·, t)∥C(Ω) ¯ ≤ K1 for all t ≥ 0. ∫ ∫ For the estimate of v(x, t), let U (t) = Ω u(x, t)dx, V (t) = Ω v(x, t)dx; then, ∫ ∫ ∫ ∫ u dU u = ut dx = d1 [ru(1 − ) − f (u, v)v]dx, ∆udx + [ru(1 − ) − f (u, v)v]dx = dt K K Ω Ω Ω Ω ∫ ∫ ∫ ∫ dV = vt dx = d2 ∆vdx + (ef (u, v)v − mv)dx = (ef (u, v)v − mv)dx. dt Ω Ω Ω Ω Multiplying Eq. (2.4) by e and adding the result to Eq. (2.5), we have ∫ u (eU + V )t = −mV + er u(1 − )dx K Ω

(2.4) (2.5)

≤ −m(eU + V ) + e(m + r)U. By using lim supt→+∞ maxΩ¯ u(x, t) ≤ K as proven above, we have lim supt→+∞ maxΩ¯ U (t) ≤ |Ω |. Thus, for small ε > 0, there exists t1 > 0 such that (eU + V )t ≤ −m(eU + V ) + e(m + r)(K + ε)|Ω |, t > t1 . The integration of the inequality leads to ∫ e(m + r) (K + ε)|Ω |, t > t1 , v(x, t)dx = V (t) < eU (t) + V (t) ≤ m Ω which means that

(2.6)

e(m + r) (K + ε)|Ω | := K2 . m ≤ K3 , where K3 depends on K2 and ∥ψ(x)∥L∞ (Ω) . Therefore, there exists

∥v(·, t)∥L1 (Ω) ≤

By Lemma 2.1, we have ∥v(·, t)∥L∞ (Ω) a C ∗ such that ∥v(·, t)∥L∞ (Ω) ≤ C ∗ . Proof of (C): If d1 = d2 , multiplying the first equation by e and adding it to the second equation in (1.1), we have ⎧ ∂w u ⎪ ⎪ = d1 ∆w + eru(1 − ) − mv (x, t) ∈ Ω × (t1 , +∞), ⎪ ⎪ ∂t K ⎨ ∂w (2.7) = 0, x ∈ ∂Ω , t > t1 , ⎪ ⎪ ∂n ⎪ ⎪ ⎩ w(x, 0) = eu(x, t1 ) + v(x, t1 ), x ∈ Ω, where w(x, t) = eu(x, t) + v(x, t). Since eru(1 −

u er ) − mv = e(r + m)u − u2 − mw ≤ e(r + m)(K + ε) − mw K K

in (t1 , +∞) × Ω , the comparison argument shows that lim supt→+∞ maxΩ¯ v(x, t) ≤ lim supt→+∞ maxΩ¯ w(x, t) ≤ Ke(m+r) , which implies the last part of (c). □ m The following theorem provides some sufficient conditions for the permanence of (1.1). Theorem 2.2. Based on Theorem 2.1, assume further that (e − mh)λ(0)K > m holds; then, system (1.1) is permanent with φ(x) ̸≡ 0, ψ(x) ̸≡ 0.

(2.8)

272

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Proof . Let σ1 , σ2 and ϕ1 , ϕ2 be the largest eigenvalue and the corresponding eigenfunctions, respectively, of ⎧ ⎨d1 ∆ϕ + rϕ = σϕ x ∈ Ω , (2.9) ⎩ ∂ϕ = 0 x ∈ ∂Ω , ∂n and ⎧ ⎪ ⎪ ⎨d2 ∆ϕ + (−m +

eλ(0)K )ϕ = σϕ 1 + hλ(0)K

⎪ ∂ϕ ⎪ ⎩ =0 ∂n

x ∈ Ω, (2.10) x ∈ ∂Ω ,

Clearly, σ1 = r > 0, and if (e − mh)λ(0)K > m, then σ2 = −m +

eλ(0)K 1+hλ(0)K ¯ 2

> 0. Therefore, Theorem 5.7 and

Theorem 5.8 in [19] show that system (1.1) generates a semiflow on [C(Ω )] that is permanent. □ Remark 2.1. From the viewpoint of biology, Theorem 2.2 implies that the prey and predator always coexist at any time and in any location of the inhabited domain regardless of the diffusion coefficients. Remark 2.2. Theorem 2.2 implies that system (1.1) is permanent, that is, there exists a t0 > 0, m1 and M1 such that m1 ≤ u(x, t), v(x, t), |λ′ (v)|, λ(v) ≤ M1 ,

(2.11)

for all x ∈ Ω and all t > t0 . 3. Stability and bifurcation In this subsection, we discuss the stability and bifurcation of the equilibria of system (1.1). Let 0 = µ0 < µ1 < µ2 < · · · < µi < · · · → ∞ be the eigenvalues of −∆ on Ω under the homogeneous Neumann boundary condition. We define the following: (i) S(µi ) is the space of eigenfunctions corresponding to µi for i = 0, 1, 2, . . .; (ii) Xij := {c · ϕij : c ∈ R2 }, where {ϕij } are the orthonormal basis of S(µi ) for j = 1, 2, . . . , dim[S(µi )]; and } { ¯ )]2 : ∂u = ∂v = 0 , and thus, X = ⨁∞ Xi , where Xi = ⨁dim[S(µi )] Xij . (iii) X := u = (u, v) ∈ [C 1 (Ω j=1 i=1 ∂n ∂n The linearization of system (1.1) at the positive steady state E ∗ = (u∗ , v ∗ ) can be expressed by ut = (D∆ + J)u,

(3.1)

where D = diag(d1 , d2 ), u = (u(x, t), v(x, t))T , ru∗ − − v ∗ u∗ gu J =⎝ K ev ∗ g + eu∗ v ∗ gu ⎛

∗

∗ ∗

⎞

−u g − u v gv ⎠ , eu∗ v ∗ gv

λ(v) and g(u, v) = 1+hλ(v)u . We can induce the eigenvalues of system (3.1) confined on the subspace Xi . If σ is an eigenvalue of system (3.1) on Xi , then this term must be an eigenvalue of the matrix −µi D + J for each i ≥ 0. Clearly, λ needs to satisfy the characteristic equation

σ 2 − Ti σ + Di = 0,

(3.2)

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

where

273

Ti = −(d1 + d2 )µi + α2 − α1 , Di = d1 d2 µ2i − (d1 α2 − d2 α1 )µi + α3 , ru∗ + v ∗ u∗ gu , α2 = eu∗ v ∗ gv , K e − hm 2C2 λ′ (v ∗ ) ) 2 ∗ + ). α3 = mgv ∗ (r(1 − ∗ λ(v ) λ (v ) e α1 =

(3.3)

E0 is always unstable; thus, prey and predators cannot simultaneously go extinct. For boundary equilibrium E1 , we have the following conclusion. Theorem 3.1. If e < mh, then E1 = (K, 0) is globally asymptotically stable. Proof . For E1 = (K, 0), the characteristic equation (3.12) becomes (σ + d1 µi + r)(σ + d2 µi + m − ef (k, 0)) = 0. Its eigenvalues are σ1i = −r − d1 µi , σ2i = ef (k, 0) − m − d2 µi , where f (K, 0) =

Kλ(0) 1+λ(0)hK .

If e < mh, then ef (K, 0) − m =

(e − mh)Kλ(0) − m < 0. 1 + hλ(0)K

Therefore, σ1i and σ2i are all negative. Thus, E1 is locally asymptotically stable. Indeed, E1 is globally asymptotically stable. On account of Theorem 2.1, we have lim supt→+∞ maxΩ¯ u(·, t) ≤ K, and thus, there exists T1 ∈ (0, +∞) such that for an arbitrary constant ε > 0, u(·, t) ≤ K + ε, t ≥ T1 . It follows from the second equation of (1.1) that vt − d2 ∆v ≤ v(

e − m), t ≥ T1 . h

e < hm, so lim supt→+∞ maxΩ¯ v(·, t) ≤ 0, and there exists T2 > T1 such that v(·, t) ≤ ε, t ≥ T2 . It follows from the first equation of (1.1) that ut − d1 ∆u ≥ ru(1 −

u 1 ) − ε, t > T2 , x ∈ Ω . K h

On account of the arbitrariness of ε > 0, we have lim inf t→+∞ minΩ¯ u(·, t) ≥ K. This aspect, combined with lim supt→+∞ maxΩ¯ u(·, t) ≤ K, allows us to derive lim max|u(·, t) − K| = 0. ¯ t→+∞Ω

Hence, E1 is globally asymptotically stable when e < mh. □ Remark 3.1. Theorem 3.1 shows that if the handling time h is above a certain threshold, the predator becomes extinct.

274

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Theorem 3.2. If conditions 1 + K(

C2 λ(v ∗ )(e − hm) λ′ (v ∗ ) − 1) (h + e 2 ∗ ) > 0, ∗ λ(v ) e λ (v )

2C2 λ′ (v ∗ ) e − hm r(1 − ) + > 0, λ(v ∗ ) λ2 (v ∗ ) e

(3.4)

and

e − mh gv ∗ mλ′ (v ∗ ) 2rC2 − gv ∗ ) + d1 >0 λ e λ2 (v ∗ ) hold, then the positive steady state (u∗ , v ∗ ) is asymptotically stable. d2 (r −

Proof . Denote L=

( d1 ∆ − α1 α5

α4 d2 ∆ + α2

(3.5)

) ,

where α4 = −u∗ g − u∗ v ∗ gv , α5 = ev ∗ g + eu∗ v ∗ gu . Then, for each i (i = 0, 1, 2, . . .), Xi is invariant under the operator L, and ξ is an eigenvalue of L on Xi if and only if ξ is an eigenvalue of the matrix ( ) d1 µi − α1 α4 Ai = . α5 d2 µi + α2 The direct calculation shows that if conditions (3.4) hold, then TrAi = −(d1 + d2 )µi + α2 − α1 = −(d1 + d2 )µi −

C2 λ(v ∗ )(e − hm) λ′ (v ∗ ) ru∗ (1 + K( − 1) (h + e 2 ∗ )) < 0, ∗ K λ(v ) e λ (v )

detAi = d1 d2 µ2i − (d2 α1 − d1 α2 )µi − α1 α2 − α4 α5 = d1 d2 µ2i + (d2 (r − + mgv ∗ (r(1 −

2rC2 e − mh gv ∗ mλ′ (v ∗ ) − gv ∗ ) + d1 )µi λ e λ2 (v ∗ )

2C2 λ′ (v ∗ ) e − hm ) 2 ∗ + ) > 0. ∗ λ(v ) λ (v ) e

Thus, the two roots ξi± of Ai have negative real parts. Moreover, for any i ≥ 0, the following relations hold: (i) If (TrAi )2 − 4detAi ≤ 0, then Reξi± =

1 1 TrAi = (−(d1 + d2 )µi + α2 − α1 ) < 0. 2 2

(ii) If (TrAi )2 − 4detAi > 0, then √ 1 1 1 Reξi− = (TrAi − (TrAi )2 − 4detAi ) ≤ TrAi = (−(d1 + d2 )µk + α2 − α1 ) < 0, 2 2 2 √ 2detAi 1 detAi ˆ √ < −δ, Reξi+ = (TrAi + (TrAi )2 − 4detAi ) = ≤ 2 2 TrAi TrAi − (TrAi ) − 4detAi for some δˆ > 0 that is independent of i. ˆ which does not depend on i, such that It is shown that there exists a positive constant δ, ˆ ∀i ≥ 0. Reξi± < −δ, ˆ According to Theorem 5.1.1 of [20], the conclusion follows, and Therefore, the spectrum of L lies in {ξ|Reξ < −δ}. thus, the proof is completed. □

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

275

3.1. Global stability Theorem 3.3. If ∗

r K ∗

+ ∗

v ∗ λ(v ∗ )hm1 (1+hM12 )2

−

KC2 M1 m2 1

> 0 and

m e

− M1 hKC2 > 0 hold, then the unique positive constant

steady state E = (u , v ) is globally asymptotically stable. Proof . Define

∫

∫ 1 u v )dx + (3.6) (v − v ∗ − v ∗ ln ∗ )dx. ∗ u e v Ω Ω Calculating the derivatives along the positive solution of system (1.1) yields ∫ ∫ dV u∗ 1 ∂v v∗ ∂u = (1 − )dx + (1 − )dx dt u e Ω ∂t v Ω ∂t ∫ ∫ ∫ λ(v)v λ(v ∗ )v ∗ m r (u − u∗ )2 dx − (u − u∗ )( (v − v ∗ )2 dx − )dx − = −I(t) − K Ω 1 + hλ(v)u 1 + hλ(v ∗ )u∗ e Ω Ω ∫ λ(v)u λ(v ∗ )u∗ + (v − v ∗ )( − )dx 1 + hλ(v)u 1 + hλ(v ∗ )u∗ Ω ∫ ∫ ∫ m λ(v) r ∗ 2 ∗ 2 (u − u ) dx − (v − v ) dx − (u − u∗ )(v − v ∗ )dx = −I(t) − K Ω e Ω 1 + hλ(v)u Ω ∫ ∫ v ∗ (λ(v) − λ(v ∗ )) v ∗ λ(v ∗ )hλ(v) ∗ − (u − u )dx − (u − u∗ )2 dx ∗ ∗ ∗ ∗ Ω (1 + hλ(v)u)(1 + hλ(v )u ) Ω (1 + hλ(v)u)(1 + hλ(v )u ) ∫ ∫ λ(v) λ(v ∗ )hλ(v) ∗ ∗ + (u − u∗ )(v − v ∗ )dx, (u − u )(v − v )dx + ∗ )u∗ ) 1 + hλ(v)u (1 + hλ(v)u)(1 + hλ(v Ω Ω ∫ |∇u|2 ∫ |∇v|2 ∗ ∗ where I(t) = −d1 u Ω u2 dx − d2 v Ω v2 dx. In the above derivation, note that (2.11). We obtain ∫ ∫ ∫ dV r m KC2 M1 ≤ −I(t) − (u − u∗ )2 dx − (v − v ∗ )2 dx + (u − u∗ )2 dx dt K Ω e Ω m21 Ω ∫ ∫ ′ λ (ξ)hKC v ∗ λ(v ∗ )hm1 2 (u − u∗ )2 dx − (v − v ∗ )2 dx − ∗ ∗ (1 + hM12 )2 Ω Ω (1 + hλ(v)u)(1 + hλ(v )u ) ∫ ∫ r v ∗ λ(v ∗ )hm1 KC2 M1 m ∗ 2 ≤ −I(t) − ( + − ) (u − u ) dx − ( − M hKC ) (v − v ∗ )2 dx. 1 2 K (1 + hM12 )2 m21 e Ω Ω V (t) =

(u − u∗ − u∗ ln

Under the assumption, we have dVdt(t) ≤ 0, which implies the desired assertion since the equality holds only when (u, v) = (u∗ , v ∗ ). The proof is completed. □ Remark 3.2. Theorem 3.3 implies that under these conditions, the prey and the predator coexist and are evenly distributed in space. 3.2. Hopf bifurcation In this section, we discuss the Hopf bifurcations of system (1.1). Let the parameters r, K, h, e, m and d1 be fixed, and consider d2 > 0 as the bifurcation parameters. Then, according to [21], if there exists an integer i0 ≥ 0 and d∗2 > 0 such that Ti0 (d∗2 ) = 0, Di0 (d∗2 ) > 0, Ti0 ̸= 0, Di0 (d∗2 ) ̸= 0, for all

i ̸= i0 .

(3.7)

and for the unique pair of complex eigenvalues α(d2 ) ± ω(d2 ) of (3.12) near the imaginary axis, then the derivative α′ (d∗2 ) ̸= 0. Then, a Hopf bifurcation occurs at (d∗2 , E ∗ ), which implies that close to the critical magnitude d∗2 , a family of spatially homogeneous or inhomogeneous periodic solutions of the system (1.1) emanates from the constant steady state E ∗ .

276

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Theorem 3.4. (i) Assume that 1 + K( r(1 −

λ(v ∗ )(e − hm) λ′ (v ∗ ) C2 − 1) (h + e ) = 0, λ(v ∗ ) e λ2 (v ∗ )

2C2 λ′ (v ∗ ) e − hm ) + > 0, λ(v ∗ ) λ2 (v ∗ ) e

d2 (r −

(3.8)

e − mh gv ∗ mλ′ (v ∗ ) 2rC2 − gv ∗ ) + d1 >0 λ e λ2 (v ∗ )

hold; then, system (1.1) undergoes a spatially homogeneous Hopf bifurcation. (ii) If d1 µ1 < α2 − α1 , let k0 be the largest positive integer such that α2 − α1 − d1 µk0 > 0. In addition, we assume that d2k1 ̸= d2k2 whenever k1 ̸= k2 , 1 ≤ k1 , k2 < k0 , and { } { } √ √ α2 − α1 α1 max −α1 − α1 α2 + α3 , < µ1 < min −α1 + α1 α2 + α3 , , (3.9) 2d1 d1 ∗ Then, for 1 ≤ k ≤ k0 , system (1.1) undergoes a spatially inhomogeneous Hopf bifurcation occurring at (dH 2k , E ) for 1 ≤ k ≤ k0 , where α2 − α1 − d1 µk dH . 2k = µk

Proof . Proof of (i): If the conditions (3.8) hold, then T0 = 0, Ti ̸= 0(i ≥ 1), D0 = 0, and Di = 0(i ≥ 1). Therefore, Eq. (3.12) has the conjugate complex eigenvalue } 1{ √ ± D0 i , λ= 2 which implies that the system undergoes a spatially homogeneous Hopf bifurcation. H Proof of (ii): By the assumption, we have T1 (dH 21 ) = 0, Tk (d21 ) ̸= 0 for k ≥ 0 and T0 (d20 ) = α2 − α1 > 0. In addition, D0 (d20 ) = α3 > 0 for any d2 > 0. Clearly, 2 2 D1 (dH 21 ) = −d1 µ1 − 2d1 α1 µ1 − α1 + α1 α2 + α3 .

If condition (3.9) holds, then D1 (dH 21 ) > 0. Moreover, if µ1 >

α2 −α1 2d1 ,

then

dDk = 2d1 d2 − (d1 α2 − d2 α1 ) > 0. µk H H Therefore, Dk (dH 21 ) is nondecreasing with respect to k ≥ 1. Hence, when k ≥ 2, Dk (d21 ) ≥ D1 (d21 ) > 0. Therefore, Eq. (3.12) has the conjugate complex eigenvalue } { √ 1 −T1 (d2 ) ± T12 (d2 ) − 4D12 (d2 ) . λ= 2

Clearly, Re’(λ) = − µ21 ̸= 0. ∗ As a consequence, a Hopf bifurcation occurs at (dH 21 , E ), which also implies that the system (1.1) has a family of ∗ inhomogeneous periodic solutions near E . □ 3.3. Turing instability It is well known that the Turing instability occurs only if the prey population diffuses more slowly than the predator population. In the following, we provide the necessary condition for the emergence of a Turing instability.

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

277

First, we consider the spatially homogeneous system corresponding to model (1.1), which is as follows: ⎧ u λ(v)uv du ⎪ ⎪ = ru(1 − ) − , ⎨ dt K 1 + hλ(v)u (3.10) dv eλ(v)uv ⎪ ⎪ ⎩ = − mv. dt 1 + hλ(v)u According to the work by Turing, a positive steady state E ∗ is a Turing instability, meaning that this state is an asymptotically stable steady state of model (3.10) but is unstable with respect to the solutions of spatial model (1.1). Given this result, we can obtain the following results. Theorem 3.5. If the conditions (3.4) and 2rC2 gv ∗ mλ′ (v ∗ ) e − mh − gv ∗ ) + d1 < 0, λ e λ2 (v ∗ ) √ 2α3 − α1 α2 − (α1 α2 − 2α3 )2 − α22 α12 d1 < 0< d2 α22 d2 (r −

(3.11)

hold, then the Turing instability occurs. Proof . According to the discussions above, the characteristic equation of model (3.10) at the positive steady state is as follows: σ 2 + (α1 − α2 )σ + α3 = 0. (3.12) The direct computations shows that if conditions (3.4) hold, then α1 − α2 =

C2 λ(v ∗ )(e − hm) λ′ (v ∗ ) ru∗ (1 + K( − 1) (h + e 2 ∗ )) > 0, ∗ K λ(v ) e λ (v )

2C2 λ′ (v ∗ ) e − hm + ) ) > 0. α3 = mgv (r(1 − λ(v ∗ ) λ2 (v ∗ ) e

(3.13)

∗

Therefore, the positive steady state of model (3.10) is asymptotically stable if condition (3.4) holds. For spatial model (1.1), if conditions (3.4) hold, then Ti > 0. Therefore, (3.12) has no imaginary root with positive real parts; thus, E ∗ is unstable if (3.12) has at least a positive real root. That is, Di < 0 has two real roots, one of which is positive and the other is negative. Note that Di = d1 d2 µ2i − (d1 α2 − d2 α1 )µi + α3 . Therefore, Di takes the minimum value at µi = µmin = d1 α2 − d2 α1 = −(d2 (r − This condition implies that

d1 α2 −d2 α1 2d1 d2

(3.14)

when

2rC2 e − mh gv ∗ mλ′ (v ∗ ) − gv ∗ ) + d1 ) > 0. λ e λ2 (v ∗ )

α22 d21 2α1 α2 − 4α3 d1 + + 1 > 0. α12 d22 α12 d2

(3.15)

(3.16)

Hence, Di is negative when (3.16) is satisfied and for the wave numbers close to µmin , the growth rate of perturbation σ can be positive. Thus, (3.16) is the criterion for Turing instability of (1.1). From (3.16), we have √ 2α3 − α1 α2 − (α1 α2 − 2α3 )2 − α22 α12 d1 0< < (3.17) d2 α22 which completes the proof. □

278

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

4. Nonconstant positive steady states In this section, we analytically derive the condition for the nonexistence and existence of a nonconstant steady state. To this end, we consider the following elliptic system corresponding to the system (1.1): ⎧ λ(v)u u ⎪ ⎪ v, x ∈ Ω , −d1 ∆u = ru(1 − ) − ⎪ ⎪ K 1 + hλ(v)u ⎪ ⎪ ⎨ λ(v)u (4.1) −d2 ∆v = e v − mv, x ∈ Ω, ⎪ 1 + hλ(v)u ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂u = ∂v = 0, x ∈ ∂Ω . ∂ν ∂ν Throughout the remaining part of this paper, the solutions refer to the classical solutions, by which we ¯ ). First, we discuss a priori upper and lower bounds for positive mean the solutions being in C 2 (Ω ) ∩ C 1 (Ω solutions of system (4.1). 4.1. A priori estimates To present our results, we need the following technical lemma [22]. ¯ × R). Lemma 4.1 (Maximum Principle). Suppose that g ∈ C(Ω 2 1 ¯ (i) Assume that ω ∈ C (Ω ) ∩ C (Ω ) satisfies ∆ω(x) + g(x, ω(x)) ≥ 0 in Ω , and

∂ω ≤ 0 on ∂Ω . ∂n

If ω(x0 ) = maxΩ¯ ω, then g(x0 , ω(x0 )) ≥ 0. ¯ ) satisfies (ii) Assume that ω ∈ C 2 (Ω ) ∩ C 1 (Ω ∆ω(x) + g(x, ω(x)) ≤ 0 in Ω , and

∂ω ≥ 0 on ∂Ω . ∂n

If ω(x0 ) = minΩ¯ ω, then g(x0 , ω(x0 )) ≤ 0. ¯ ) be a positive solution Lemma 4.2 (Harnack Inequality). Assume that c(x) ∈ C(Ω ), and let ω ∈ C 2 (Ω ) ∩ C 1 (Ω of ∆ω(x) + c(x)ω(x) = 0, x ∈ Ω , ∂η ω = 0, x ∈ ∂Ω Then, there exists a positive constant C∗ = C∗ (∥c(x)∥α , Ω ) such that max ω(x) ≤ C∗ min ω(x) ¯ Ω

¯ Ω

For convenience, we use the notation ∧ = ∧(r, K, m, h, e). Theorem 4.1 (Upper Bounds). For any positive solution (u, v) of system (4.1), there holds 0 < max u(x) ≤ K, 0 < max v(x) ≤ Ω

Ω

e d1 (r + m )K. m d2

(4.2)

Proof . Let (u(x), v(x)) be a nonnegative solution of (1.1). If there exists x0 ∈ Ω such that v(x0 ) = 0, then v(x) ≡ 0 from the strong maximum principle, and u(x) satisfies ⎧ u ⎨−d1 ∆u = ru(1 − ), x ∈ Ω , K ⎩ ∂u = 0, x ∈ ∂Ω . ∂n

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

279

Similarly, if u(x0 ) = 0 for some x0 ∈ Ω , we also have u(x) ≡ 0, which also implies that v ≡ 0. Otherwise, u(x) > 0 ¯. and v(x) > 0 for x ∈ Ω ¯ . Multiplying the From Lemma 4.1, u(x) ≤ K, and from the strong maximum principle, u(x) < K for all x ∈ Ω first equation of (4.1) by e and adding it to the second equation, we have u ) − mv K emd1 m ≤ erK + K − (ed1 u + d2 v). d2 d2

−(ed1 ∆u + d2 ∆v) = eru(1 −

Then, the maximum principle implies that ed1 u + d2 v ≤ Hence, v(x) <

e m (r

e (d2 r + md1 )K. m

+ m dd12 )K := VK for all x ∈ Ω . □

Theorem 4.2 (Lower Bounds). Assume that m + (mh − e)Kλ(0) > 0. There then exists a positive constant C = C(Λ) such that any solutions (u, v) of system (4.1) satisfy min u(x) ≥ C, min v(x) ≥ C. Ω

Proof . Let c1 (x) = d−1 1 (r(1 −

(4.3)

Ω

u λ(v)v eλ(v)u )− ), c2 (x) = d−1 − m). 2 ( K 1 + hλ(v)u 1 + hλ(v)u

ˆ Λ) such that ∥c1 ∥∞ , ∥c2 ∥∞ ≤ C(d, ˆ Λ) provided that d1 , Then, in view of (4.2), there exists a positive constant C(d, d2 ≥ d. Thus, as u and v satisfy ⎧ ∆u + c1 (x)u = 0, x ∈ Ω , ⎪ ⎪ ⎨ ∆v + c2 (x)v = 0, x ∈ Ω , (4.4) ⎪ ∂u ∂v ⎪ ⎩ = = 0, x ∈ ∂Ω . ∂ν ∂ν the Harnack inequality in Lemma 4.2 shows that there exists a positive constant C2 = C2 (Λ, d) such that max u(x) ≤ C2 min u(x), ¯ Ω

¯ Ω

max v(x) ≤ C2 min v(x). ¯ Ω

¯ Ω

In contrast, suppose that the result is false. There then exists a sequence {(d1,n , d2,n )}∞ n=1 with d1,n , d2,n ≥ d such that the corresponding positive solutions {(un , vn )} of system (4.1) satisfy max un → 0 or max vn → 0 ¯ Ω

¯ Ω

as n → +∞.

(4.5)

By the regularity theory for elliptic equations, there exists a subsequence of {(un , vn )}, which is denoted again by ¯ ) as n → +∞. Observe that u0 ≤ K, and from (4.5), either {(un , vn )} such that {(un , vn )} → (u0 , v0 ) in C 2 (Ω u0 ≡ 0 or v0 ≡ 0. Therefore, we have the following two cases: (i) u0 ≡ 0, v0 ̸≡ 0; or u0 ≡ 0, v0 ≡ 0; (ii) u0 ̸≡ 0, v0 ≡ 0. Since {(un , vn )} is a positive solution of (4.1), one can obtain the following integral equation by integrating Eq. (4.1) for un and vn over Ω , we obtain that ⎧∫ un ⎪ ⎪ ⎨ (r1 un (1 − K ) − f (un , vn )vn )dx = 0, ∫Ω (4.6) ⎪ ⎪ ⎩ (ef (un , vn )vn − mvn )dx = 0. Ω

280

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

(i) In this case, u0 ≡ 0, then ef (un , vn ) − m → −m < 0 uniformly as n → ∞ and vn > 0. For sufficiently large n, we have ∫ (ef (un , vn )vn − mvn )dx < 0, Ω

which is a contradiction. (ii) If u0 ̸≡ 0, v0 ≡ 0, then this result implies that u0 satisfies (4.1). Thus, u0 ≡ K for large n, and ef (un , vn ) − m → ef (K, 0) − m > 0 for large n since f (K, 0) < m e for m + (mh − e)Kλ(0) > 0, which produces a contradiction again for the second integral equation of (4.6), which completes the proof. □ Remark 4.1. By Theorems 4.1 and 4.2, u(x) and v(x) are bounded on Ω . Therefore, there exists m2 and M2 such that m2 ≤ λ(v), |λ′ (v)| ≤ M2 , x ∈ Ω . (4.7) 4.2. Nonexistence of nonconstant positive steady states The following theorem gives the conditions for the nonexistence of nonconstant positive solutions to (4.1). Theorem 4.3. (a) If e < mh, then E1 = (K, 0) is the only nonzero solution of system (4.1); (b) For any fixed r, K, m and e, there exists a positive constant d∗ such that if min{d1 , d2 } > d∗ , then ∗ E = (u∗ , v ∗ ) is the unique positive solution of system (4.1). Proof . Proof of (a): Since e < mh, note that system (1.1) has no positive constant solution in this case. Suppose that (1.1) has a nonconstant positive solution (u, v); then, by Theorem 4.1, we obtain that u(x) ≤ K in Ω , and therefore, we have e eλ(v)u ) ≤ v(−m + ) ≤ 0, x ∈ Ω −d2 ∆v = v(−m + 1 + hλ(v)u h from the assumption. Hence, the maximum principle yields v ≡ 0 in Ω , which is a contradiction. Proof of (b): Let u(x) and v(x) be nonnegative solutions of (4.1). Denote ∫ ∫ −1 −1 u = |Ω | u(x)dx ≥ 0, v = |Ω | v(x)dx ≥ 0. Ω

Then, we have

∫

Ω

∫ (v − v)dx = 0.

(u − u)dx = 0 and Ω

Ω

Multiplying both sides of the first equation of (4.1) by u − u ¯ and using Theorem 4.1, we obtain ∫ ∫ r 2 d1 |∇(u − u)| dx = (u − u)(r(u − u) − (u − u)(u + u) − (f (u, v)v − f (¯ u, v¯)¯ v ))dx K Ω Ω ∫ ∫ λ(v)u ≤3r (u − u)2 dx − (u − u)( (v − v¯) 1 + hλ(v)u Ω Ω λ(v)(u − u ¯) + u ¯λ′ (ξ)(v − v¯) )dx (1 + hλ(v)u)(1 + hλ(¯ v )¯ u) ∫ ∫ 1 KVk M2 1 KVK M2 2 ≤(3r + M2 VK + + ) (u − u) dx + ( + ) (v − v)2 dx 2h 2 2h 2 Ω Ω (4.8) where ξ ∈ (v, v¯) or (¯ v , v). + v¯

281

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Similarly, multiplying both sides of the second equation of (4.1) by v − v¯ and applying Theorem 4.1 gives ∫ ∫ 2 u, v¯)¯ v )dx d2 |∇(v − v)| dx = (v − v)(−m(v − v) + f (u, v)v − f (¯ Ω

Ω

≤(m + ≤(m +

1 + KVK M2 ) h

∫

(v − v)2 dx + M

∫ ∥u − u∥∥v − v∥dx,

(4.9)

Ω

Ω

1 M2 V K + KVk M2 + ) h 2

∫

(v − v)2 dx +

Ω

VK M 5 2

∫

(u − u)2 dx.

Ω

Using the P oincar´ e inequality ∫ ∫ ∫ ∫ 2 2 µ1 (u − u)2 dx ≤ |∇(u − u)| dx, µ1 (v − v)2 dx ≤ |∇(v − v)| dx. Ω

Ω

Ω

Ω

Combining (4.8) and (4.9) leads to ∫

∫

2

Ω

Ω

∫ ≤A

2

A = 3r +

∫

(u − u) dx + B Ω

where

2

(v − v) dx

(u − u) dx + d2 µ1

d1 µ1

(4.10) 2

(v − v¯) dx, Ω

3M2 VK 1 KVK M2 3 3KVK M2 M2 V K + + ,B = m + + + . 2 2h 2 2h 2 2

This result shows that if min{d1 , d2 } >

1 max{A, B}, µ1

then ∇(u − u) = ∇(v − v) = 0, and (u, v) must be a constant solution. □ 4.3. Existence of nonconstant positive steady states In this subsection, we discuss the existence of nonconstant positive solutions to system (4.1) when the diffusion coefficients d1 and d2 vary while the parameters r, K, m, e, h, b, and ω are fixed by using the Leray–Schauder degree theory. For simplicity, denote u = (u, v), and let E ∗ = (u∗ , v ∗ ) be the possible positive constant solutions of system (4.1). ⎛ ⎞ u λ(v)uv ( ) ru(1 − ) − ru∗ ∗ ∗ ∗ ∗ ∗ ⎜ ⎟ − − v u g −u g − u v g K 1 + hλ(v)u u v K ⎟,B = Φ(u) = ⎜ , D = Diag{d1 , d2 }, ⎝ ⎠ eλ(v)uv ev ∗ g + eu∗ v ∗ gu eu∗ v ∗ gv −mv + 1 + hλ(v)u { } ¯ )]2 : ∂u = ∂v = 0 on ∂Ω , X + = {(u, v) : u, v ≥ 0, (u, v) ∈ X}, X = (u, v) ∈ [C 1 (Ω ∂n ∂n ¯ }. Γ = {(u, v) ∈ X : C −1 < u, v < C on Ω Thus, (4.1) can be written as ⎧ ⎨−∆u = D−1 Φ(u), x ∈ Ω , ⎩ ∂u = 0, x ∈ Ω. ∂n

(4.11)

282

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Therefore, u is a positive solution of (4.11) if and only if it satisfies F(u) := u − (I − ∆)−1 (D−1 Φ(u) + u)

in X + ,

(4.12)

where (I − ∆)−1 is the inverse of I − ∆ with the homogeneous Neumann boundary condition. As F(·) is a compact perturbation of the identity operator, the Leray–Schauder degree deg(F(·), Γ , 0) is well-defined because of F(u) ̸= 0 on ∂Γ . By direct computation, we have Fu (E ∗ ) = I − (I − ∆)−1 (D−1 B + I). If Fu (E ∗ ) is invertible, the index of F at E ∗ is defined as index(F(·), E ∗ ) = (−1)γ , where γ is the number of negative eigenvalues of Fu (E ∗ ). Note that ξ is an eigenvalue of Fu (E ∗ ) on Xj if and only if ξ(1 + µj ) is an eigenvalue of the matrix Mj := µj I − D−1 B. Thus, Fu (E ∗ ) is invertible if and only if for all j ≥ 0, the matrix Bj is nonsingular. Denoting H(d1 , d2 ; µ) = d1 d2 det[µI − D−1 B] = det[µD − B],

(4.13)

we have that if H(d1 , d2 ; µ) ̸= 0; then, H(d1 , d2 ; µ) < 0 if and only if the number of negative eigenvalues of Fu (E ∗ ) in Xj is odd. The following lemma provides the explicit formula of calculating the index. Lemma 4.3. If H(d1 , d2 ; µj ) ̸= 0 for all j ≥ 0, then ∑

index(F(·), E ∗ ) = (−1)γ , γ =

m(µj ),

j≥0,H(d1 ,d2 ;µ)<0

where m(µj ) is the algebraic multiplicity of µj . To facilitate our computation of deg(F(·), E ∗ ), we need only consider the sign of det[µD −Φu (E ∗ )]. Direct calculation gives H(d1 , d2 ; µ) = det[µD − B] = d1 d2 µ2 − (α2 d1 − α1 d2 )µ + α3 = 0,

(4.14)

where α1 , α2 and α3 are defined in (3.3). Simple computations indicate that if Υ := (α2 d1 − α1 d2 )2 − 4d1 d2 α3 > 0,

(4.15)

then Eq. (4.14) has exactly two different positive roots µ+ and µ− , which are, respectively, expressed as √ 1 µ− = (α2 d1 − α1 d2 − Υ ), 2d1 d2 √ 1 µ+ = (α2 d1 − α1 d2 + Υ ). 2d1 d2 Here, we remark that for any fixed d1 > 0, (4.15) must be true for large d2 . Set Sp = {µ0 , µ1 , µ2 , . . .} and let m(µj ) be the multiplicity of µj . We find that lim µ− = 0, lim µ+ =

d2 →∞

d2 →∞

α1 . d1

(4.16)

Clearly, H(d1 , d2 ; µ) < 0 if and only if µ ∈ (µ− , µ+ ). We have the following main result of this subsection.

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

283

Theorem 4.4. Assume that Υ > 0, and let µ− < µ+ be the two positive solutions of Φu (E ∗ ). If there exist i, j ∈ N such that 0 ≤ µi < µ− < µi+1 ≤ µj < µ+ < µj+1 ∑j and k=i+1 m(µk ) is odd, then (4.1) has at least one nonconstant positive solution. Proof . For t ∈ [0, 1], we define (1 − t

⎞ t ) f (u, v) ⎟ d∗ d1 ⎜ At (u) ≜ (−∆ + I)−1 ⎝ ⎠, (1 − t t ) + g(u, v) v+ d∗ d2 ⎛

u+

+

where d∗ is defined in Theorem 4.3. The positive solutions of the problem At (u) = u

in Ω ,

∂u = 0 on ∂n

∂Ω

(4.17)

are contained in Γ . Note that u is a positive solution of system (4.1) if and only if it is a positive solution of (4.17) with t = 1. u∗ is the unique positive constant solution of (4.17) for any t ∈ [0, 1]. According to the choice of d∗ in Theorem 4.3, we find that E ∗ is the only fixed point of A0 . deg(I − A0 , Γ , 0) = index(I − A0 , Γ , E ∗ ) = 1.

(4.18)

Since F = I − H(·, 1), and if (4.1) has no other solutions except the constant one E ∗ , then we have ∑i m(µk ) deg(I − A1 , Γ , (0, 0)) = index(F, E ∗ ) = (−1) k=j+1 = −1 On the other hand, by the homotopy invariance of the topological degree, deg(I − A0 , Λ, 0) = deg(I − A1 , Λ, 0),

(4.19)

which is a contradiction. Therefore, there exists at least one nonconstant solution of (4.1). □ 5. Numerical results and discussions In this section, we present some results of the numerical simulations to illustrate our mathematical findings in the previous sections. 5.1. Steady-state bifurcation λ0 We choose λ(v) = (b+v) ω and the parameters r = 2, K = 8, e = 1, m = 1, h = 0.25, λ0 = 5, b = 0.5, ω = 1.5, d1 = 0.008 and d2 = 0.2. System (1.1) is in a state of oscillation (see Fig. 1), which is in accord with the conclusion in [23] for the ODE system. However, if we set d2 = 1, then system (1.1) has a stationary spatial pattern, which means that system (1.1) undergoes steady-state bifurcation. (Fig. 2).

5.2. Turing instability λ0 Let λ(v) = (b+v) ω , and take the parameters to be same as those in [23], i.e., r = 2, K = 8, e = 1, m = 1, h = 0.25, λ0 = 5, b = 0.5, and ω = 1.5. Let the diffusion coefficients d1 = 0.008, d2 = 0.2 and Ω = (0, π). By simple calculation, we obtain that system (1.1) has a unique positive steady state E ∗ = (0.2925, 0.5637).

284

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Fig. 1. System (1.1) is in a state of oscillation with d2 = 0.2.

Fig. 2. System (1.1) has a stationary spatial pattern with d2 = 1.

Fig. 3. The positive steady state E ∗ = (0.29250.5637) for the ODE system is locally asymptotically stable.

The positive steady state E ∗ for the ODE system is stable (see Fig. 3). However, for the PDE system, the positive steady state becomes unstable, and system (1.1) has a stable nonconstant equilibrium solution, which means that a Turing instability occurs (see Fig. 4). The implication is that the population is unevenly distributed in space. In addition, we take the values of the parameters to be the same as in [2], i.e., taking r = 1, K = 4, e = 2, m = 4, h = 0.25, λ0 = 6.3, and b = 1. According to the conclusion in [2], the positive constant steady state for the ODE system is unstable with ω < 0.5 and becomes locally asymptotically stable with ω > 0.5. However, in our system, if we set Ω = (0, 100) and the diffusion coefficients d1 = 0.008 and d2 = 1, then we

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

285

Fig. 4. The positive steady state E ∗ = (0.29250.5637) is unstable, and there exists a stable nonconstant equilibrium solution.

Fig. 5. The positive steady state of system (1.1) is unstable, and there exists a stationary spatial pattern.

can find that the positive constant steady state of system (1.1) becomes unstable and that the system (1.1) has a stationary spatial pattern (see Fig. 5) with ω = 0.6). 5.3. Spatial chaotic patterns λ0 Letting λ(v) = (b+v) w and choosing r = 0.7, K = 8, e = 1, m = 1, h = 0.25, λ0 = 5, b = 0.5, ω = 3, d1 = 0.008, d2 = 0.2 and Ω = (0, 60), we find that under the effect of a diffusion rate of d2 , system (1.1) shows different spatial patterns that can be stationary, periodic and chaotic, as shown in Fig. 6. λ0 Choose λ(v) = L − (b+v) w and the parameters r = 3, K = 17, e = 1, m = 1, h = 0.25, λ0 = 5.1783, b = 2, ω = 0.06, and L = 5, where the function λ(v) and the value of the parameters are same as the ones in [3]. According to the discussions in [3], the ODE system has a stable limit cycle. For the PDE system (1.1), take Ω = (0, 100) and the diffusion coefficients d1 = 0.008 and d2 = 3; then, we find that the system (1.1) presents chaotic behavior and has a spatial pattern (see Fig. 7).

6. Concluding remarks In this paper, we have studied a class of predator–prey population model under homogeneous Neumann boundary conditions that covers many existing well-known examples. We have investigated the global boundedness, stability of the constant steady state, and existence and nonexistence of nonconstant steady states in addition to the Hopf and Turing bifurcations.

286

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

Fig. 6. Examples of spatial patterns (snapshots at the exhibited final time) and spatially averaged population dynamics for d2 = 0.2 (left), d2 = 1 (middle), d2 = 4 (right).

Fig. 7. System (1.1) presents a spatial chaotic pattern.

In Section 2, we showed that our system is dissipative and uniformly persistent. One of the most fundamental questions for biological models concerns the long-term behavior (i.e., survival or extinction) of each species. In particular, the long-term survival of all species is equivalent to the persistence analysis of the related mathematical models. Clearly, for a dissipative system, the uniform persistence is equivalent to permanence. By Theorem 2.1, we know that model (1.1) is bounded. Theorem 2.2 shows that if condition (2.8) holds, then system (1.1) is permanent. In Section 3, the local/global stability for all feasible constant steady states was investigated, as were the Hopf and Turing bifurcations. Theorem 3.1 shows that if e < mh, then the equilibrium E1 (K, 0) is globally asymptotically stable. This fact implies that the predator becomes extinct under this condition. By Theorem 3.2, we know that if conditions (3.4) and (3.5) hold, then the equilibrium E ∗ = (u∗ , v ∗ ) is locally

X. Zhang and H. Zhu / Nonlinear Analysis: Real World Applications 48 (2019) 267–287

287

asymptotically stable, which implies that if the initial state nears E ∗ , two species coexist for system (1.1). Theorem 3.3 shows that E ∗ = (u∗ , v ∗ ) is globally asymptotically stable, which implies that the prey and predator coexist. Theorem 3.4 tells us that diffusion can create a Hopf bifurcation. In Section 4, we obtain results on the nonexistence and existence of the nonconstant positive steady-state solution of system (4.1). Theorem 4.3 implies nonexistence of the nonconstant positive steady-state solution when the diffusion coefficients are large. From Theorem 4.4, it follows that the existence of the nonconstant positive steady state depends upon certain parametric restrictions. Acknowledgments We would like to express our gratitude to the referees for their valuable comments and suggestions that led to a truly significant improvement of the manuscript. The work is supported by the Natural Science Foundation of Jiangsu Province, China under grant no. BK20150420 and by the Startup Foundation for Introducing Talent of NUIST. References [1] L. Berec, Impacts of foraging facilitation among predators on predator–prey dynamics, Bull. Math. Biol. 72 (1) (2010) 94–121. [2] L. Pˇribylov´a, L. Berec, Predator interference and stability of predator–prey dynamics, J. Math. Biol. 71 (2) (2015) 301–323. [3] L. Pˇribylov´a, A. Peniaˇskov´a, Foraging facilitation among predators and its impact on the stability of predator–prey dynamics, Ecol. Complex. 29 (2017) 30–39. [4] C. Cosner, D.L. DeAngelis, J.S. Ault, D.B. Olson, Effects of spatial grouping on the functional response of predators, Theoret. Popul. Biol. 56 (1) (1999) 65–75. [5] C.B. Huffaker, Experimental studies predators: disperson factors and predator–prey oscillations, Hilgardia 27 (1958) 343–383. [6] L.S. Luckinbill, Coexistence of laboratory populations of paramecium aurelia and its predator didinium nautum, Ecology 54 (1973) 1320– 1327. [7] R.S. Cantrell, C. Cosner, Spatial Ecology via Reaction–Diffusion Equations, Wiley, London, 2004. [8] E.O. Wilson, R.H. MacArthur, The Theory of Island Biogeography, Princeton: Princeton University Press. [9] L. Yuan, W. Ni, S. Yotsutani, Pattern formation in a cross-diffusion system, Discrete Contin. Dyn. Syst. Ser. A 35 (2015) 1589–1607. [10] Y. Cai, Y. Kang, M. Banerjee, W. Wang, Complex dynamics of a host–parasite model with both horizontal and vertical transmissions in a spatial heterogeneous environment, Nonlinear Anal. RWA 40 (2018) 444–465. [11] Y. Cai, W. Wang, Fish-hook bifurcation branch in a spatial heterogeneous epidemic model with cross-diffusion, Nonlinear Anal. RWA 30 (2016) 99–125. [12] L. Zhang, J. Liu, M. Banerjee, Hopf and steady state bifurcation analysis in a ratio-dependent predator–prey model, Commun. Nonlinear Sci. Numer. Simul. 44 (2017) 52–73. [13] X. Tian, Global dynamics for a diffusive predator–prey model with prey-taxis and classical Lotka-Volterra kinetics, Nonlinear Anal. RWA 39 (2018) (2018) 278–299. [14] W. Ni, M. Wang, Dynamics and patterns of a diffusive Leslie-Gower prey-predator model with strong allee effect in prey, J. Differential Equations 261 (7) (2016) 4244–4274. [15] R. Ruiz-Baier, C. Tian, Mathematical analysis and numerical simulation of pattern formation under cross-diffusion, Nonlinear Anal. RWA 14 (1) (2013) 601–612. [16] J. Wang, J. Wei, J. Shi, Global bifurcation analysis and pattern formation in homogeneous diffusive predator–prey systems, J. Differential Equations 260 (4) (2016) 3495–3523. [17] G. Gambino, M.C. Lombardo, M. Sammartino, Pattern formation driven by cross-diffusion in a 2d domain, Nonlinear Anal. RWA 14 (3) (2013) 1755–1779. [18] R. Peng, F.Q. Yi, X.Q. Zhao, Spatiotemporal patterns in a reaction–diffusion model with the Degn-Harrison reaction scheme, J. Differential Equations 254 (6) (2013) 2465–2498. [19] R.S. Cantrell, C. Cosners, V. Hutson, Permanence in ecological systems with spatial heterogeneity, Proc. Roy. Soc. Edinburgh Sect. A Math. 123 (3) (1993) 533–559. [20] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Springer, 1981. [21] F. Yi, J. Wei, J. Shi, Bifurcation and spatiotemporal patterns in a homogeneous diffusive predator–prey system, J. Differential Equations 246 (5) (2009) 1944–1977. [22] Y. Lou, W.M. Ni, Diffusion, self-diffusion and cross-diffusion, J. Differential Equations 131 (1) (1996) 79–131. [23] L. Berec, Impacts of foraging facilitation among predators on predator–prey dynamics, Bull. Math. Biol. 72 (1) (2010) 94–121.