Abstract: This paper concerned with a diffusive predator–prey model with fear effect. First, some basic dynamics of system is analyzed. Then based on stability analysis, we derive some conditions for stability and bifurcation of constant steady state. Furthermore, we derive some results on the existence and nonexistence of nonconstant steady states of this model by considering the effect of diffusion. Finally, we present some numerical simulations to verify our theoretical results. By mathematical and numerical analyses, we find that the fear can prevent the occurrence of limit cycle oscillation and increase the stability of the system, and the diffusion can also induce the chaos in the system.
Keywords: diffusion, stability, fear effect, predator–prey model.
Spatiotemporal dynamics of a diffusive predator–prey model with fear effect*
Vilniaus Universitetas
Esta obra está bajo una Licencia Creative Commons Atribución 4.0 Internacional.
Recepción: 08 Diciembre 2021
Aprobación: 14 Mayo 2022
Since Lotka [11] and Volterra [17] proposed famous Lotka–Volterra equations, the construction and study of models for the population dynamics of predator–prey interactions has been an important topic in theoretical ecology. According to different background, researchers have proposed many types of predator–prey models, and rich dynamics of these systems have been investigated extensively [6, 8, 18, 21]. In the wild, it is easy to observe that the reduction of prey is due to the direct killing of predators, which is reflected by functional responses in the predator–prey model such as Holling type and Beddington–DeAngelis [1, 7, 9, 10, 16, 24].
However, a new study suggested that the behavior of the prey can be changed by the predator, and it could have a greater impact than direct killing. In fact, Zanette et al. [22] found that the offspring production of the song sparrows reduced by 40% because of the fear from predator. To model the fear effect in predator–prey interactions, Wang et al. [19] proposed a predator–prey model as follows:
where 
is the birth rate of the prey, 
is the natural death rate of the prey, a represents the death rate due to intraspecies competition. The parameter 
refers to the level of fear, which reflects the reduction of prey growth rate due to the antipredator behavior. With the increase of 
, the growth rate of prey decreases. In [19], the authors consider that the functional response 
is the linear 
or the Holling type II 
. Their theoretical results show that fear effect could improve the stability of the predator–prey system.
It is considered that the trait effect has reduced the growth rate of the prey due to fear of predators, and the prey has been subjected to a strong Allee effect caused by mating during reproduction. Inspired by this idea, [14] considered a predator–prey model with the trait effect that reduced the growth rate of the prey due to fear of predators, and the prey has been subjected to a strong Allee effect caused by mating during reproduction. Their results showed that the fear effect does not affect the stability of the equilibria, but with the increasing of the cost of fear, the equilibrium density of predator decreases. Sasmal and Takeuchi [15] studied a predator–prey model that incorporates fear effect due to the presence of predator and group defense. Wang et al. [23] investigated a predator– prey model incorporating the fear of predators and a prey refuge, and they found that the fear effect can not only reduce the population density of predator, but also stabilize the system by excluding the existence of periodic solutions. Here, we remark that all models in these papers did not consider the factor of diffusion.
It should be pointed out that in real life, species are heterogeneous in space, so individuals tend to migrate to areas with low population density, which will increase the possibility of survival. Hence, some researchers considered reaction–diffusion predator– prey model by incorporating the fear effect into prey population. Niu et al. [4] investigated a diffusive predator–prey model with the fear effect. Taking the mature delay as bifurcation parameter, they found that the delay can induce Hopf and Hopf–Hopf bifurcations. Wang and Zou [20] proposed and analyzed a reaction–diffusion–advection predator–prey model. [3] investigated a diffusive predator–prey model with fear effect. Their results show that for the Holling type II predator functional response case, there exist no nonconstant positive steady states for large conversion rate.
Motivated by these pioneer work and note that none of the above mentioned models considered the Holling III functional response, we are led to consider a diffusive predator-prey model as follows:
(1)where 
denote the density of the prey and the predator at location 
and time
, respectively. r is the birth rate of prey, 
is the natural death rate of prey, 
represents the death rate due to intraspecies competition. The parameter 
reflects the level of fear, which drives antipredator behaviours of the prey. 
is Holling type-III function (see [5]). The parameter 
is the death rate of predator. 
is a bounded region with smooth boundary 
, and 
denotes the outward normal vector to the boundary 
. The homogeneous Neumann boundary condition indicates that there is no population flow across the boundary.
We also assume that
, and 
is defined by
In this paper, our goal is to investigate the dynamical properties of (1) such as global existence of the solutions, stability and bifurcation of the constant steady state. In addition, we will use energy estimates to obtain of the dynamic and steady state solutions and so to discuss the nonexistence and existence of spatial patterns.
Our paper is organized as follows. In Section 2, we study some basic dynamics of the system. In Section 3, we obtain the stability and bifurcation of the equilibria. In Section 4, we investigate the nonexistence and existence of the nonconstant steady state. In Section 5, numerical results are presented to verify the theoretical results.
In this section, we discuss some basic dynamics of system (1) including the existence of solution and the priori bound of the solution.
First, we let 
be the Lebesgue measure of 
and denote
Theorem 1. For system (1), the following conclusions are true:
-
If
then system (1) admits a unique solution 
and 
-
If
, then 
-
If
, then any solution 
of (1) satisfies
In addition,
where
, and 
depends on 
, 
and
; -
If
then
Proof. (i) Define
then 
and 
in 
Consequently, system (1) is a mixed quasimonotone system. Consider the following ordinary differential equation model:
(2)where 
Let 
be the unique solution of system (2). Then 
and 
are the lower solution and upper solution of system (1). Thus, according to the [13, Thm. 8.3.3], system (1) has a unique globally defined solution 
, which satisfies
The strong maximum principle ensures that 
(ii) The first equation of system (2) implies that
Obviously, 
leads 
Consequently, 
(iii) It is noted that
Thus, by the comparison principle, one have
The maximum principle ensures that 
for all 
Let 
then
(3)
Multiplying both sides of Eq. (3) by
then combining with Eq. (4), we obtain
Noting that 
proved above, we have 
Thus
where 
The integration of inequity (5) results in
Figure 1
Invariant region Rα for system (1).
which means that
From [2, Thm. 3.1] we have 
where 
depends on 
and 
As a result, there is a 
such that 
(iv) Obviously, 
proved above that if 
then 
uniformly on 
Theorem 2. The trapezoidal region
is a positively invariant region for system (1) (see Fig. 1).
Proof. The reaction kinetics do not point out of 
along 
and 
Setting
and denoting the outward normal to 
alon the line 
then denoting 
one obtain
as 
Consequently, 
is an invariant region.
Theorem 3. For system (1), the following conclusions hold:
-
If
, then system (1) has only the trivial constant solution 
; -
If
, then system (1) has a predator-free constant steady state solution 
denoting the extinction of predator, while system has no positive constant steady state solution -
If
then system (1) has a unique positive constant steady state solution.
Proof. Obviously, (i) and (ii) hold. The positive constant steady state solution 
satisfies
(6)From the second equation of (6) we have 
Then according to the first equation of (6), we obtain
(7)where
Clearly, 
Therefore, Eq. (7) has at most a positive root as that 
, which means that 
Therefore, we have the conclusion.
Recall tha t
are the eigenvalues of the Laplace operator 
on 
under homogeneous Neumann boundary condition, and 
is the space of eigenfunctions corresponding to 
in 
Let 
be an orthonormal basis of 
and 
Then
Assume that 
is a constant solution of system (1), then we have
with domain 
where
and
For each 
is invariant under the operator 
, and 
is an eigenvalue of 
on 
if and only if 
is an eigenvalue of 
for all 
The direct calculation shows
(8)where
Theorem 4.
-
If
, then 
is globally asymptotically stable. -
If
then 
is globally asymptotically stable. -
If
then 
is locally asymptotically stable.
Proof. (i) For 
, the corresponding characteristic equation is
Clearly, we obtain
Hence, if 
, then 
is locally asymptotically stable. Note that there is no other constant steady states in this case. This means that 
is indeed globally asymptotically stable.
(ii) For 
, the corresponding characteristic equation is
Obviously,
Consequently, if 
then 
is locally asymptotically stable. In fact, 
is globally asymptotically stable.
It follows from Theorem 1 that 
so for 
By the second equation of (1) we have
Therefore, 
and there exists 
such that
Then by first equation of (1), one have
Then we obtain that 
Combining with 
allows us to derive
Hence, 
is globally asymptotically stable when 
(iii) For the positive steady state 
Hence, the corresponding characteristic equation is
Obviously,
(9)All roots of (9) have negative real parts if
Therefore, the positive constant steady state 
is locally asymptotically stable when condition (10) holds.
Remark 1. Theorems 3 and 4 show that when 
, system has only trivial constant solution 
, and it is globally asymptotically stable; when 
increases and enter the interval 
, 
lose s its stability to a predator-free constant to a positive steady state 
and when 
further passes 
loses its stabilityto a positive steady state
. We can conclude that as the parameter r increases, the model experiences two bifurcations of constant steady state.
Remark 2. Obviously, the conditions of Theorem 4 are independent of the diffusion. Consequently, the conclusions of Theorem 4 are still valid for the corresponding ODE model. In addition, we can also conclude that the diffusion cannot destabilize the positive steady state
. Therefore, the PDE system (1) cannot occur Turing instability/bifurcation.
In this subsection, we will discuss the bifurcation of system (1). Let the parameters 
and 
be fixed, and be fixed, take 
as a bifurcation parameter.
Theorem 5.
-
If
holds, then spatially homogeneous Hopf bifurcation occurs. -
If
let 
be the largest positive integer such that 
In addition, we assume that 
whenever 
, and
(11)
Then system (1)undergoes spatially inhomogeneous Hopf bifurcation at 
for 
where 
Proof. (i) If 
holds, then 
and
Therefore, Eq. (8) has a pair of pure imaginary roots 
which means that spatially homogeneous Hopf bifurcation occurs.
(ii) From the assumption it follows that
for 
and 
In addition, 
for any
Clearly,
Obviously, if condition (11) holds, then 
Moreover, if 
then
Therefore, 
is nondecreasing with respect to 
Hence, when
Therefore, when 
is near 
Eq. (8) has a pair of conjugate eigenvalues
Clearly Ré 
As a result, Hopf bifurcation occurs at 
, which also means that system (1) has a family of inhomogeneous periodic solutions near 
.
In this section, we will discuss nonexistence and existence of nonconstant steady state of system (1). To this end, we consider the following elliptic system:
(12)To derive some priori estimates for nonnegative solutions of system (12), we need the following technical lemma [12].
Lemma 1 [Maximum principle]. Suppose that 
is a bounded domain in 
and 
is a weak solution of the inequalities
and if there is a constant 
such that 
then 
Lemma 2 [Harnack inequality]. Suppose that 
and 
is a positive classical solution to 
subject to the homogenuous Neumann boundary condition. Then there exists a positive constant 
such that
For the sake of discussion, we shall write 
Theorem 6 [Upper bounds]. Suppose that 
is a nonnegative solution of (12), then either 
is one of constant solutions 
and 
or for 
satisfies
(13)where
Proof. If there exists 
satisfying 
then by the strong maximum principle, 
and
Thus,
Otherwise, 
Further, from Lemma 1 we obtain that 
and by the strong maximum principle, we have 
for all 
Then
It can be obtained from the maximum principle that
Therefore,
Theorem 7. Let 
be a positiveconstant. Then for 
there exists two positive constants 
with 
depending possibly on 
such that any solutions 
of system (12)satisfies
Proof. We choose
, for any 
Next, we shall prove 
. Let
Thus,
Lemma 2 shows that there exists a positive constant 
such that
Hence, now it remains to prove that there exists 
such that
(14)Contrariwise, let us assume that (14) does not hold. Then there exists a sequence 
such that
By the regularity theory for elliptic equations, there exists a subsequence of 
, which will be denoted again by 
such that 
in 
as 
Note that 
and from (15) either 
Therefore, we have that
Also, 
satisfy (13), so do 
and 
. Letting 
we get that (u., v.) is a positive solution of (12). Therefore, by integrating Eq. (12) for 
and 
over 
, we have
(i) In this case,
then
and 
then
for sufficiently large 
So, we obtain a contradiction.
(ii) If 
using the first equation of (12). So, 
for large 
Thus
So, we have
for a sufficiently large 
, which is a contradiction. This completes the proof.
In this subsection, we show the nonexistence of positive steady state solutions when the diffusion coefficients 
and 
are large.
Theorem 8. For any fixed
there exists a positive constant 
such that if 
then (12) has no nonconstant solutions.
Proof. Assume that 
is nonnegative solution of (12). Denote
Obviously,
and
. For the purpose of discussions, let 
By the mean value theorem of bivariate functions, we have
Obviously, 
where 
Multiplying both sides of the first equation of (12) by 
and using Theorem 6, we get
(16)Applying Theorem 6 and by multiplying 
to the second equation in (12) and then integrating on 
, we have
(17)Using the Poincaré inequality,
where 
is the second eigenvalue of the Laplace operator 
under homogeneous Neumann boundary condition.
Combining (16) and (17) leads to
where
This implies that
then we can conclude that
To study the existence of nonconstant positive solutions, we use Leray–Schauder degree theory. Let 
and
Thus, (12) can be rewritten as
or equivalently,
(18)where 
represents the inverse of 
with the homogeneous Neumann boundary condition. From (18), by a direct computation, we have
Clearly,
where
Obviously, if 
then 
for all 
If
then 
has two positive roots as follows:
Set 
Obviusly,
Theorem 9. Assume that
and there exist 
such that
, and 
is odd. Then there exists at least one nonconstant solution of (12).
Proof. Let 
be defined in Theorem 8, and for 
Consider the following problem:
(19)It is easy to see that solving (12) is equivalent to find a fixed point of 
with 
is the unique constant solution of (19) for any 
By the definition of 
in Theorem 8, one have that 
is the only fixed point of 
.
Since 
and if (12) has no nonconstant positive solution, then we have
In addition, by the homotopy invariance of the topological degree,
which is a contradiction
In this section, we take some numerical simulations to discuss the effect of diffusion and the cost of fear.
In order to discuss the effect of diffusion, in 
we assume the parameters values:
A direct calculation shows that system (1) has a positive steady state
According to Theorem 4, the positive steady state is unstable. Figure 2 shows that system (1) has a stable limit cycle around the positive steady state 
with the initial conditions 
However, if we change the diffusion rate of 
to be 
we find that the stable limit cycle is broken with the occurrence of spatial pattern (see Fig. 3), where the periodic pattern disappears and a strip pattern appears.
Figure 2
The positive steady state E∗ = (0.2, 0.0528) is unstable, and there exists a stable limit cycle with the initial values u0(x) = 0.2 + 4 10−4 cos(2x), v0(x) = 0.0528 + 5 10−4 cos(2x) and the diffusion rate d2 = 0.5.
Furthermore, if we vary the diffusion coefficient 
from 
then we find that system (1) is chaotic (see Fig. 4).
We further find that different initial conditions with the same diffusion rate 
can lead to different spatial patterns that can be stationary or periodic (Fig. 5).
Figure 3
The emergent stationary spatial pattern with the initial values u0(x) = 0.2 + 4 · 10−4 cos(2x), v0(x) = 0.0528 + 5 · 10−4 cos(2x) and the diffusion rate d2 = 1.
Figure 4
System (1) is chaotic with the initial values u0(x) = 0.2 + 4 · 10−4 cos(2x), v0(x) = 0.0528 + 5 · 10−4 cos(2x) and the diffusion rate d2 = 0.002.
Figure 5
Spatial patterns and spatially averaged population dynamics for different random perturbed initial conditions with the same diffusion rate d2 = 1.2.
Choose 
Calculations show that system (1) has a unique positive steady state 
According to Theorem 4, we observe that the positive steady state 
of system (1) is locally asymptotically stable, and the dynamic behaviors of system (1) is illustrated graphically in Fig. 6.
From above discussions we can obtain that fear can affect the stability of the positive steady state, and it can induce the Hopf bifurcation, which is different from the results found in [14, 19] with linear functional response (see Fig. 9). Figure 9 shows that there exists a threshold value k0 such that when 
system (1) has a periodic solution. But when k passes the threshold value, then system becomes stable.
Figure 6
The positive steady state E∗ = (0.0325, 0.2000) is locally asymptotically stable.
Figure 7
Hopf bifurcation of system (1) with k = 20.37.
Figure 8
The positive steady state v∗ with varying the cost of fear k.
Figure 9
The maximum and minimum of prey u and predator v with the cost of fear k varying in [0, 50].
If we choose 
, while other parameters do not change, according to Theorem 5, system (1) undergoes spatial homogeneous Hopf bifurcation (see Fig. 7). It is shown that system (1) has spatially homogeneous periodic solutions emerged from the positive steady state 
.
In addition, we find that the positive steady state can be changed by the different value of the cost of fear. Figure 8 shows that the positive steady state 
decreases with increasing of the cost of fear.
A diffusive predator–prey model with the fear effect is studied in our paper. We derive some basic dynamics of the system and give condition for the existence of the positive steady state. According to eigenvalue analysis method, we investigate the stability and bifurcation of the positive constant steady state. We also give some conditions for the nonexistence and existence of nonconstant solutions of the system.
Theorems 3 and 4 show that the birth rate of prey r can not only induce the static bifurcation, but also can induce saddle-node bifurcation.
Theorem 4 indicates that the diffusion can not induce the Turing instability/bifurcation. However, Theorem 5 provides that the diffusion can induce the inhomogeneous Hopf bifurcation, which can lead to the formation of spatial patterns. Furthermore, Theorem 9 shows that system (12) has at least one nonconstant positive solution under the effect of the diffusion. From Section 5.1 we can obtain that the different diffusion rate 
can lead to different spatial patterns, which can be periodic (Fig. 2), stationary (Fig. 3) and chaotic (Fig. 4). In addition, we also find that system has different spatial patterns with the different initial conditions (Fig. 5).
We further obtain that the fear effect can reduce the density of predator: with increasing the cost of fear, the density of predator population decreases at the positive steady state (see Fig. 8). From Section 5.2 it is obtained that the fear can prevent the occurrence of limit cycle oscillation and increase the stability of the system (see Fig. 9).
Figure 1
Invariant region Rα for system (1).
Figure 2
The positive steady state E∗ = (0.2, 0.0528) is unstable, and there exists a stable limit cycle with the initial values u0(x) = 0.2 + 4 10−4 cos(2x), v0(x) = 0.0528 + 5 10−4 cos(2x) and the diffusion rate d2 = 0.5.
Figure 3
The emergent stationary spatial pattern with the initial values u0(x) = 0.2 + 4 · 10−4 cos(2x), v0(x) = 0.0528 + 5 · 10−4 cos(2x) and the diffusion rate d2 = 1.
Figure 4
System (1) is chaotic with the initial values u0(x) = 0.2 + 4 · 10−4 cos(2x), v0(x) = 0.0528 + 5 · 10−4 cos(2x) and the diffusion rate d2 = 0.002.
Figure 5
Spatial patterns and spatially averaged population dynamics for different random perturbed initial conditions with the same diffusion rate d2 = 1.2.
Figure 6
The positive steady state E∗ = (0.0325, 0.2000) is locally asymptotically stable.
Figure 7
Hopf bifurcation of system (1) with k = 20.37.
Figure 8
The positive steady state v∗ with varying the cost of fear k.
Figure 9
The maximum and minimum of prey u and predator v with the cost of fear k varying in [0, 50].
