Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 70, pp. 1–31.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT

IN SPATIALLY HETEROGENEOUS ENVIRONMENT

SHANBING LI, YANNI XIAO, YAYING DONG

Abstract. This article concerns diffusive predator-prey models incorporating

the cost of fear and environmental heterogeneity. Under homogeneous Neu-

mann boundary conditions, we establish the uniform boundedness of global
solutions and global stability of the trivial and semi-trivial solutions for the

parabolic system. For the corresponding steady-state problem, we obtain suf-

ficient conditions for the existence of positive steady states, and then study
the effects of functional responses and the cost of fear on the existence, sta-

bility and number of positive steady states. We also discuss the effects of

spatial heterogeneity and spatial diffusion on the dynamic behavior and estab-
lish asymptotic profiles of positive steady states as the diffusion rate of prey or

predator individuals approaches zero or infinity. Our theoretical results sug-

gest that fear plays a very important role in determining the dynamic behavior
of the models, and it is necessary to revisit existing predator-prey models by

incorporating the cost of fear.

1. Introduction

In 2011, Zanette et al. [30] conducted a manipulation on female song sparrows
(Melospiza melodia) during an entire breeding season. This is the first experimental
evidence demonstrating that the perceived predation risk can affect the populations
of terrestrial vertebrates although many biologists early realized that the cost of
fear should be considered except for direct killing in the predator-prey interactions
(for example, see [5, 16, 19] and references therein). Since then, fear effect in the
predator-prey interactions has attracted considerable attention (for example, see
[7, 18, 21, 24, 25, 26, 27] and references therein).

In 2016, Wang, Zanette and Zou [25] initially proposed and analyzed the follow-
ing ODE model:

du

dt
= r0uf(k, v)− du− au2 − g(u)v,

dv

dt
= −mv + cg(u)v,

(1.1)

which models the fear effect in predator-prey interactions. Here u(t) and v(t) denote
the population densities of respective species at time t > 0; r0 > 0 is the birth rate
of the prey; d > 0 is the natural death rate of the prey; a > 0 is the intra-specific

2010 Mathematics Subject Classification. 35B30, 35B32, 92D25.
Key words and phrases. Predator-prey model; fear cost; spatial diffusion;

environmental heterogeneity.
c©2021. This work is licensed under a CC BY 4.0 license.

Submitted November 16, 2020. Published August 23, 2021.

1



2 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

pressure of the prey; m > 0 is the natural death rate of the predator; c > 0 is
the conversion rate; g(u) represents the predator functional response which is a
continuously differentiable function of u ∈ R+; f(k, v) represents the cost of fear,
where k ≥ 0 accounts for the level of fear, and f(k, v) satisfies

f(0, v) = 1,
∂f(k, v)

∂k
< 0, lim

k→∞
f(k, v) = 0,

f(k, 0) = 1,
∂f(k, v)

∂v
< 0, lim

v→∞
f(k, v) = 0.

The following 3 functions satisfy all the above hypotheses: f(k, v) = 1/(1 + kv),
f(k, v) = e−kv, and f(k, v) = 1/(1 + kv + k′v2). These functions describe different
decreasing rates. By theoretical and numerical analyses, Wang, Zanette and Zou
[25] have demonstrated that high levels of fear can exclude the existence of periodic
solutions, while the low levels of fear can induce multiple limit cycles. These results
quantitatively reveal the effect of the cost of fear on the dynamics of (1.1).

When the spatial distribution of respective species and the intra-specific pressure
of the predator are considered, Wang and Zou [27] proposed and analyzed the
following reaction-diffusion-advection predator-prey model incorporating the cost
of fear and avoidance behaviors of the prey:

∂u

∂t
= du∆u+ α∇ · (β(u)u∇v) + f0(k0α, v)r0u

− du− au2 − up(u, v)v, x ∈ Ω, t > 0,

∂v

∂t
= dv∆v + v(−m(v) + cup(u, v)), x ∈ Ω, t > 0,

∂u

∂µ
=
∂v

∂µ
= 0, x ∈ ∂Ω, t > 0,

u(x, 0) = u0(x) ≥ 0, v(x, 0) = v0(x) ≥ 0, x ∈ Ω.

(1.2)

Here Ω is a bounded domain in Rn (n ≥ 1) with smooth boundary ∂Ω; u(x, t) and
v(x, t) are the population densities of respective species at location x ∈ Ω and time
t > 0; du > 0 and dv > 0 are the random diffusion coefficients of respective species;
m(v) = m1 +m2v, where m1 > 0 is the death rate of the predator and m2 ≥ 0 is the
intra-specific pressure of the predator; up(u, v) is the predator functional response;
µ(x) is the outer unit normal vector at x ∈ ∂Ω and ∂u/∂µ = µ(x) · ∇u is the out-
flux of u; α∇ · (β(u)u∇v) is the diffusion of the prey which represents a directed
movement towards lower density of the predator (i.e., predator-taxis), where α
measures the strength of predator-taxis, and β(u) = 1− u/M for 0 ≤ u ≤ M , = 0
for u > M , where M > 0 measures the maximum number of the prey that a unit
volume can accommodate; f0(k0α, v) = 1/(1 + k0αv) satisfies the same hypotheses
as f(k, v) with k0 as a nonnegative constant. By theoretical and numerical analyses,
Wang and Zou [27] have established the sufficient and necessary conditions of spatial
pattern formation for different functional response, and showed that the cost of fear
and functional responses play an important role in spatial pattern formation.

For most biological species, the natural environment where they live is usually
spatially heterogeneous. Therefore, it is reasonable to expect the dynamic behavior
to be influenced by the environmental heterogeneity, apart from the direct effect
(through predation) and indirect effect (fear effect) between the species. Moreover,



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 3

the impact of environmental change on the dynamics of the predator-prey inter-
actions is increasingly recognized (for example, see [2, 6, 10, 9, 15] and references
therein). Therefore, it seems imperative to include such the environmental het-
erogeneity while modelling the predator-prey interactions incorporating fear effect.
This constitutes our first motivation of the present paper.

Note that in particular that when modeling the predator-prey interactions incor-
porating fear effect in the literatures mentioned above, the predator species is always
assumed to be specialist predators. However, most predator-prey interactions are
generally comprised of a wide variety of predators, both specialists and generalists.
Moreover, ecologists have studied the impacts of generalist versus specialist preda-
tors separately along with the potential outcome of interactions between them in
various environments, and have shown many different dynamic behaviors (for ex-
ample, see [12, 13, 23] and references therein). Therefore, the predator species
should not only be assumed to be specialists, but should include both specialists
and generalists. This constitutes our second motivation of the present paper.

As far as we know, although lots of mathematical models (mostly ODE models)
have been proposed and studied to quantitatively investigate the effect of fear cost
in the predator-prey interactions, there are few works to study a diffusive predator-
prey model incorporating fear cost in spatially heterogeneous environment. There-
fore, based on these considerations, we consider the following reaction-diffusion
predator-prey model incorporating fear effect in spatially heterogeneous environ-
ment:

∂u

∂t
= du∆u+

ru

1 + kv
− du− au2 − b(x)p(u, v)v, x ∈ Ω, t > 0,

∂v

∂t
= dv∆v +mv − v2 + cb(x)p(u, v)v, x ∈ Ω, t > 0,

∂u

∂µ
=
∂v

∂µ
= 0, x ∈ ∂Ω, t > 0,

u(x, 0) = u0(x) ≥ 0, v(x, 0) = v0(x) ≥ 0, x ∈ Ω.

(1.3)

Note that we extend model (1.2) by assuming the predator-prey interaction b(x) to
be spatially dependent function instead of constant, and assuming the growth rate
of the predator to be either positive or negative. In particular, m > 0 means that
predator individuals are generalists and m < 0 means that predator individuals
are specialists. The function b(x) is positive and Hölder continuous in Ω. The
functional response p(u, v) is a non-negative C1-function of (u, v) ∈ [0,∞)× [0,∞)
such that

p(0, v) = 0, 0 ≤ pu(u, v) <∞, puu(u, v) ≤ 0,

pv(u, v) ≤ 0,
∂(p(u, v)v)

∂v
≥ 0

for (u, v) ∈ [0,∞)×[0,∞). The most widely used forms of p(u, v)v in the literatures
are p(u, v)v = uv (Linear functional response), p(u, v)v = uv/(1 + qu) (Holling-type
II functional response) and p(u, v)v = uv/(1 + qu+ fv) (Beddington-DeAngelis
functional response). Here q, f are positive constants. Note in particular that
compared with (1.2), we compromise a little bit in the diffusion term for the prey
by considering the random movement instead of the directed movement towards
lower density of the predator. Indeed, in spatially heterogeneous environment,



4 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

such directed movement for the prey is interesting, biologically important but yet
mathematically challenging, and we have to leave it for future research project.

The main purpose of this paper is to reveal the impact of fear cost, spatial
diffusion and environmental heterogeneity on the dynamics of (1.3). In order to
better demonstrate their influence, we also need to consider the non-negative steady
states of (1.3) which satisfy the following nonlinear elliptic equations:

−du∆u =
ru

1 + kv
− du− au2 − b(x)p(u, v)v, x ∈ Ω,

−dv∆v = mv − v2 + cb(x)p(u, v)v, x ∈ Ω,

∂u

∂µ
=
∂v

∂µ
= 0, x ∈ ∂Ω.

(1.4)

Here u(x) and v(x) denote the density of the prey and predator individuals at equi-
librium, respectively. It is, naturally, the dynamics in the biologically meaningful
region u ≥ 0, v ≥ 0 that are of interest. It is clear that (1.4) admits a trivial
solution (u, v) = (0, 0), two semi-trivial solutions (u, v) = ((r − d)/a, 0) with r > d
and (u, v) = (0,m) with m > 0, and positive solutions (u, v) with no component
identically zero. From now on, (u, v) will be called a positive solution of (1.4) or
a positive steady state of (1.3) if (u, v) is a classical solution satisfying u > 0 and
v > 0 in Ω.

The rest of this paper is organized as follows. In Sect. 2, we establish the long-
time behavior of solutions to (1.3) and the sufficient conditions for the existence of
positive solutions to (1.4). In Sect. 3, we obtain some insights on how fear cost
affects the population dynamics of (1.3) by choosing different functional responses
and different sets of parameters. Section 4 is dedicated to the effect of spatial
diffusion and environmental heterogeneity on the dynamics of (1.3). Finally, the
significance of current studies is outlined in Sect. 5.

2. Solution and equilibria of system (1.3)

The purpose of this section is to investigate the long-time behavior of solutions
to (1.3) and the existence of positive solutions to (1.4).

2.1. Long-time behavior of solutions to (1.3). For any given function φ ∈
C(Ω), we denote

φ∗ =: max
Ω

φ, φ∗ =: min
Ω
φ, φ =

1

|Ω|

∫
Ω

φdx,

where |Ω| represents the volume of the region Ω. Let X = {φ ∈W 2,p
µ (Ω) : ∂φ/∂µ =

0} and Y = Lp(Ω) with p > n. Define a closed linear operator D in Y × Y by
D(u, v) = (−du∆u,−dv∆v), where (u, v) ∈ D(D) =: X ×X. Here D(D) represents
the domain of the operator D. By [29], the closed linear operator −D generates an
analytic semigroup {e(−tD)}t≥0 in Y × Y . Then we can use the similar argument
to that of Proposition 3.1 in [29] to prove the following global solvability theorem.

Theorem 2.1. Assume that the initial values u0(x) and v0(x) are non-negative
functions of class C(Ω). There exists a unique solution (u(x, t), v(x, t)) of (1.3) in
C([0,∞);Y × Y ) ∩ C1((0,∞);D(D)), and it satisfies

0 ≤ u(x, t) ≤ max {(r − d)/a, u∗0} =: M1,

0 ≤ v(x, t) ≤ max {m+ cb∗p(M1, 0), v∗0}



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 5

in Ω× [0,∞).

For later discussion, we now collect some results on linear eigenvalue problem
(see e.g., [1, 3]). For any q(x) ∈ Cν(Ω) with ν ∈ (0, 1), the eigenvalue problem

− d∆φ+ q(x)φ = λφ, x ∈ Ω,
∂φ

∂µ
= 0, x ∈ ∂Ω (2.1)

has an infinite number of eigenvalues. Let them be {λi(d, q(x))}∞i=1 satisfying
λi(d, q(x)) ≥ λj(d, q(x)) for i > j ≥ 1, where λ1(d, q(x)) is the least eigenvalue and
is called the principal eigenvalue. In particular, λ1(d, 0) = 0. Moreover, λ1(d, q(x))
is a simple eigenvalue and the corresponding eigenfunction does not change sign in
Ω. It follows from the variational characterization that λ1(d, q(x)) is given by

λ1(d, q(x)) = inf
φ∈H1(Ω),φ6=0

∫
Ω

(
d|∇φ|2 + q(x)φ2

)
dx
/∫

Ω

φ2dx.

Next we call some properties of λ1(d, q(x)).

Proposition 2.2. The following assertions hold.

(1) λ1(d, q(x)) is continuous and monotone increasing with respect to q(x) in
the sense that q1 ≤ q2 and q1 6≡ q2 implies λ1(d, q1(x)) < λ1(d, q2(x)).

(2) λ1(d, q(x)) is strictly monotone increasing with respect to d > 0 such that
λ1(d, q(x))→ q∗ as d→ 0+ and λ1(d, q(x))→ q as d→∞.

Following [29], we say that any non-negative solution of (1.4) is locally asymp-
totically stable provided that the spectrum of the corresponding linearized operator
lies in the right-hand side of the imaginary axis and it is unstable provided that
there are some points in the spectrum with negative real parts. Then the (local)
stability of trivial and semi-trivial solutions reads as follows.

Theorem 2.3. The following assertions hold.

(1) (0, 0) is locally asymptotically stable if r < d and m < 0; unstable if r > d
or m > 0.

(2) ( r−da , 0) is locally asymptotically stable if λ1

(
dv,−m− cb(x)p

(
r−d
a , 0

))
>

0; unstable if λ1

(
dv,−m− cb(x)p

(
r−d
a , 0

))
< 0.

(3) (0,m) is locally asymptotically stable if λ1

(
du,− r

1+km+d+b(x)pu(0,m)m
)
>

0; unstable if λ1

(
du,− r

1+km + d+ b(x)pu(0,m)m
)
< 0.

Proof. Since the proofs of (1)-(3) are similar, we only prove (3) here. The linearized
operator of (1.4) at (u, v) = (0,m) is

L(0,m) =

(
−du∆− r

1+km + d+ b(x)pu(0,m)m 0

−cb(x)pu(0,m) −dv∆ +m

)
.

It follows from the Riesz-Schauder theory that the spectrum σ(L(0,m)) of the oper-
ator L(0,m) is composed of real eigenvalues, moreover

σ(L(0,m)) = σ
(
− du∆− r

1 + km
+ d+ b(x)pu(0,m)m

)
∪ σ (−dv∆ +m) .

Since λi(dv,m) ≥ λ1(dv,m) = m > 0, there is no point with a negative real part in
the spectrum σ(−dv∆ +m). In addition, the spectrum

σ
(
− du∆− r

1 + km
+ d+ b(x)pu(0,m)m

)



6 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

lies on the real axis and the least eigenvalue is given by

λ1

(
du,−

r

1 + km
+ d+ b(x)pu(0,m)m

)
.

Hence, the stability semi-trivial solution (0,m) is determined by the sign of this
eigenvalue. This completes the proof. �

By applying the comparison principle, the global attractivity of trivial and semi-
trivial solutions can be established.

Theorem 2.4. The following assertions hold.

(1) Assume that m < 0 and r ≤ d. Then any non-negative solution of (1.3)
converges to (0, 0) uniformly in Ω as t→∞.

(2) Assume that m < 0, r > d and λ1 (dv,−m− cb(x)p ((r − d)/a, 0)) > 0.
Then any non-negative solution of (1.3) with u0(x) ≥ ( 6≡)0 converges to
((r − d)/a, 0) uniformly in Ω as t→∞.

(3) Assume that m > 0. Then any non-negative solution of (1.3) with v0(x) ≥
( 6≡)0 converges to (0,m) uniformly in Ω as t → ∞ if one of the following
hypotheses holds:
(i) r

1+km ≤ d;

(ii) r
1+km > d and λ1

(
du,− r

1+km + d+ b(x)pu

((
r

1+km − d
)
/a,m

)
m
)
>

0.

Proof. Since the proofs of (1)–(3) are similar and the proof of (3) is a little more
complicated, we only give prove (3). It follows from equation (1.3) for v(x, t) that

∂v

∂t
= dv∆v +mv − v2 + cb(x)p(u, v)v ≥ dv∆v +mv − v2

in Ω× (0,∞). Let V (x, t) be the solution of

Vt = dv∆V +mV − V 2, x ∈ Ω, t > 0,

∂V

∂µ
= 0, x ∈ ∂Ω, t > 0,

V (x, 0) = v0(x) ≥ (6≡)0, x ∈ Ω.

It is well known that V (x, t)→ 0 uniformly in Ω as t→∞ if m ≤ 0, and V (x, t)→
m uniformly in Ω as t→∞ if m > 0. Thus, we apply the comparison principle to
derive that v(x, t) ≥ V (x, t) in Ω× (0,∞). This ensures that for any ε1 > 0, there
exists Tε1 > 0 such that v(x, t) ≥ m− ε1 in Ω× [Tε1 ,∞) since m > 0. In addition,
it follows from the equation for u(x, t) of (1.3) that

∂u

∂t
= du∆u+

ru

1 + kv
−du−au2−b(x)p(u, v)v ≤ du∆u+

( r

1 + k(m− ε1)
−d
)
u−au2

in Ω× [Tε1 ,∞). As above, we apply the comparison principle to derive that for any
ε2 > 0, there exists Tε2 > Tε1 such that

u(x, t) ≤ max
{( r

1 + k(m− ε1)
− d
)
/a, 0

}
+ ε2



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 7

in Ω× [Tε2 ,∞). This and (1.3) show that

∂u

∂t
= du∆u+

ru

1 + kv
− du− au2 − b(x)p(u, v)v

≤ du∆u+
ru

1 + k(m− ε1)
− du− au2 − b(x)p(u,m− ε1)(m− ε1)

= du∆u+
ru

1 + k(m− ε1)
− du− au2 − b(x)pu(u′,m− ε1)(m− ε1)u

≤ du∆u+ Λu− au2

(2.2)

in Ω× [Tε2 ,∞), where u′ ∈ [0, u] and

Λ :=
r

1 + k(m− ε1)
− d

− b(x)pu

(
max

{( r

1 + k(m− ε1)
− d
)
/a, 0

}
+ ε2,m− ε1

)
(m− ε1).

Hence, the comparison principle shows that u(x, t) ≤ U(x, t) in Ω× [Tε2 ,∞), where
U(x, t) is the solution of

Ut = du∆U + ΛU − aU2, x ∈ Ω, t > Tε2 ,

∂U

∂µ
= 0, x ∈ ∂Ω, t > Tε2 ,

U(x, Tε2) = u(x, Tε2), x ∈ Ω.

We consider two possibilities:
(i) If r/(1 + km) ≤ d, then it is clear that

λ1

(
du,−

r

1 + km
+ d+ b(x)pu(0,m)m

)
> 0.

By choosing enough small ε1 > 0 and ε2 > 0, we find that λ1(du,−Λ) ≥ 0, and
hence U(x, t) → 0 uniformly in Ω as t → ∞. Thus, for any ε3 > 0, there exists
Tε3 > Tε2 such that u(x, t) ≤ ε3 in Ω × [Tε3 ,∞). Letting ε3 → 0, we obtain that
limt→∞ u(x, t) = 0 uniformly in Ω. This, together with the equation for v(x, t) of
(1.3), shows that

∂v

∂t
= dv∆v +mv − v2 + cb(x)p(u, v)v ≤ dv∆v + (m+ cb∗p(ε3,m− ε1)) v − v2

in Ω × [Tε3 ,∞). As above, the comparison principle shows that for any ε4 > 0,
there exists Tε4 > Tε3 such that v(x, t) ≤ m+ cb∗p(ε3,m− ε1) + ε4 in Ω× [Tε4 ,∞).
Letting εi → 0 with i = 1, 2, 3, 4, we have that limt→∞ v(x, t) = m uniformly in Ω.

(ii) If r/(1 + km) > d, then for any ε1 > 0 small, we have

Λ =
r

1 + k(m− ε1)
− d− b(x)pu

(( r

1 + k(m− ε1)
− d
)
/a+ ε2,m− ε1

)
(m− ε1).

Moreover, by choosing enough small ε1 > 0 and ε2 > 0 if necessary, we obtain
λ1 (du,−Λ) ≥ 0 by the hypothesis

λ1

(
du,−

r

(1 + km)
+ d+ b(x)pu

(( r

(1 + km)
− d
)
/a,m

)
m
)
> 0.

The rest of the argument is the same as that of the case (i) and hence is omitted. �



8 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

Remark 2.5. (i) Choosing p(u, v) = u/(1 + qu). By a straightforward modification
to the expression of (2.2), that is,

∂u

∂t
= du∆u+

ru

1 + kv
− du− au2 − b(x)uv

1 + qu

≤ du∆u+
ru

1 + k(m− ε1)
− du− au2 − b(x)(m− ε1)u

1 + qu

≤ du∆u+
ru

1 + k(m− ε1)
− du− au2

− b(x)(m− ε1)u

1 + q
(

max{( r
1+k(m−ε1) − d)/a, 0}+ ε2

)
≤ du∆u+

( r

1 + k(m− ε1)
− d

− b(x)(m− ε1)

1 + q
(

max{( r
1+k(m−ε1) − d)/a, 0}+ ε2

))u− au2,

the hypothesis

λ1

(
du,−

r

1 + km
+ d+ b(x)pu

(
(

r

1 + km
− d)/a,m

)
m
)
> 0

becomes

λ1

(
du,−

r

1 + km
+ d+

amb(x)

a+ q(r/(1 + km)− d)

)
> 0.

(ii) Choosing p(u, v) = u/(1 + qu+ fv). By a similar modification to the ex-
pression of (2.2) to that of (i), the hypothesis

λ1

(
du,−

r

1 + km
+ d+ b(x)pu

(( r

1 + km
− d
)
/a,m

)
m
)
> 0

can be rewritten as

λ1

(
du,−

r

1 + km
+ d+

amb(x)

a+ q(r/(1 + km)− d) + afm

)
> 0.

2.2. Existence of positive solutions to (1.4). In this subsection, we establish
the sufficient conditions for the existence of positive solutions to (1.4) by applying
degree theory in cones. The following lemma gives the L∞(Ω)-estimate for any
positive solution which is independent of the diffusion coefficients du and dv.

Lemma 2.6. Assume that r > d and (u, v) is any positive solution of (1.4). Then
0 ≤ u ≤ (r − d)/a and max{m, 0} ≤ v ≤ m+ cb∗p ((r − d)/a, 0) in Ω. Moreover, if
m > 0 and r/(1 + km) > d, then 0 ≤ u ≤ (r/(1 + km)− d)/a in Ω.

The proof of the above lemma is standard by a simple comparison argument,
hence is omitted. When the diffusion coefficients du and dv are away from 0, we
provice the W 2,p(Ω)-estimate with p ∈ (1,∞) for any positive solution of (1.4).

Lemma 2.7. Assume that r > d and (u, v) is any positive solution of (1.4) with
du ≥ ε and dv ≥ ε, where ε > 0 is any small number. Then there is a positive
number C = C(ε, r, d, a,m, c, b(x), |Ω|) such that ‖u‖W 2,p(Ω) ≤ C and ‖v‖W 2,p(Ω) ≤
C for any p ∈ (1,∞).



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 9

Proof. Suppose that (u, v) = (u(x), v(x)) is any positive solution of (1.4) with du ≥
ε and dv ≥ ε. Let Ci =: Ci(ε, n, r, d, a,m, c, b(x), |Ω|) for simplicity. Multiplying
(1.4) by u and integrate the resulting expression we obtain∫

Ω

|∇u|2dx =
1

du

∫
Ω

( ru2

1 + kv
− du2 − au3 − b(x)p(u, v)uv

)
dx

≤ 1

ε

∫
Ω

(r − d− au)u2dx.

Note that (r − d − au)u2 ≤ 4(r − d)3/27a2 for 0 ≤ u ≤ (r − d)/a. It follows

from Lemma 2.6 that
∫

Ω
|∇u|2dx ≤ 4(r−d)3

27εa2 |Ω|. Moreover, Lemma 2.6 shows that
ru

1+kv − du − au
2 − b(x)p(u, v)v is bounded for any p ∈ (1,∞). By the standard

Lp-estimate for elliptic equation, we have

‖u‖H2(Ω) ≤ C2

(
‖u‖H1(Ω) +

1

ε

∥∥ ru

1 + kv
− du− au2 − b(x)p(u, v)v

∥∥
L2(Ω)

)
≤ C3.

Therefore, we apply the Sobolev embedding theorem to have

u ∈ H2(Ω) ⊂


C1(Ω) if n = 1;

W 1,p1(Ω), ∀p1 ∈ [1,∞) if n = 2;

W 1,p1(Ω), ∀p1 ∈ [1, 2n
n−2 ] if n > 2.

Here n is the dimension of space. For the case that n > 2, the standard Lp-estimate
for elliptic equation yields

‖u‖W 2,p1 (Ω) ≤ C4

(
‖u‖W 1,p1 (Ω) +

1

ε

∥∥ ru

1 + kv
−du−au2− b(x)p(u, v)v

∥∥
Lp1 (Ω)

)
≤ C5.

Thus, applying the Sobolev embedding theorem we have

u ∈ ‖u‖W 2,p1 (Ω) ⊂


C1(Ω) if n < p1;

W 1,p2(Ω), ∀p2 ∈ [1,∞) if n = p1;

W 1,p2(Ω), ∀p2 ∈ [1, p1n
n−p1 ] if n > p1.

Consequently, for any p ∈ (1,∞), we can repeat the above argument if necessary
to find a positive constant C6 such that ‖u‖W 2,p(Ω) ≤ C6.

Similarly, we can obtain the desired estimate for ‖v‖W 2,p(Ω). �

Assume that E is a real Banach space. For a closed convex set W in E, it is said
to be a wedge if αW ⊂W for all α ≥ 0. Moreover, when W∩{−W} = 0, W is called
a cone. Define Wy =: cl{x ∈ E : y+ωx ∈W}, where y ∈W , ω > 0 is some constant
and “cl” represents the closure of the set. Clearly, Wy is a wedge, and its maximal
linear subspace is denoted by Sy. Let A be a Fréchet differentiable compact operator
in E such that y ∈ W is a fixed point of A and A(W ) ⊆ W . Then Wy and Sy are
both invariant under L(y), where L(y) is the Fréchet derivative of A at y. Let U be
an open subset of W . Define indexW (A, U) = index(A, U,W ) = degW (I −A, U, 0),
where I is the identity map. Assume that y is an isolated fixed point of A. Then
the index of A at y in W is given by indexW (A, y) = index(A, U(y),W ), where
U(y) is a small open neighborhood of y in W .

Assume that E has the decomposition E = Ey ⊕ Sy, where Ey is a closed linear
subspace of E, and Wy is generating. Then the index of A at y can be calculated
by the following proposition (see [20, Theorems 2.2 and 2.3]).



10 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

Proposition 2.8. Suppose that Q : E → Ey is a projection operator of Ey along
Sy. If L(y) has no non-zero fixed point on Wy, then index(A, y) exists. Moreover,
the following assertions hold.

(1) indexW (A, y) = 0 if Q ◦ L(y) has an eigenvalue greater than 1 on Wy;
(2) indexW (A, y) = indexSy

(L(y), 0) = (−1)σ if Q ◦ L(y) has no eigenvalue
greater than 1 on Wy, where σ is the sum of algebraic multiplicities of the
eigenvalues of L(y) restricted in Sy which are greater than 1.

Let E = C(Ω)⊕C(Ω) and W = P ⊕P , where P = {φ ∈ C(Ω) : φ(x) ≥ 0 in Ω}.
Set D =

{
(u, v) ∈ W : u ≤ (r − d)/a + 1, v ≤ m + 1 + cb∗p ((r − d)/a, 0)

}
. For

each t ∈ [0, 1], we define a Fréchet differentiable compact operator At : D → E by

At
(
u
v

)
= (−∆ +M)−1

( 1
du

(
ru

1+kv − du− au
2 − tb(x)p(u, v)v

)
+Mu

1
dv

(
mv − v2 + tcb(x)p(u, v)v

)
+Mv

)
.

Here M is a suitably large number such that

(du)−1
( ru

1 + kv
− du− au2 − tb(x)p(u, v)v

)
+Mu > 0

and (dv)
−1
(
mv − v2 + tcb(x)p(u, v)v

)
+Mv > 0 for any (u, v) ∈ D.

Assume that r > d and m > 0. Then (0, 0), ((r − d)/a, 0) and (0,m) are the only
non-negative fixed points of A1 which are not positive. We can apply Proposition
2.8 to calculate their indices.

Lemma 2.9. Assume that r > d and m > 0. Then

(1) indexW (A1, (0, 0)) = 0;
(2) indexW (A1, ((r − d)/a, 0)) = 0;

(3) indexW (A1, (0,m)) =

{
0, if λ1

(
du,− r

1+km + d+ b(x)pu(0,m)m
)
< 0,

1, if λ1

(
du,− r

1+km + d+ b(x)pu(0,m)m
)
> 0;

(4) indexW (A1, D) = 1 if r/(1 + km) 6= d.

Proof. Since the proofs of (1)–(4) are similar and the proof of (3) is a little more
complicated, we only prove (3). By the definitions of W(0,m) and S(0,m), it is easy

to check that W(0,m) = P ⊕ C(Ω) and S(0,m) = {0} ⊕ C(Ω). Choosing E(0,m) =

C(Ω)⊕ {0}. Then E = S(0,m) ⊕ E(0,m). Define

L1(0,m) = (−∆ +M)−1

( 1
du

(
r

1+km − d− b(x)pu(0,m)m
)

+M 0
1
dv

(cb(x)pu(0,m)m) − 1
dv
m+M

)
.

Then L1(0,m) is the Fréchet derivative of A1 with respect to (u, v) at (0,m). As-
sume that (φ, ψ) ∈W(0,m) is an eigenfunction of L1(0,m) and κ is the corresponding
eigenvalue. Then

(−∆ +M)−1
{ 1

du

( r

1 + km
− d− b(x)pu(0,m)m

)
+M

}
φ = κφ, x ∈ Ω,

(−∆ +M)−1
{ 1

dv

(
cb(x)pu(0,m)m

)}
φ

+ (−∆ +M)−1
{
− 1

dv
m+M

}
ψ = κψ, x ∈ Ω,

∂φ

∂µ
=
∂ψ

∂µ
= 0, x ∈ ∂Ω.



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 11

When λ1 (du,−r/(1 + km) + d+ b(x)pu(0,m)m) 6= 0 and m > 0, it is easy to check
that I − L1(0,m) is invertible on W(0,m) by contradiction. Hence, L1(0,m) has no
non-zero fixed point on W(0,m). Consequently, it follows from Proposition 2.8 that
indexW (A1, (0,m)) exists.

To obtain the value of indexW (A1, (0,m)), we now analyze the eigenvalues of
Q ◦ L1(0,m), where Q : E = C(Ω) ⊕ C(Ω) → E(0,m) = C(Ω) ⊕ {0} is a projec-
tion operator of E(0,m) along S(0,m). By the definition of Q, we see that every
eigenfunction of Q ◦L1(0,m) has the form (φ, 0), where φ is a non-zero solution of

(−∆ +M)−1
{ 1

du

( r

1 + km
− d− b(x)pu(0,m)m

)
+M

}
φ = κφ, x ∈ Ω,

∂φ

∂µ
= 0, x ∈ ∂Ω.

By [14, Lemma 2.4], it is well known that

r
[
(−∆ +M)−1

{ 1

du

( r

1 + km
− d− b(x)pu(0,m)m

)
+M

}]
> 1 (resp. < 1)

if λ1

(
du,− r

1+km + d + b(x)pu(0,m)m
)
< 0 (resp.> 0). Therefore, Q ◦ L1(0,m)

has an eigenvalue larger than 1 when λ1

(
du,− r

1+km + d + b(x)pu(0,m)m
)
< 0,

and thus indexW (A1, (0,m)) = 0 by Proposition 2.8(1). On the other hand, if
λ1

(
du,− r

1+km +d+b(x)pu(0,m)m
)
> 0, then Q◦L1(0,m) has no eigenvalue larger

than or equal to 1. Hence, by Proposition 2.8(2), we see that indexW (A1, (0,m)) =
(−1)σ, where σ is the sum of algebraic multiplicities of the eigenvalues of L1(0,m)
restricted in S(0,m) which are greater than 1. Indeed, we can prove σ = 0. Assume
that (φ, ψ) ∈ S(0,m) is an eigenfunction of L1(0,m) and κ is the corresponding
eigenvalue. Then φ = 0 and ψ is a non-zero solution of

(−∆ +M)−1
{
− m

dv
+M

}
ψ = κψ, ∈ Ω,

∂ψ

∂µ
= 0, x ∈ ∂Ω.

Since λ(dv,m) = m > 0, we derive from [14, Lemma 2.4] that r[(−∆+M)−1(−m
dv

+

M)] < 1. This implies that L1(0,m) has no eigenvalue larger than or equal to 1 in
S(0,m), and thus σ = 0. Consequently, indexW (A1, (0,m)) = 1. �

Assume that r > d and m < 0. Then (0, 0) and ((r − d)/a, 0) are the only non-
negative fixed points of A1 which are not positive. There is a result corresponding
to Lemma 2.9 for the case m < 0.

Lemma 2.10. Assume that r > d and m < 0. Then

(1) indexW (A1, (0, 0)) = 0;

(2) indexW (A1, ((r − d)/a, 0)) =

{
0, if λ1(dv,−m− cb(x)p( r−da , 0)) < 0,

1, if λ1(dv,−m− cb(x)p( r−da , 0)) > 0;

(3) indexW (A1, D) = 1.

With the help of Lemmas 2.9 and 2.10, we use the excision property of fixed
point index to obtain the existence of positive solutions to (1.4) as follows.

Theorem 2.11. The following assertions hold.

(1) Assume that r/(1 + km) > d and m > 0. Then system (1.4) admits at least
one positive solution if λ1 (du,−r/(1 + km) + d+ b(x)pu(0,m)m) < 0.

(2) Assume that r > d and m < 0. Then system (1.4) admits at least one
positive solution if and only if λ1 (dv,−m− cb(x)p ((r − d)/a, 0)) < 0.



12 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

Proof. We argue indirectly. For the case that r/(1 + km) > d and m > 0. Assume
that (1.4) has no positive solution when

λ1

(
du,−

r

1 + km
+ d+ b(x)pu(0,m)m

)
< 0.

Then from Lemma 2.9, the additivity property of fixed point indices yields

1 = indexW (A1, D)

= indexW (A1, (0, 0)) + indexW (A1, ((r − d)/a, 0)) + indexW (A1, (0,m)) = 0.

This contradiction shows that when

λ1

(
du,−

r

1 + km
+ d+ b(x)pu(0,m)m

)
< 0,

(1.4) admits at least one positive solution. Similarly, we can show that (1.4) admits
at least one positive solution when λ1 (dv,−m− cb(x)p ((r − d)/a, 0)) < 0.

It remains to prove that (1.4) has no positive solution if

λ1

(
dv,−m− cb(x)p

(r − d
a

, 0
))
≥ 0.

On the contrary, assume that (1.4) has a positive solution (u, v). Then from the
equation for u and Lemma 2.6, we deduce that

0 = λ1 (dv,−m+ v − cb(x)p(u, v))

> λ1(dv,−m− cb(x)p(u, 0))

≥ λ1 (dv,−m− cb(x)p((r − d)/a, 0)) ≥ 0.

This is impossible which proves the desired result. �

3. Effect of fear

The purpose of this section is to obtain some insights on how fear affects the
dynamics of (1.3). To make the analysis more explicit, we consider three particular
forms for the functional response p(u, v).

3.1. Linear functional response. In this subsection, we choose the linear func-
tional response (i.e., p(u, v) = u) to demonstrate the effect of fear on the dynamics
of (1.3).

The zero level curve of λ1 (du,−r/(1 + km) + d+mb(x)) = 0 is given by r =
r(k;m) =: (1+km)λ1 (du, d+mb(x)). By Proposition 2.2(1), we have the following
lemma.

Lemma 3.1. For any k ≥ 0, the function r(k;m) is continuously differentiable and
monotone increasing with respect to m > 0, moreover it satisfies limm→0+ r(k;m) =
d and limm→∞ r(k;m) = ∞. Furthermore, λ1 (du,−r/(1 + km) + d+mb(x)) < 0
if r > r(k;m); > 0 if r < r(k;m).

Denote m = m(r) =: λ1 (dv,−c(r − d)b(x)/a). Then we have the following
properties of m(r).

Lemma 3.2. For r > d, the function m(r) is continuously differentiable and
monotone decreasing with respect to r, moreover it satisfies limr→d+ m(r) = 0 and
limr→∞m(r) = −∞. Furthermore, λ1 (dv,−m− c(r − d)b(x)/a) < 0 if m > m(r);
> 0 if m < m(r).



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 13

Suppose that (1.4) with p(u, v) = u has a positive solution (u, v). Then by
Lemma 2.6, we have

0 = λ1

(
du,−

r

1 + kv
+ d+ au+ b(x)v

)
> λ1

(
du,−

r

1 + km
+ d+mb(x)

)
for any m > 0. Therefore, combined with Theorems 2.4 and 2.11, we summarize
the above discussion to obtain the following result which can be described in the
rm-plane (see Figure 1).

Figure 1. Dynamics of (1.3) with linear functional response.

Theorem 3.3. Let p(u, v) = u. Then the following assertions hold.
(1) Assume that m < 0. Then any non-negative solution of (1.3) converges to

(0, 0) uniformly in Ω as t→∞ for r ≤ d, any non-negative solution of (1.3) with
u0(x) ≥ (6≡)0 converges to ((r − d)/a, 0) uniformly in Ω as t → ∞ for d < r <
m−1(r), where m−1(r) is the inverse function of m(r) =: λ1 (dv,−cb(x)(r − d)/a),
and (1.3) admits at least one positive steady state for r > m−1(r).

(2) Assume that m > 0. Then any non-negative solution of (1.3) with v0(x) ≥
(6≡)0 converges to (0,m) uniformly in Ω as t → ∞ for r < r(k;m) =: (1 +
km)λ1 (du, d+ b(x)m), and (1.3) admits at least one positive steady state for r >
r(k;m).

3.2. Holling-type II functional response. In this subsection, we choose the
Holling-type II functional response (i.e., p(u, v) = u/(1 + qu)) to demonstrate the
effect of fear on the dynamics of (1.3).

It is clear that the zero level curve of λ1 (dv,−m− c(r − d)b(x)/(a+ q(r − d))) =
0 is given by m = m(r) =: λ1 (dv,−c(r − d)b(x)/(a+ q(r − d))). By Proposition
2.2(1), one can easily obtain the following properties of m(r).

Lemma 3.4. For r > d, the function m(r) is continuously differentiable and
monotone decreasing with respect to r, moreover it satisfies limr→d+ m(r) = 0 and

limr→∞m(r) = λ1 (dv,−cb(x)/q). Furthermore, λ1

(
dv,−m− (r−d)cb(x)

a+q(r−d)

)
< 0 if

m > m(r); > 0 if m < m(r).



14 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

From the equation for u and Lemma 2.6, a necessary condition for the existence
of positive solutions to (1.4) with p(u, v) = u/(1 + qu) is

Ψ(r,m) =: λ1

(
du,−

r

1 + km
+ d+

amb(x)

a+ q(r/(1 + km)− d)

)
< 0.

The following lemma gives the zero level curve of Ψ(r,m).

Lemma 3.5. For any k ≥ 0 and r/(1 + km) > d, there exists a monotone in-
creasing function r(k;m) with respect to m such that Ψ(r(k;m),m) = 0 for any
m > 0, moreover it satisfies d(1 + km) < r(k;m) < r(k;m) for any m > 0,
limm→0+ r(k;m) = d and limm→∞ r(k;m) = ∞. Furthermore, Ψ(r,m) < 0 if
r > r(k;m); > 0 if r < r(k;m).

Proof. Define

ψ(r,m) =: − r

1 + km
+ d+

amb(x)

a+ q(r/(1 + km)− d)
.

Then Ψ(r,m) = λ1 (du, ψ(r,m)). Since ψ(r,m) is monotone decreasing with respect
to r and is monotone increasing with respect to m, it follows from Proposition 2.2(1)
that Ψ(r,m) is monotone decreasing with respect to r and is monotone increasing
with respect to m. Since ψ(r,m) → −∞ as r → ∞ and ψ(r,m) → mb(x) > 0 as
r → (d(1 + km))+, we deduce from the variational characterization that

lim
r→∞

Ψ(r,m) = −∞

and

lim
r→(d(1+km))+

Ψ(r,m) =

∫
Ω

(
du|∇φ|2 +mb(x)φ2

)
dx
/∫

Ω

φ2dx > 0

for some φ ∈ H1(Ω), where φ does not change sign in Ω. By the intermediate value
theorem, there exists a unique number r(k;m) such that Ψ(r(k;m),m) = 0 for any
m > 0. Moreover, it follows from implicit function theorem that r(k;m) ∈ C1(0,∞)
is monotone increasing with respect to m.

Since

λ1

(
du,−

r

1 + km
+ d+

amb(x)

a+ q(r/(1 + km)− d)

)
< λ1

(
du,−

r

1 + km
+ d+mb(x)

)
we have λ1 (du,−r(k;m)/(1 + km) + d+mb(x)) > 0. Thus, it follows from Lemma
3.1 that r(k;m) < r(k;m) for any m > 0. On the other hand,

λ1

(
du,−

r

1 + km
+ d+

amb(x)

a+ q(r/(1 + km)− d)

)
> λ1

(
du,−

r

1 + km
+ d
)

= − r

1 + km
+ d,

we deduce from this inequality to obtain r(k;m) > d(1+km) for any m > 0. Hence,
we proved that d(1 + km) < r(k;m) < r(k;m) for any m > 0. This, together with
Lemma 3.1, implies that limm→0+ r(k;m) = d and limm→∞ r(k;m) = ∞. This
completes the proof. �

Therefore, combining Theorems 2.4 with 2.11, we obtain the following result
which can be described in the rm-plane (see Figure 2).



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 15

Figure 2. Dynamics of (1.3) with Holling-type II functional response.

Theorem 3.6. Let p(u, v) = u/(1 + qu). The following assertions hold.
(1) Assume that m ≤ λ1 (dv,−cb(x)/q). Then any non-negative solution of

(1.3) converges to (0, 0) uniformly in Ω as t→∞ for r ≤ d, and any non-negative
solution of (1.3) with u0(x) ≥ (6≡)0 converges to ((r − d)/a, 0) uniformly in Ω as
t→∞ for r > d.

(2) Assume that λ1 (dv,−cb(x)/q) < m < 0. Then any non-negative solution
of (1.3) converges to (0, 0) uniformly in Ω as t → ∞ for r ≤ d, any non-negative
solution of (1.3) with u0(x) ≥ (6≡)0 converges to ((r − d)/a, 0) uniformly in Ω as
t → ∞ for d < r < m−1(r), where m−1(r) is the inverse function of m(r) =:
λ1 (dv,−c(r − d)b(x)/(a+ q(r − d))), and (1.3) admits at least one positive steady
state for r > m−1(r).

(3) Assume that m > 0. Then any non-negative solution of (1.3) with v0(x) ≥
( 6≡)0 converges to (0,m) uniformly in Ω as t→∞ for r < r(k;m) given in Lemma
3.5, and (1.3) admits at least one positive steady state for r > r(k;m) given in
Lemma 3.1.

For any fixed m > 0, Theorem 3.6 provides no information on the existence
and non-existence of positive steady states to (1.3) with p(u, v) = u/(1 + qu) for
r ∈ (r(k;m), r(k;m)). A further study is therefore necessary in order to better
understand the dynamics of (1.3) with p(u, v) = u/(1 + qu).

We now apply the bifurcation result of Crandall-Rabinowitz in [4] to obtain a
branch of positive solutions to (1.4) with p(u, v) = u/(1 + qu) emanating from
(r, u, v) = (r(k;m), 0,m) with m > 0. Recall that X = {u ∈W 2,p

µ (Ω) : ∂u/∂µ = 0}
and Y = Lp(Ω). Let B : R×X ×X → Y × Y be given by

B(r, u, v) =

(
du∆u+ ru

1+kv − du− au
2 − b(x)uv

1+qu

dv∆v +mv − v2 + cb(x)uv
1+qu

)
.

It follows from the Krein-Rutman theorem that B(u,v)(r, 0,m)(φ, ψ) = 0 has a

solution with φ > 0 in Ω if and only if r = r(k;m). Moreover, the kernel



16 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

N (B(u,v)(r(k;m), 0,m)) = span{(φ1, ψ1)}, where φ1 is a positive solution of

−du∆φ+
(
− r(k;m)

1 + km
+ d+ b(x)m

)
φ = 0, x ∈ Ω,

∂φ

∂µ
= 0, x ∈ ∂Ω,

and ψ1 = (−dv∆ +mI)−1(cmb(x)φ1). Moreover, it follows from Fredholm alterna-
tive theorem that the range R(B(u,v)(r(k;m), 0,m)) = {(φ, ψ) ∈ Y 2 :

∫
Ω
φφ1dx =

0}. In addition, a simple calculation yields

Br(u,v)(r(k;m), 0,m)(φ1, ψ1) = (φ1/(1 + km), 0) 6∈ R(B(u,v)(r(k;m), 0,m)).

Therefore, by applying the local bifurcation theorem in [4], we conclude that posi-
tive solutions of (1.4) with p(u, v) = u/(1 + qu) near (r(k;m), 0,m) lie in a smooth
curve Σε = {(r(s), u(s), v(s)) : s ∈ (0, ε)}, where (r(s), u(s), v(s)) is continuously
differentiable with respect to s and (r(0), u(0), v(0)) = (r(k;m), sφ1 + o(s),m +
sψ1 + o(s)). Furthermore, by the formula (4.5) in [22], we see that

r′(0) =
(∫

Ω

φ2
1

1 + km
dx
)−1[ ∫

Ω

(a− qmb(x))φ3
1dx

+

∫
Ω

(kλ1(du, d+ b(x)m)

1 + km
+ b(x)

)
φ2

1ψ1dx
]
.

Consequently, the bifurcation of Σε at (r(k;m), 0,m) is subcritical (r′(0) < 0) if
q > q0(k) and it is supercritical (r′(0) > 0) if 0 ≤ q < q0(k), where

q0(k) =
(∫

Ω

mb(x)φ3
1dx
)−1

∫
Ω

(
aφ1 +

(
b(x) +

kλ1(du, d+ b(x)m)

1 + km

)
ψ1

)
φ2

1dx.

Remark 3.7. When q ∈ (q0(0), q0(k)), the level of fear k alters the direction of
bifurcation from subcritical (r′(0) < 0) to supercritical (r′(0) > 0). This, together
with Theorem 3.6(3), indicates that the cost of fear will not only affect the exis-
tence of positive solutions to (1.4), but also change the direction of steady-state
bifurcation.

A global bifurcation consideration, together with the maximum principle, proves
that the local curve Σε is contained in a global branch of positive solutions to (1.4)
with p(u, v) = u/(1 + qu) which is denoted by Σ = {(r, u, v)}. Furthermore, by the
similar argument to that of Theorem 2.4 in [15], we can show that the global branch
Σ is unbounded in R×X×X. It follows from Lemma 2.6 that any positive solution
(u, v) of (1.4) with p(u, v) = u/(1 + qu) satisfies 0 ≤ u ≤ 1

a (r/(1 + km)− d) and
m ≤ v ≤ m + cb∗(r − d)/(a+ q(r − d)) for any m > 0. Thus, the only possible is
(r(k;m),∞) ⊂ {r : (r, u, v) ∈ Σ}. Summarizing the above discussion, we obtain
the following result.

Proposition 3.8. Assume that r/(1 + km) > d and m > 0. Then there ex-
ists a global unbounded continuum Σ of positive solutions to (1.4) with p(u, v) =
u/(1 + qu) which contains Σε such that Σ emanates from (r, u, v) = (r(k;m), 0,m)
and tends to infinity as r goes to infinity. Moreover, the bifurcation of Σε at
(r(k;m), 0,m) is subcritical (r′(0) < 0) if q > q0(k) and it is supercritical (r′(0) > 0)
if 0 ≤ q < q0(k).

Based on Lemma 2.9 and Proposition 3.8, we establish the multiplicity of positive
solutions to (1.4) with p(u, v) = u/(1 + qu) by the similar argument to that of
Theorem 2.3 in [15].



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 17

Theorem 3.9. Assume that r/(1 + km) > d and m > 0. If q > q0(k), then there
exists a constant r ∈ (r(k;m), r(k;m)) such that (1.4) with p(u, v) = u/(1 + qu)
has at least two positive solutions for r ∈ (r, r(k;m)) and has at least one positive
solution for r ∈ [r(k;m),∞).

3.3. Beddington-DeAngelis functional response. In this subsection, we choose
the Beddington-DeAngelis functional response (i.e., p(u, v) = u/(1 + qu+ fv)) to
demonstrate the effect of fear on the dynamics of (1.3).

The zero level curve of λ1 (du,−r/(1 + km) + d+mb(x)/(1 + fm)) = 0 is given
by r = r̃(k;m) =: (1 + km)λ1 (du, d+mb(x)/(1 + fm)). We have the following
lemma.

Lemma 3.10. For any k ≥ 0, the function r̃(k;m) is continuously differentiable
and monotone increasing with respect to m > 0, moreover it satisfies d(1 + km) <
r̃(k;m) < r(k;m) for any m > 0, and

lim
m→0+

r̃(0;m) = d, lim
m→∞

r̃(0;m) = λ1 (du, d+ b(x)/f) ,

lim
m→0+

r̃(k;m) = d, lim
m→∞

r̃(k;m) =∞ with k > 0.

Furthermore, λ1 (du,−r/(1 + km) + d+mb(x)/(1 + fm)) < 0 if r > r̃(k;m); > 0
if r < r̃(k;m).

From the equation for u and Lemma 2.6, a necessary condition for the exis-
tence of positive solutions to (1.4) with p(u, v) = u/(1 + qu+ fv) is Φ(r,m) =:

λ1

(
du,− r

1+km + d + amb(x)
a+q(r/(1+km)−d)+afm

)
< 0. By the similar argument to that

of Lemma 3.5, we have the following lemma.

Lemma 3.11. For any k ≥ 0 and r/(1 + km) > d, there exists a monotone in-
creasing function r̂(k;m) with respect to m such that Φ(r̂(k;m),m) = 0 for any
m > 0, moreover it satisfies d(1 + km) < r̂(k;m) < r̃(k;m) for any m > 0, and

lim
m→0+

r̂(0;m) = d, lim
m→∞

r̂(0;m) = λ1 (du, d+ b(x)/f) ,

lim
m→0+

r̂(k;m) = d, lim
m→∞

r̂(k;m) =∞ with k > 0.

Furthermore, Φ(r,m) < 0 if r > r̂(k;m); and Φ(r,m) > 0 if r < r̂(k;m).

Combined this with Theorems 2.4 and 2.11, we obtain the following result which
can be described in the rm-plane (see Figure 3).

Theorem 3.12. Let p(u, v) = u/(1 + qu+ fv).
(1) Assume that m ≤ λ1 (dv,−cb(x)/q). Then any non-negative solution of

(1.3) converges to (0, 0) uniformly in Ω as t→∞ for r ≤ d, and any non-negative
solution of (1.3) with u0(x) ≥ (6≡)0 converges to ((r − d)/a, 0) uniformly in Ω as
t→∞ for r > d.

(2) Assume that λ1 (dv,−cb(x)/q) < m < 0. Then any non-negative solution
of (1.3) converges to (0, 0) uniformly in Ω as t → ∞ for r ≤ d, any non-negative
solution of (1.3) with u0(x) ≥ (6≡)0 converges to ((r − d)/a, 0) uniformly in Ω as
t → ∞ for d < r < m−1(r), where m−1(r) is the inverse function of m(r) =:
λ1 (dv,−c(r − d)b(x)/(a+ q(r − d))), and (1.3) admits at least one positive steady
state for r > m−1(r).

(3) Assume that m > 0. Then any non-negative solution of (1.3) with v0(x) ≥
( 6≡)0 converges to (0,m) uniformly in Ω as t→∞ for r < r̂(k;m) given in Lemma



18 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

3.11, and (1.3) admits at least one positive steady state for r > r̃(k;m) given in
Lemma 3.10.

(a) k = 0 (b) k > 0

Figure 3. Dynamics of (1.3) with Beddington-DeAngelis func-
tional response.

Our next goal is to investigate the uniqueness and stability of positive solutions
to (1.4) with p(u, v) = u/(1 + qu+ fv) when the competition of the predator is
strong (i.e., f is large). Since the analyses are similar to those of subsection 2.3.2
in [15], we leave the details of the proof to the interested reader and only state the
main results. For the case m < 0, the following result holds.

Theorem 3.13. Assume that r > d and m ∈
(
λ1

(
dv,− c(r−d)b(x)

a+q(r−d)

)
, 0
)
. There exists

F large such that if f > F , then (1.4) with p(u, v) = u/(1 + qu+ fv) has a unique
and asymptotically stable positive solution (uf , vf ), moreover (uf , vf ) is close to
((r − d)/a, 0) and fvf is close to V , where V is the unique positive solution of

− dv∆V =
(
m+

c(r − d)b(x)

a+ q(r − d) + aV

)
V, x ∈ Ω,

∂V

∂µ
= 0, x ∈ ∂Ω. (3.1)

We want to point out that the cost of fear has no qualitative impact on the
behavior of positive solutions (1.4) for the case m < 0. However, for the case m > 0,
the cost of fear decreases the densities of the prey population. More precisely, we
have the following result.

Theorem 3.14. Assume that m > 0 and r/(1+km) > d. There exists F large such
that if f > F , then (1.4) with p(u, v) = u/(1 + qu+ fv) has no positive solution for
r < r̃(k;m) and has a unique and asymptotically stable positive solution (uf , vf )
for r > r̃(k;m) which is close to ((r/(1 + km)− d)/a,m).

To conclude this section, we state the effects of fear on the population dynamics
of (1.3). In view of Theorem 3.3, 3.6 and 3.12, one finds that a common feature
is that the cost of fear has no qualitative impact on the dynamics of (1.3) for
m < 0, while the cost of fear has obvious impact on the dynamics of (1.3) for
m > 0. More precisely, when the predator individuals are generalists, the cost
of fear makes the predator-only steady state (0,m) more stable by excluding the



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 19

existence of positive steady states under certain conditions. A special numerical
simulation example is presented in Figure 4 to verify this result. In particular,
the cost of fear seems to have a more profound effect on the dynamics of (1.3)
by choosing p(u, v) = u/(1 + qu+ fv). To be more specific, when the cost of
fear is ignored, the predator-only steady state (0,m) is never a stable one for any
r > λ1 (du, d+ b(x)/f) and all m > 0, and a positive steady state always exists, this
implies that two species can coexist; while when the cost of fear is considered, the
predator-only steady state (0,m) is globally asymptotically stable for any fixed r ∈
R and large m, and hence no positive steady state exists, this shows that the prey
will become extinct unconditionally and the predator will persist unconditionally
for any initial population distribution. Additionally, Theorem 3.14 shows that when
the competition of the predator is strong (i.e., f is large), the cost of fear makes the
total population of the prey drop form (r−d)/a · |Ω| to ((r/(1+km)−d)/a,m) · |Ω|,
where |Ω| represents the area of the region Ω. This implies that the cost of fear
not only makes it more difficult for the prey to survive, but also reduce the total
population of the prey even though the prey population survives.

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

u 
or

 v

u
v

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

u 
or

 v

u
v

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0

0.02

0.04

0.06

0.08

0.1

0.12

u 
or

 v

u
v

0 0.5 1 1.5 2 2.5 3 3.5 4
x

-2

0

2

4

6

8

10

12

u 
or

 v

10-3

u
v

0 0.5 1 1.5 2 2.5 3 3.5 4
x

-2

0

2

4

6

8

10

12

u 
or

 v

10-3

u
v

0 0.5 1 1.5 2 2.5 3 3.5 4
x

-2

0

2

4

6

8

10

12

u 
or

 v

10-3

u
v

(a) q = 0 and f = 0 (b) q = 0.1 and f = 0 (c) q = 0.1 and f = 1

Figure 4. Fix Ω = (0, 4), b(x) = 1 + x2, u0(x) = 0.1 + 0.1 sin(5x)
and v0(x) = 0.01 + 0.01 sin(5x). By choosing du = 11, dv = 1,
r = 6, d = 5.17, a = 10, m = 0.01, c = 0.6, the effect of the cost
of fear is revealed for k = 0 in the first line and for k = 16 in the
second line.

4. Effect of spatial diffusion and environmental heterogeneity

The purpose of this section is to investigate the impact of spatial diffusion and
environmental heterogeneity on the dynamics of (1.3). To make the analysis more
explicit, we choose the linear functional response (i.e., p(u, v) = u) in the equations
(1.3) and (1.4).

By Proposition 2.2, we can easily derive the following properties of

λ1

(
du,−

r

1 + km
+ d+ b(x)m

)



20 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

with respect to the diffusion rate of the prey individuals.

Lemma 4.1. Assume that r/(1 + km) > d and m > 0. Then

(1) λ1(du,−r/(1 + km)+d+b(x)m) is strictly monotone increasing with respect
to du > 0.

(2) If d < r/(1 + km) ≤ d + b∗m, then λ1(du,−r/(1 + km) + d + b(x)m) > 0
for all du > 0.

(3) If d + b∗m < r/(1 + km) < d + b(x)m, then there is a unique value Du =
Du(r,m, d, k, b(x)) ∈ (0,∞) such that λ1(du,−r/(1 + km)+d+b(x)m) < 0
for each du ∈ (0, Du) and λ1(du,−r/(1 + km) + d + b(x)m) > 0 for each
du ∈ (Du,∞).

(4) If r/(1 + km) ≥ d+ b(x)m, then λ1(du,−r/(1 + km) + d+ b(x)m) < 0 for
all du > 0.

Based on Theorem 2.4(3), Theorem 2.11(1) and Lemma 4.1, the dynamic be-
havior of (1.3) with p(u, v) = u with respect to the diffusion rate du is as follows.

Theorem 4.2. Let p(u, v) = u, and assume that r/(1 + km) > d and m > 0. Then

(1) If d < r/(1 + km) ≤ d + b∗m, then (0,m) is globally asymptotically stable
for all du > 0.

(2) If d+b∗m < r/(1 + km) < d+b(x)m, then (0,m) is globally asymptotically
stable for each du ∈ (Du,∞); (0,m) is unstable and (1.3) admits at least
one positive steady state for each du ∈ (0, Du).

(3) If r/(1 + km) ≥ d+b(x)m, then (0,m) is unstable and (1.3) admits at least
one positive steady state for all du > 0.

Next we investigate the asymptotic profiles of positive steady states to (1.3) with
p(u, v) = u as as the diffusion rate of the prey or predator individuals approaches
zero or infinity. The following theorem gives the asymptotic behavior of any positive
solution to (1.4) with p(u, v) = u as du goes to∞, and a special numerical simulation
example is presented in Figure 5.

Theorem 4.3. Assume that m > 0 and r/(1 + km) > d + mb(x). Then for fixed
dv > 0, any positive solution (u, v) of (1.4) with p(u, v) = u satisfies

(u, v)→ (u∞, v∞) in C1(Ω)× C1(Ω) as du →∞,
where

u∞ =
1

a|Ω|

∫
Ω

( r

1 + kv∞
− d− b(x)v∞

)
dx

and v∞ is a positive solution of

−dv∆v∞ =
(
m+

cb(x)

a|Ω|

∫
Ω

( r

1 + kv∞
− d− b(x)v∞

)
dx
)
v∞ − v2

∞, x ∈ Ω,

∂v∞
∂µ

= 0, x ∈ ∂Ω.

(4.1)

Proof. From Theorem 4.2(3), our assumption conditions ensure that (1.4) with
p(u, v) = u admits at least one positive solution for all du > 0. Hence, we may
assume that (udu,n

, vdu,n
) is a positive solution of (1.4) with p(udu,n

, vdu,n
) = udu,n

and du = du,n, where du,n → ∞ as n → ∞. By Lemma 2.7 and the embedding
theorem, we can assume that

(udu,n
, vdu,n

)→ (u∞, v∞) in C1(Ω)× C1(Ω) as n→∞.



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 21

As n → ∞, it follows from the equation for udu,n
and Lemma 2.6 that u∞ is a

non-negative constant. If we assume that u∞ = 0 in Ω, then from the equation for
vdu,n

, we find that v∞ satisfies

−dv∆v∞ = mv∞ − v2
∞, x ∈ Ω,

∂v∞
∂µ

= 0, x ∈ ∂Ω.

By Lemma 2.6, vdu,n
≥ m > 0 for all n ∈ N+, and thus v∞ ≥ m > 0 in Ω. Hence,

we have v∞ = m in Ω. Let ũdu,n = udu,n/‖udu,n‖L∞(Ω). Then from the equation
for udu,n it follows that

−du,n∆ũdu,n
=
( r

1 + kvdu,n

− d− audu,n
− b(x)vdu,n

)
ũdu,n

, x ∈ Ω,

∂ũdu,n

∂µ
= 0, x ∈ ∂Ω.

By the standard elliptic regularity theory, we can assume that ũdu,n
→ 1 in C1(Ω)

since du,n →∞ as n→∞. Moreover, for all n ∈ N+, we integrate the equation for
ũdu,n

over Ω to obtain∫
Ω

( r

1 + kvdu,n

− d− audu,n
− b(x)vdu,n

)
ũdu,n

dx = 0.

Letting n → ∞, we obtain r/(1 + km) = d + mb(x). This contradiction implies
that u∞ is a positive constant. Hence, we deduce from the equation for udu,n that

u∞ =
1

a|Ω|

∫
Ω

( r

1 + kv∞
− d− b(x)v∞

)
dx.

This implies that v∞ is a non-negative solution of (4.1).
It remains to prove that v∞ is a positive solution of (4.1). Otherwise, we apply

Lemma 2.6 and the Harnack inequality to (4.1) to derive that v∞ ≡ 0 in Ω. This
shows that u∞ = (r − d)/a. Let ṽdu,n

= vdu,n
/‖vdu,n

‖L∞(Ω). Then it follows from
the equation for vdu,n

that

−dv∆ṽdu,n
=
(
m− vdu,n

+ cb(x)udu,n

)
ṽdu,n

, x ∈ Ω,
∂ṽdu,n

∂µ
= 0, x ∈ ∂Ω.

By the standard elliptic regularity theory, we can assume that ṽdu,n
→ ṽ∞ in C1(Ω)

as n→∞. Moreover, ṽ∞ satisfies

−dv∆ṽ∞ = (m+ cb(x)(r − d)/a) ṽ∞, x ∈ Ω,
∂ṽ∞
∂µ

= 0, x ∈ ∂Ω.

Since ‖ṽ∞‖L∞(Ω) = 1, we apply the Harnack inequality to obtain ṽ∞ > 0 in Ω.
Thus, it follows from the Krein-Rutman theorem that λ1(dv,−m−cb(x)(r−d)/a) =
0. This is a contradiction, and thus we derive the desired result. This completes
the proof. �

We next investigate the asymptotic behavior of any positive solution to (1.4)
with p(u, v) = u as du goes to 0, and a special numerical simulation example is
presented in Figure 6.

Theorem 4.4. Assume that r/(1 + km) > d, m > 0 and

r

1 + k(m+ cb∗(r − d)/a)
> d+ b∗

(
m+

cb∗(r − d)

a

)
. (4.2)



22 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0

0.005

0.01

0.015

0.02

0.025

u

du=105

du=106

du=107

du=108

du=109

du=1010

du=1011

du=1012

0 0.5 1 1.5 2 2.5 3 3.5 4
x

1.01

1.02

1.03

1.04

1.05

1.06

1.07

1.08

1.09

v

du=105

du=106

du=107

du=108

du=109

du=1010

du=1011

du=1012

Figure 5. Fix Ω = (0, 4), b(x) = 1+x2, u0(x) = 0.01+0.01 sin(5x)
and v0(x) = 1 + sin(5x). By choosing dv = 1, r = 24, k = 1, d = 5,
a = 1, m = 1, c = 0.6, each curve represents the solution (u, v) of
(1.4) with p(u, v) = u for large du.

Then for fixed dv > 0, any positive solution (u, v) of (1.4) with p(u, v) = u satisfies

(u, v)→
(1

a

( r

1 + kv∞
− d− b(x)v∞

)
, v∞

)
in L∞(Ω)× C1(Ω) as du → 0,

where v∞ is the unique positive solution of

−dv∆v∞ =
(
m+

c(r − d)b(x)

a
−
( crb(x)k

a(1 + kv∞)
+
cb2(x)

a
+ 1
)
v∞

)
v∞, x ∈ Ω,

∂v∞
∂µ

= 0, x ∈ ∂Ω.

(4.3)

Proof. By (4.2), it is clear that r/(1 + km) > d + mb∗. Hence, from Theorem
4.2(2) and (3), we see that (1.4) with p(u, v) = u admits at least one positive
solution for small du > 0. Let (udu,n

, vdu,n
) be any positive solution of (1.4) with

p(udu,n
, vdu,n

) = udu,n
and du = du,n, where du,n → 0 as n → ∞. Then for fixed

dv > 0, Lemma 2.7 shows that ‖vdu,n
‖W 2,p(Ω) is bounded for any p > 1. Hence,

it follows from the Sobolev embedding theorem that the sequence {vdu,n
}∞n=1 is

compact in C1(Ω). We may assume that vdu,n → v∞ ≥ m > 0 in C1(Ω) since
vdu,n ≥ m > 0 for all n ∈ N+. This result ensures that there exists some small
ε > 0 such that

0 < v∞ − ε ≤ vdu,n
≤ v∞ + ε in Ω (4.4)

for all large n.
By (4.4) and the equation for udu,n , we find that

−du,n∆udu,n ≤
( r

1 + k(v∞ − ε)
− d− audu,n − b(x)(v∞ − ε)

)
udu,n , x ∈ Ω.

A standard comparison argument yields udu,n
≤ Udu,n

in Ω for all large n, where

Udu,n is the unique positive solution of

−du,n∆Udu,n =
( r

1 + k(v∞ − ε)
− d− b(x)(v∞ − ε)− aUdu,n

)
Udu,n , x ∈ Ω,



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 23

∂Udu,n

∂µ
= 0, x ∈ ∂Ω.

Indeed, when (4.2) holds, we use Proposition 2.2(1) and Lemma 2.6 to obtain

λ1

(
du,n,−

r

1 + k(v∞ − ε)
+ d+ b(x)(v∞ − ε)

)
< 0.

Thus, the existence and uniqueness of Udu,n
is clear by [11, Theorems 3.5 and 3.7].

Moreover, from the above analysis, we find that the hypothesis (4.2) ensures that
r/(1 + k(v∞ − ε))− d− b(x)(v∞− ε) > 0 in Ω. Hence, by [17, Theorem 1.1(c)], we
see that

Udu,n
→ 1

a

( r

1 + k(v∞ − ε)
− d− b(x)(v∞ − ε)

)
in L∞(Ω) as n→∞.

This shows that

lim sup
n→∞

udu,n ≤
1

a

( r

1 + k(v∞ − ε)
− d− b(x)(v∞ − ε)

)
. (4.5)

On the other hand, by (4.4) and the equation for udu,n
, we find that

−du,n∆udu,n ≥
( r

1 + k(v∞ + ε)
− d− audu,n − b(x)(v∞ + ε)

)
udu,n , x ∈ Ω.

It follows from the standard comparison argument that udu,n
≥ Udu,n

in Ω for all
large n, where Udu,n

is the unique positive solution of

−du,n∆Udu,n
=
( r

1 + k(v∞ + ε)
− d− b(x)(v∞ + ε)− aUdu,n

)
Udu,n

, x ∈ Ω,

∂Udu,n

∂µ
= 0, x ∈ ∂Ω.

Note that the existence and uniqueness of Udu,n
can be verified as above. From

[17, Theorem 1.1(c)] it follows that

Udu,n
→ 1

a

( r

1 + k(v∞ + ε)
− d− b(x)(v∞ + ε)

)
in L∞(Ω) as n→∞.

This shows that

lim inf
n→∞

udu,n ≥
1

a

( r

1 + k(v∞ + ε)
− d− b(x)(v∞ + ε)

)
. (4.6)

Letting ε→ 0 in (4.5) and (4.6), we derive

udu,n
→ 1

a

( r

1 + kv∞
− d− b(x)v∞

)
in L∞(Ω) as n→∞.

This means that v∞ satisfies (4.3).
Since λ1(dv,−m−c(r−d)b(x)/a) < 0, [11, Theorems 3.5 and 3.7] imply that (4.3)

has a unique positive solution. Moreover, by (4.4), v∞ is a positive solution of (4.3).
Therefore, v∞ is the unique positive solution of (4.3). The proof is complete. �

Now we study the effect of the slow or fast movement of the predator individuals
on the profiles of positive solutions of (1.4) with p(u, v) = u. The following theorem
gives the asymptotic profiles of positive solutions as the predator diffusion coefficient
dv goes to 0, and a special numerical simulation example is presented in Figure 7.



24 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0.3

0.35

0.4

0.45

u

du=1 10

du=1

du=1 10 -1

du=8 10 -2

du=5 10 -3

du=1 10 -3

du=8 10 -4

du=5 10 -4

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0.105

0.11

0.115

0.12

0.125

0.13

0.135

0.14

0.145

0.15

0.155

v

du=1 10

du=1

du=1 10 -1

du=8 10 -2

du=5 10 -3

du=1 10 -3

du=8 10 -4

du=5 10 -4

Figure 6. Fix Ω = (0, 4), b(x) = 1 + x2, u0(x) = 0.01 +
0.01 sin(5x) and v0(x) = 1 + sin(5x). By choosing dv = 0.01,
r = 10, k = 1, d = 0.1, a = 20, m = 0.1, c = 0.01, each curve
represents the solution (u, v) of (1.4) with p(u, v) = u for small du.

Theorem 4.5. Assume that m > 0 and either d+ b∗m < r/(1 + km) < d+ b(x)m

and du ∈ (0, Du) or r/(1 + km) ≥ d + b(x)m and du > 0. Then for fixed du > 0,
any positive solution (u, v) of (1.4) with p(u, v) = u satisfies

(u, v)→ (u∞,m+ cb(x)u∞) in C1(Ω)× L∞(Ω) as dv → 0,

where u∞ is the unique positive solution of

−du∆u∞ =
ru∞

1 + k (m+ cb(x)u∞)
− (d+mb(x))u∞ − (a+ cb2(x))u2

∞, x ∈ Ω,

∂u∞
∂µ

= 0, x ∈ ∂Ω.

(4.7)

Proof. From Theorem 4.2(2) and (3), our assumption conditions ensure that (1.4)
with p(u, v) = u admits at least one positive solution for all dv > 0. Hence, we may
assume that (udv,n

, vdv,n
) is a positive solution of (1.4) with p(udv,n

, vdv,n
) = udv,n

and dv = dv,n, where dv,n → 0 as n → ∞. Then by Lemma 2.7, we find that
‖udv,n

‖W 2,p(Ω) is bounded for any p > 1. By the Sobolev embedding theorem,

{udv,n
}∞n=1 is compact in C1(Ω). Thus, we may assume that udv,n

→ u∞ ≥ 0 in

C1(Ω). This ensures that there exists some small ε > 0 such that

u∞ − ε ≤ udv,n
≤ u∞ + ε in Ω (4.8)

for all large n.
From the equation for vdv,n and (4.8), it follows from

−dv,n∆vdv,n ≤
(
m+ cb(x)(u∞ + ε)− vdv,n

)
vdv,n , x ∈ Ω,

∂vdv,n

∂µ
= 0, x ∈ ∂Ω,

−dv,n∆vdv,n
≥
(
m+ cb(x)(u∞ − ε)− vdv,n

)
vdv,n

, x ∈ Ω,
∂vdv,n

∂µ
= 0, x ∈ ∂Ω.



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 25

It follows from a comparison argument that V dv,n
≤ vdv,n

≤ V dv,n
in Ω for all large

n, where V dv,n
is the unique positive solution of

−dv,n∆V dv,n =
(
m+ cb(x)(u∞ + ε)− V dv,n

)
V dv,n , x ∈ Ω,

∂V dv,n

∂µ
= 0, x ∈ ∂Ω,

and V dv,n
is the unique positive solution of

−dv,n∆V dv,n
=
(
m+ cb(x)(u∞ − ε)− V dv,n

)
V dv,n

, x ∈ Ω,
∂V dv,n

∂µ
= 0, x ∈ ∂Ω.

From [11, Theorems 3.5 and 3.7 ], the uniqueness of V dv,n
or V dv,n

is clear since

λ1(dv,n,−m− cb(x)u∞± εcb(x)) < 0 for m > 0 and small ε > 0. We can apply [17,
Theorem 1.1(c)] to conclude that

V dv,n → m+ cb(x)(u∞ + ε) in L∞(Ω) as n→∞,
V dv,n

→ m+ cb(x)(u∞ − ε) in L∞(Ω) as n→∞.

Therefore, we let ε→ 0 to derive

vdv,n
→ m+ cb(x)u∞ in L∞(Ω) as n→∞.

This shows that u∞ is a non-negative solution of (4.7).
We next claim that u∞ is a positive solution of (4.7). Argue by a contradiction.

We may assume that there exists some point x0 ∈ Ω such that u∞(x0) = 0. Thus,
we apply the Harnack inequality to obtain that u∞ ≡ 0 in Ω. Hence, we find that
vdv,n → m uniformly in Ω as n → ∞. Denote ũdv,n = udv,n/‖udv,n‖L∞(Ω). From
the equation for udv,n

, it follows that

−du∆ũdv,n =
( r

1 + kvdv,n

− d− audv,n − b(x)vdv,n

)
ũdv,n , x ∈ Ω,

∂ũdv,n

∂µ
= 0, x ∈ ∂Ω.

Then we apply the standard elliptic regularity theory to derive ũdv,n
→ ũ∞ in

C1(Ω) as n→∞, where ũ∞ is a non-negative function and satisfies

−du∆ũ∞ =
( r

1 + km
− d− b(x)m

)
ũ∞, x ∈ Ω,

∂ũ∞
∂µ

= 0, x ∈ ∂Ω.

Since ‖ũ∞‖L∞(Ω) = 1, the Harnack inequality implies that ũ∞ > 0 in Ω. Conse-
quently, we apply the Krein-Rutman theorem to obtain that

λ1

(
du,−

r

1 + km
+ d+ b(x)m

)
= 0,

a contradiction. Therefore, u∞ > 0 in Ω.
Since our assumption ensure that λ1 (du,−r/(1 + km) + d+ b(x)m) < 0, and

(4.7) is equivalent to

−du∆u∞ =
( r

1 + km
− d−mb(x)

)
u∞

−
( crkb(x)

(1 + km)(1 + km+ ckb(x)u∞)
+ a+ cb2(x)

)
u2
∞, x ∈ Ω,

∂u∞
∂µ

= 0, x ∈ ∂Ω,



26 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

we derive from [11, Theorems 3.5 and 3.7] that the above equation has a unique
positive solution. Thus (4.7) has a unique positive solution u∞. The proof is
complete. �

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0

0.5

1

1.5

2

2.5

u

10-3

dv=1 10 -2

dv=1 10 -3

dv=1 10 -4

dv=1 10 -5

dv=1 10 -6

dv=1 10 -7

dv=1 10 -8

dv=1 10 -9

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0.01

0.01005

0.0101

0.01015

0.0102

0.01025

0.0103

0.01035

0.0104

v

dv=1 10 -2

dv=1 10 -3

dv=1 10 -4

dv=1 10 -5

dv=1 10 -6

dv=1 10 -7

dv=1 10 -8

dv=1 10 -9

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0.32

0.34

0.36

0.38

0.4

0.42

0.44

0.46

u

dv=1 10 -1

dv=1 10 -2

dv=1 10 -3

dv=1 10 -4

dv=1 10 -5

dv=1 10 -6

dv=1 10 -7

dv=1 10 -8

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0.32

0.34

0.36

0.38

0.4

0.42

0.44

0.46

u

dv=1 10 -1

dv=1 10 -2

dv=1 10 -3

dv=1 10 -4

dv=1 10 -5

dv=1 10 -6

dv=1 10 -7

dv=1 10 -8

Figure 7. Fix Ω = (0, 4), b(x) = 1+x2, u0(x) = 0.01+0.01 sin(5x)
and v0(x) = 1 + sin(5x), each curve represents the solution (u, v)
of (1.4) with p(u, v) = u for small dv by choosing du = 0.01, r = 5,
k = 2, d = 4.84, a = 20, m = 0.01, c = 0.1 in the first line, and by
choosing du = 1, r = 10, k = 1, d = 0.1, a = 20, m = 0.1, c = 0.01
in the second line.

As dv goes to ∞, the asymptotic behavior of positive solutions of (1.4) with
p(u, v) = u reads as follows, and a special numerical simulation example is presented
in Figure 8.

Theorem 4.6. Assume that m > 0 and either d+ b∗m < r/(1 + km) < d+ b(x)m

and du ∈ (0, Du) or r/(1 + km) ≥ d + b(x)m and du > 0. Then for fixed du > 0,
any positive solution (u, v) of (1.4) with p(u, v) = u satisfies

(u, v)→ (u∞, v∞) in C1(Ω)× C1(Ω) as dv →∞,

where

v∞ =
1

|Ω|

∫
Ω

(m+ cb(x)u∞)dx



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 27

and u∞ is a positive solution of

−du∆u∞ =
( r

1 + k
|Ω|
∫

Ω
(m+ cb(x)u∞)dx

− d

− b(x)

|Ω|

∫
Ω

(m+ cb(x)u∞)dx
)
u∞ − au2

∞, x ∈ Ω,

∂u∞
∂µ

= 0, x ∈ ∂Ω.

Since the proof of the above theorem is similar to that of Theorem 4.3, we omit
it.

0 0.5 1 1.5 2 2.5 3 3.5 4
x

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

u

dv=1 105

dv=1 106

dv=1 107

dv=1 108

dv=1 109

dv=1 1010

dv=1 1011

dv=1 1012

0 0.5 1 1.5 2 2.5 3 3.5 4
x

-0.5

0

0.5

1

1.5

2

2.5

v

dv=1 105

dv=1 106

dv=1 107

dv=1 108

dv=1 109

dv=1 1010

dv=1 1011

dv=1 1012

0 0.5 1 1.5 2 2.5 3 3.5 4
x

-0.5

0

0.5

1

1.5

2

2.5

3

u

dv=1 106

dv=1 107

dv=1 108

dv=1 109

dv=1 1010

dv=1 1011

dv=1 1012

dv=1 1013

0 0.5 1 1.5 2 2.5 3 3.5 4
x

0

0.5

1

1.5

2

2.5

v

dv=1 106

dv=1 107

dv=1 108

dv=1 109

dv=1 1010

dv=1 1011

dv=1 1012

dv=1 1013

Figure 8. Fix Ω = (0, 4), b(x) = 1 + x2, u0(x) = 0.01 +
0.01 sin(5x) and v0(x) = 1 + sin(5x), each curve represents the
solution (u, v) of (1.4) with p(u, v) = u for large dv by choosing
du = 1, r = 24, k = 1, d = 5, a = 1, m = 1, c = 0.6 in the first
line, and by choosing du = 1, r = 24, k = 1, d = 9, a = 1, m = 1,
c = 0.6 in the second line.

Finally, we consider the case m < 0 to investigate the effects of spatial diffusion
and environmental heterogeneity on the dynamics of (1.3). By Proposition 2.2,
the following properties of λ1 (dv,−m− c(r − d)b(x)/a) hold with respect to the
diffusion rate of the predator individuals.

Lemma 4.7. Assume that r > d and m < 0. Then

(1) λ1 (dv,−m− c(r − d)b(x)/a) is strictly monotone increasing with respect to
dv > 0.



28 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

(2) If m ≤ −cb∗(r− d)/a, then λ1(dv,−m− c(r− d)b(x)/a) > 0 for all dv > 0.

(3) If −cb∗(r−d)/a < m < −cb(x)(r−d)/a, then there is a unique value Dv =
Dv(r,m, d, a, c, b(x)) ∈ (0,∞) such that λ1(dv,−m − c(r − d)b(x)/a) < 0
for dv ∈ (0, Dv) and λ1(dv,−m− c(r − d)b(x)/a) > 0 for dv ∈ (Dv,∞).

(4) If −cb(x)(r − d)/a ≤ m < 0, then λ1(dv,−m− c(r − d)b(x)/a) < 0 for all
dv > 0.

In view of Theorems 2.4(2) and 2.11(2) and Lemma 4.7, the dynamic behavior
of (1.3) with p(u, v) = u with respect to the diffusion rate dv is as follows.

Theorem 4.8. Let p(u, v) = u. Assume that r > d and m < 0. Then

(1) If m ≤ −cb∗(r − d)/a, then ((r − d)/a, 0) is globally asymptotically stable
for all dv > 0.

(2) If −cb∗(r − d)/a < m < −cb(x)(r − d)/a, then ((r − d)/a, 0) is globally
asymptotically stable for each dv ∈ (Dv,∞); ((r − d)/a, 0) is unstable and
(1.3) admits at least one positive steady state for each dv ∈ (0, Dv).

(3) If −cb(x)(r−d)/a ≤ m < 0, then ((r−d)/a, 0) is unstable and (1.3) admits
at least one positive steady state for all dv > 0.

5. Summary and discussion

Motivated by some recent experimental field study and mathematical model anal-
ysis on the fear effect of prey, we proposed a reaction-diffusion system to demon-
strate the impact of fear cost on the population dynamics of prey. The novelty lies
in the incorporation of fear cost, spatial diffusion and environmental heterogene-
ity. We have theoretically analyzed the dynamics of the model and obtained some
insights on how fear cost, spatial diffusion and environmental heterogeneity affects
the population dynamics.

Our theoretical results established in section 3 indicate that when the predator
individuals are generalists, the cost of fear makes the prey more likely to become
extinct. Moreover, our results also indicate that the impact of fear cost on the
dynamics of the model is closely related to the functional response. For the case
with the linear functional response, it seems that only the effect described above is
observed. For the case with the Holling-type II functional response, we find that
the cost of fear will not only affect the existence of positive steady states but also
change the direction of steady-state bifurcation. For the case with the Beddington-
DeAngelis functional response, a more profound effect is observed: when the cost of
fear is ignored, a positive steady state always exists for any r > λ1 (du, d+ b(x)/f)
even though the predator is very strong; while when the cost of fear is considered,
(0,m) will become globally asymptotically stable, and no positive steady state
exists. Another effect of fear cost is that when the competition of the predator is
strong, the cost of fear makes the total population of the prey drop form (r−d)/a·|Ω|
to ((r/(1 + km)− d)/a,m) · |Ω|. This implies that the cost of fear not only makes
the prey more difficult to survive, but also reduce the total population of the prey
even though the prey population survives.

Our theoretical results established in section 4 indicate that spatial diffusion and
environmental heterogeneity significantly influence the dynamics of (1.3). When the
predator individuals are generalists, the predator can always invade (i.e., (0,m) is
globally asymptotically stable) for any diffusion rate of the prey if the birth rate
of the prey is small; the predator can invade for large diffusion rate of the prey



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 29

and will coexist with the prey for small diffusion rate of the prey if the birth rate
of the prey is intermediate; the predator and the prey can always coexist for any
diffusion rate of the prey if the birth rate of the prey is large. When the predator
individuals are specialists, the predator can rarely invade (i.e., ((r − d)/a, 0) is
globally asymptotically stable) for any diffusion rate of the predator if the death
rate of the predator is large; the predator can rarely invade for large diffusion rate of
the predator and will coexist with the prey for small diffusion rate of the predator if
the death rate of the predator is intermediate; the predator and the prey can always
coexist for any diffusion rate of the predator if the death rate of the predator is
small. These results shows that for the generalist predator, predator invasion or
coexistence with the prey mainly depends on the diffusion rate and the birth rate
of the prey; while for the specialist predator, predator invasion or coexistence with
the prey mainly depends on the diffusion rate and the death rate of the predator.

We next compare the findings here with those in [27] where system (1.2) were
studied where the coefficients are positive constants. We find that the analytical
methods and main results are quiet different from those in [27]. In particular,
spatially dependent coefficients bring the theoretical analysis lots of difficult and
challenging. In terms of results, we summarize and highlight according to various
functional response. When the linear functional response is adopted, consideration
of spatial heterogeneity can lead to spatial pattern formation. When the Holling-
type II or Beddington-DeAngelis functional response is adopted, we explicitly es-
tablish the sufficient conditions of spatial pattern formation, rather than just by
numerical simulation as in [27]. More importantly, the multiplicity and uniqueness
of non-constant positive steady states to (1.3) are established by choosing differ-
ent sets of parameters. In addition, by comparing our results with those results
obtained in [27], we find that the cost of fear has obvious impact on the dynamics
of (1.2) for m < 0, but has no qualitative impact on the dynamics of (1.3). This
indicates that whether the cost of fear affects the dynamics (especially for m < 0)
may be dependent on spatial heterogeneity.

Our theoretical results suggest that factors such as the level of fear, functional
responses, environmental heterogeneity and movement of individuals play vital but
subtle roles in the dynamics of (1.3). Therefore, a good understanding of these
factors could be helpful in understanding the dynamic behavior of the model and
be important in designing effective species conservation measures.

Acknowledgments. S. Li was supported by the NSF of China (11901446, 11971369,
12071270), by the Postdoctoral Science Foundation of China (2021T140530), by the
Young Talent fund of University Association for Science and Technology in Xi’an
(095920201325), and by the Fundamental Research Funds for the Central Universi-
ties (JB210705). Y. Xiao was supported by the NSF of China (11631012,11571273).
Y. Dong was supported by the NSF of China (11801431), by the Postdoctoral
Science Foundation of China (2018M631133), and by the Young Talent fund of
University Association for Science and Technology in Shaanxi China (20190509).

References

[1] L. J. S. Allen, B. M. Bolker, Y. Lou, A. L. Nevai; Asymptotic profiles of the steady states for

an SIS epidemic reaction-diffusion model, Discrete Contin. Dyn. Syst., 21 (2008) 1-20.
[2] I. Averill, K. Y. Lam, Y. Lou; The role of advection in a two-species competition model: a

bifurcation approach, Mem. Amer. Math. Soc., 245 (2017), 1161.



30 S. LI, Y. XIAO, Y. DONG EJDE-2021/70

[3] R. S. Cantrell, C. Cosner; Spatial Ecology via Reaction-Diffusion Equations, Series in Math-

ematical and Computational Biology, Wiley, Chichester, UK, 2003.

[4] M. G. Crandall, P. H. Rabinowitz; Bifurcation from simple eigenvalues, J. Funct. Anal. 8
(1971), 321-340.

[5] S. Creel, D. Christianson; Relationships between direct predation and risk effects, Trends

Ecol. Evol., 23 (2008), 194-201.
[6] E. N. Dancer, Y. H. Du; Effects of certain degeneracies in the predator-prey model, SIAM.

J. Math. Anal., 34 (2002), 292-314.

[7] A. Das, G. P. Samanta; Modeling the fear effect on a stochastic prey-predator system with
additional food for the predator, J. Phys. A Math. Theor., 51 (2018) 465601.

[8] Y. H. Du, Y. Lou; Some uniqueness and exact multiplicity results for a predator-prey model,

Trans. Amer. Math. Soc., 349 (1997), 2443-2475.
[9] Y.H. Du, R. Peng, M. X. Wang; Effect of a protection zone in the diffusive Leslie predator-

prey model, J. Differential Equations, 246 (2009), 3932-3956.
[10] Y. H. Du, J. P. Shi; Allee effect and bistability in a spatially heterogeneous predator-prey

model, Trans. Amer. Math. Soc., 359 (2007), 4557-4593.

[11] J. M. Fraile, P. Koch Medina, J. López-Gómez, S. Merino; Elliptic eigenvalue problems and
unbounded continua of positive solutions of a semilinear elliptic equation, J. Differential

Equations, 127 (1996), 295-319.

[12] I. Hanski, L. Hansson, H. Henttonen; Specialist predators, generalist predators, and the mi-
crotine rodent cycle, J. Anim. Ecol., 60 (1991), 353-367.

[13] M. P. Hassell, R. M. May; Generalist and specialist natural enemies in insect predator-prey

interactions, Br. Ecol. Soc., 55 (1986), 923-940.
[14] L. Li; Coexistence theorems of steady-states for predator-prey interacting systems, Trans.

Amer. Math. Soc. 305 (1988) 143-166.

[15] S. B. Li, J. H. Wu, Y. Y. Dong; Uniqueness and stability of positive solutions for a diffusive
predator-prey model in heterogeneous environment, Calc. Var. Partial Differential Equations,

(2019) 58: 110.
[16] S. L. Lima; Predators and the breeding bird: behavioural and reproductive flexibility under

the risk of predation, Biol. Rev., 84 (2009), 85-513.

[17] Y. Lou; On the effects of migration and spatial heterogeneity on single and multiple species,
J. Differential Equations, 223 (2006), 400-426.

[18] S. Mondal, A. Maiti, G. P. Samanta; Effects of fear and additional food in a delayed predator-

prey model, Biophys. Rev. Lett., 13 (2018), 157-177.
[19] E. L. Preisser, D. I. Bolnick; The many faces of fear: comparing the pathways and impacts

of nonconsumptive predator effects on prey populations, PloS One 3 (2008) e2465.

[20] W. H. Ruan, W. Feng; On the fixed point index and multiple steady-state solutions of reaction-
diffusion systems, Differential Integral Equations 8 (1995) 371-391.

[21] S. K. Sasmal; Population dynamics with multiple Allee effects induced by fear factors-A

mathematical study on prey-predator interactions, Appl. Math. Model., 64 (2018), 1-14.
[22] J. P. Shi; Persistence and bifurcation of degenerate solutions, J. Funct. Anal. 169 (1999)

494-531.
[23] W. E. Snyder, A. R. Ives; Interactions between specialist and generalist natural enemies:

parasitoids, predators, and pea aphid biocontrol, Ecol. Soc. Am., 84 (2003), 91-107.

[24] J. Wang, Y. L. Cai, S. M. Fu, W. M. Wang; The effect of the fear factor on the dynamics of
a predator-prey model incorporating the prey refuge, Chaos 29 (2019), 083109.

[25] X. Y. Wang, L. Y. Zanette, X. F. Zou; Modelling the fear effect in predator-prey interactions,
J. Math. Biol., 73 (2016), 1179-1204.

[26] X. Y. Wang, X. F. Zou; Modeling the fear effect in predator-prey interactions with adaptive

avoidance of predators, Bull. Math. Biol., 79 (2017), 1325-1359.

[27] X. Y. Wang, X. F. Zou; Pattern formation of a predator-prey model with the cost of anti-
predator behaviors, Math. Biosci. Eng., 15 (2018), 775-805.

[28] X. Y. Wang, X. F. Zou; On a predator-prey system with digestion delay and anti-predation
strategy, J. Nonlinear Sci., 30 (2020), 1579-1605.

[29] Y. Yamada; Stability of steady-states for prey-predator diffusion equations with homogeneous

Dirichlet conditions, SIAM J. Math. Anal., 21 (1990), 327-345.

[30] L. Y. Zanette, A. F. White, M. C. Allen, M. Clinchy; Perceived predation risk reduces the
number of offspring songbirds produce per year, Science, 334 (2011), 1398-1401.



EJDE-2021/70 DIFFUSIVE PREDATOR-PREY MODELS WITH FEAR EFFECT 31

Shanbing Li (corresponding author)

School of Mathematics and Statistics, Xidian University, Xi’an 710071, China

Email address: lishanbing@xidian.edu.cn

Yanni Xiao

School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China
Email address: yxiao@mail.xjtu.edu.cn

Yaying Dong

School of Science, Xi’an Polytechnic University, Xi’an 710048, China
Email address: dongyaying@xpu.edu.cn


	1. Introduction
	2. Solution and equilibria of system (1.3)
	2.1. Long-time behavior of solutions to (1.3)
	2.2. Existence of positive solutions to (1.4)

	3. Effect of fear
	3.1. Linear functional response
	3.2. Holling-type II functional response
	3.3. Beddington-DeAngelis functional response

	4. Effect of spatial diffusion and environmental heterogeneity
	5. Summary and discussion
	Acknowledgments

	References