Electronic Journal of Differential Equations, Vol. 2023 (2023), No. 37, pp. 1–20.

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

BOUNDEDNESS, STABILITY AND PATTERN FORMATION FOR

A PREDATOR-PREY MODEL WITH SIGMOID FUNCTIONAL

RESPONSE AND PREY-TAXIS

ZHIHONG ZHAO, HUANQIN HU

Abstract. This article concerns the structure of the nonconstant steady states
for a predator-prey model of Leslie-Gower type with Sigmoid functional and

prey-taxis subject to the homogeneous Neumann boundary condition. The

existence of bounded classical global solutions is discussed in bounded do-
mains with arbitrary spatial dimension and any prey-taxis sensitivity coef-

ficient. The local stability of the homogeneous steady state is analyzed to

show that the prey-taxis sensitivity coefficient destabilizes the stability of the
homogeneous steady state when prey defends. Then we study the existence

and stability of the nonconstant positive steady state of the system over 1D

domain by applying the bifurcation theory and present properties of local
branches such as pitchfork and turning direction. Moreover, we discuss global

bifurcation, homogeneous steady state solutions, nonconstant steady states

solutions, spatio-temporal periodic solutions and spatio-temporal irregular so-
lutions which demonstrate the coexistence and spatial distribution of prey and

predator species. Finally, we perform numerical simulations to illustrate and

support our theoretical analysis.

1. Introduction

In the field of biomathematics, the study of predator-prey interactions is one
of the fundamental subject. It is well known that there are a large number of
factors that influence dynamics of the predator-prey models, such as the birth
rate, mortality, interspecific competition, food, infectious disease, and functional
responses etc. One of the typical functional responses is the Holling II type:

P (u) =
bu

bhu+ 1
,

where u is the density of the prey, the positive constants b and h denote the search
rate and the processing time, respectively. It applies to invertebrates [19] (insects
and parasitic species are included here) and can be seen that the predator search
rate b is constant. As the study progressed, Hessel [8] pointed out that many
invertebrate predators search actively when the density of prey increases, and that
the search efficiency of predators decreases when density of prey falls below a certain
threshold, e.g., Coccinella septempunctata and first instar aphid, Plea atomaria and

2020 Mathematics Subject Classification. 35B32, 35J67, 35K57, 92D25, 92D40.
Key words and phrases. Predator-prey model; prey-taxis; boundedness; steady states;

global bifurcation; pattern formation.
©2023. This work is licensed under a CC BY 4.0 license.

Submitted January 6, 2023. Published May 4, 2023.

1



2 ZH. ZHAO, H. HU EJDE-2023/37

small Aedes larvae. This means that the search rate of many invertebrate predators
depends on the density of prey, i.e., b = au

u+g , where the positive constants a and

g denote the maximum achievable search rate and the half-saturation constant,
respectively. In this way we obtain the Sigmoid functional response [6], which is
the so-called generalized Holling type III:

P (u) =
au2

ahu2 + u+ g
. (1.1)

Turing [29] proposed a reaction diffusion system for the chemical basis of mor-
phogenesis, which is one of the most important mechanisms in theoretical biology.
Since the spatial heterogeneity of the environment affects predator-prey dynamics,
such system has been extensively investigated since 1970s. In particular, Turing bi-
furcation, Hopf bifurcation, pattern formation and traveling wave have been widely
studied [10, 25, 26, 28, 31, 32], where spatial diffusion refers to the variation of
species density according to some process (e.g., physical diffusion or random walk
diffusion).

Kareiva and Odell [13] pointed out that predators move toward areas of high
prey density to improve predation efficiency when the area is restricted, creating a
predator aggregation phenomenon. This behavior can be considered as the move-
ment of the predators in the direction of the prey density gradient, which is called
prey-taxis. Moreover, the predator-prey model with prey-taxis is different from the
previous models and has rich spatial-temporal dynamics. Global existence, bifurca-
tion analysis and pattern formation of prey-taxis systems have been widely studied
[2, 9, 16, 20, 24, 27, 33, 34, 36, 37, 38, 39] and the references therein.

In this article, we consider a predator-prey system with prey-taxis:

∂u

∂t
= d1∆u+ ru(1− u

K
)− P (u)v, in (0, T )× Ω,

∂v

∂t
= d2∆v − ξ∇(χ(v)∇u) + sv(1− e v

u
), in (0, T )× Ω,

∂u

∂ν
=
∂v

∂ν
= 0, on (0, T )× ∂Ω,

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

(1.2)

where u(x, t), v(x, t) are the density of the prey and predator at position x and time
t respectively, the habitat of both species Ω is a bounded domain with a smooth
boundary ∂Ω in Rn (n ≥ 1 is a positive integer), ν is the outward unit normal
vector of the boundary ∂Ω, the initial data u0(x), v0(x) are nonnegative continuous
functions and the positive constants d1, d2 are diffusion coefficients. The parameters
r and s represent the intrinsic growth rate. The prey population grows logistically
with the carrying capacity K. The predator consumes the prey according to the
sigmoid functional response P (u). The term ev

u is known as the Leslie-Gower term
[17, 18], and e is the number of prey required to support one predator when v equals
u
e . The parameter ξ denotes the prey-taxis sensitivity coefficient, χ(v) denotes the
prey-tactic cross diffusion and the term ξχ(v)∇u denotes the speed of the predator
moving in the direction of the density gradient of the prey. For the prey-taxis
sensitivity coefficient ξ, ξ > 0 indicates that the predator moves in the direction
of high prey density to improve predation efficiency, and ξ < 0 indicates that the
predator moves in the opposite direction of high prey density in order to avoid the
group defense of large amounts of prey, which is exemplified in nature by Japanese



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 3

bees that gather together to form a hot ball to kill bumblebees when faced with
bumblebees [3, 15]. More biological background can be found in [12, 30].

For convenience, we dimensionless size the model (1.2) as

ū =
u

K
, t̄ = rt, d̄1 =

d1

r
, ā = Kah, ḡ =

g

K
, ū0 =

u0

K
, v̄ =

av

r

d̄2 =
d2

r
, ξ̄ =

aKξ

r2
, χ̄(v̄) = χ(v), c̄ =

s

r
, δ̄ =

re

aK
, v̄0 =

av0

r
,

and ignore the bars on u, v and other parameters, then system (1.2) can be re-
expressed by

∂u

∂t
= d1∆u+ u(1− u)− u2v

au2 + u+ g
, in (0, T )× Ω,

∂v

∂t
= d2∆v − ξ∇(χ(v)∇u) + cv(1− v

δu
), in (0, T )× Ω,

∂u

∂ν
=
∂v

∂ν
= 0, on (0, T )× ∂Ω,

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

(1.3)

Throughout this paper we assume that

(H1) χ(v) ∈ C4[0,∞), χ(0) = 0 and there exists B > 0 such that χ(v) ≤ Bv for
any v ≥ 0 and x ∈ Ω.

In this article, we investigate the global existence of the classical solution and the
effect of the prey-taxis sensitivity coefficient on the stability, non-constant posi-
tive steady state and pattern formation of (1.3). Main results reveal that system
(1.3) has global bounded classical solutions in bounded domains, the prey-taxis
destabilizes the stability of the homogeneous steady state in the corresponding
ODE system or the PDE system with random diffusion, and a nonconstant positive
steady-state bifurcation occurs at the critical bifurcation point. Numerical simula-
tions are carried out to veirfy our theoretical results. We also find that the period
of the spatio-temporal periodic solutions becomes greater when the intrinsic growth
rate of predators increases within a certain range, which is a novel outcome of this
study.

The rest part of this paper is organized as follows. In Section 2, we prove the
global existence and boundedness of the classical solution. In Section 3, we perform
the stability of the homogeneous steady state by linearizing the model. Section 4 is
devoted to establishing the existence of a non-constant positive steady state in one-
dimensional space, the stability of this steady branch and the global bifurcation.
In Section 5, numerical simulations are presented. We discuss our results and raise
some interesting problems for the future study.

1.1. Global solutions. In this section, we discuss the existence of global solutions
of the initial-boundary value problem (1.3) of arbitrary spatial dimension with any
prey-taxis sensitivity coefficient.

Theorem 1.1. If u0(x) > 0, v0(x) ≥ 0, then the system (1.3) admits a unique
classical solution (u(x, t), v(x, t)) for all x ∈ Ω, t > 0. Furthermore, there exist
constants Ui > 0, Vi ≥ 0, i = 1, 2 depending on the initial values u0(x) and v0(x)
such that

U1 < u(x, t) < U2, V1 < v(x, t) < V2, x ∈ Ω, t > 0.



4 ZH. ZHAO, H. HU EJDE-2023/37

Proof. The existence and uniqueness of the solution to the initial-boundary value
problem (1.3) is easily obtained. For the global existence and the boundedness, we
apply the theory of invariant region in [35]. For the vector field(

u(1− u)− u2v

au2 + u+ g
, cv
(

1− v

δu

))
,

if there exists a region < = (U1, U2)×(V1, V2) in the (u, v) phase plane such that the
vector field points inward on the boundary of <, then < is a (positively) invariant
rectangle. If < is positively invariant to the vector field then the solution exists
globally and is bounded.

Let

U2 = max
{

1,max
x∈Ω

u0(x)
}
, V2 = max

{
δ,max
x∈Ω

v0(x)
}
,

U1 = min
{ g

g + δ
,min
x∈Ω

u0(x)
}
, V1 = min

{
0,min
x∈Ω

v0(x)
}
.

Clearly, the initial functions u0(x) and v0(x) are enclosed by the rectangle. And
for u = U1, V1 ≤ v ≤ V2, we have

u(1− u)− u2v

au2 + u+ g
≥ U1(1− U1)− U2

1 δ

g
≥ 0,

which means that the vector field crosses this line in the direction pointing to the
interior of <. For the other three sides, we verify using similar calculations.

Therefore, the region < is positively invariant with respect to the vector field
and the the proof is complete. �

2. Stability analysis of homogeneous steady state

In this section, we analyze the stability of the nontrivial homogeneous steady
state (u∗, v∗) by the characteristic equations. It can be seen that the corresponding
ODE system of system (1.3) is a predator-prey model of Leslie-Gower type with
Sigmoid functional response. Huang [11] investigated the types of equilibrium and
bifurcations of this ODE system for different values of parameters. By simple
calculation, system (1.3) has a boundary equilibrium (1, 0) and unique interior
equilibrium (u∗, v∗), where u∗ satisfies:

h(u) = au3 + (δ + 1− a)u2 + (g − 1)u− g = 0 (2.1)

and v∗ = δu∗. It is clear that h(0) = −g < 0, h(1) = δ > 0, then u∗ ∈ (0, 1). By
Cardano’s Method, (2.1) has a unique positive solution u∗,

u∗ = 3

√
−q

2
+
√

∆ + 3

√
−q

2
−
√

∆− δ + 1− a
3a

,

when ∆ = ( q2 )2 − (p3 )3 > 0, where

p =
3a(g − 1)− (δ + 1− a)2

3a2
, q =

2(δ + 1− a)3 − 27a2g − 9a(δ + 1− a)(g − 1)

27a3
.

We make the following assumption:

(H2) ∆ > 0 and

1− 2u∗ − (1− u∗)2(u∗ + 2g)

δu∗2
< c.



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 5

And let

f(u, v) = u(1− u)− u2v

au2 + u+ g
, g(u, v) = cv(1− v

δu
).

then we have the following results for the stability of the homogeneous steady state
(u∗, v∗) for ODE system.

Lemma 2.1. If (H2) holds, then the homogeneous steady state (u∗, v∗) is locally
asymptotically stable for system (1.3) without spatial variation.

Proof. Linearizing (1.3) without spatial variation at the steady state (u∗, v∗), then
we obtain Ut = AU , where

U =

(
u
v

)
, A =

(
f1 f2

g1 g2

)
, (2.2)

and

f1 = fu(u∗, v∗) = 1− 2u∗ − (1− u∗)2(u∗ + 2g)

δu∗2
, f2 = fv(u

∗, v∗) =
u∗ − 1

δ
< 0,

g1 = gu(u∗, v∗) = cδ > 0, g2 = gv(u
∗, v∗) = −c < 0.

Clearly, f1g2− f2g1 > 0, thus the locally asymptotically stability of (u∗, v∗) can be

guaranteed when f1 + g2 = 1− 2u∗ − (1−u∗)2(u∗+2g)
δu∗2 − c < 0 holds. �

Before developing our argument, let us set up the following notation.

• 0 = σ0 < σ1 < σ2 < . . . are the eigenvalues for the elliptic operator −∆ on
Ω under the homogeneous Neumann boundary condition.
• E(σi) and φij , j = 1, 2, . . . ,dimE(σi) be the eigenspace and eigenfunctions

corresponding to eigenvalue σi, i = 0, 1, 2 . . . , respectively.
• Xij := {c ·φij , c ∈ R2}, where {φij} are standard orthogonal bases in space
E(σi), for j = 1, 2, . . .dimE(σi).

• X :=
[
L2(Ω)

]2
, then X = ⊕+∞

i=1Xi and Xi =
⊕dimE(σi)

j=1 Xij .

Then the stability of (u∗, v∗) of system (1.3) is given in the following theorem.

Theorem 2.2. Assume (H2) holds, and denote

Qi(σi, 0) = d1d2σ
2
i − (d1g2 + d2f1)σi + f1g2 − f2g1, i ≥ 1, (2.3)

ξi :=
Qi(σi, 0)

χ(v∗)f2σi
, i ≥ 1 and ξ̂ := max

1≤i≤∞
ξi. (2.4)

Then the homogeneous steady state (u∗, v∗) is locally asymptotically stable for any

ξ > ξ̂ and (u∗, v∗) is unstable when ξ < ξ̂.

Proof. Linearizing system (1.3) at (u∗, v∗), we obtain

Ut = (A+D∆)U, D =

(
d1 0

−ξχ(v∗) d2

)
, (2.5)

where A and U are the same as (2.2). Assume that (ϕ,ψ) is the eigenfunction
corresponding to the eigenvalue µ of operator A+D∆, then they can be expressed
as (

ϕ
ψ

)
=

∑
0≤i<∞,1≤j≤dimE(σi)

(
aij
bij

)
φij(x). (2.6)

And the characteristic equation of (2.5) for an eigenvalue is

µ2 − Pi(σi)µ+Qi(σi, ξ) = 0,



6 ZH. ZHAO, H. HU EJDE-2023/37

where Pi(σi) = −(d1 + d2)σi + f1 + g2 < 0 and Qi(σi, ξ) = d1d2σ
2
i − (ξχ(v∗)f2 +

d2f1 + d1g2)σi + f1g2 − f2g1, i ≥ 1. By (2.3), we have

Qi(σi, ξ) = Qi(σi, 0)− ξχ(v∗)f2σi. (2.7)

It is easy to obtain that (u∗, v∗) is local asymptotic stability by Qi(0, ξ) =
f1g2 − f2g1 > 0 if i = 0. For i ≥ 1, we have:

(i) Qi(σi, 0) > 0 for all i ≥ 1. Then ξ̂ < 0 by (2.4). If ξ > ξ̂, then Qi(σi, 0)−
ξχ(v∗)f2σi > Qi(σi, 0)−ξ̂χ(v∗)f2σi = 0 for every i, which implies (u∗, v∗) is

local asymptotically stable by (2.7). If ξ < ξ̂, then Qi(σi, 0)−ξχ(v∗)f2σi <
0 for some i, which means (u∗, v∗) is unstable.

(ii) Qi(σi, 0) < 0 for some i ≥ 1. Then ξ̂ > 0 by (2.4). If ξ > ξ̂, then Qi(σi, 0)−
ξχ(v∗)f2σi > Qi(σi, 0)− ξ̂χ(v∗)f2σi = 0 for every i. Thus, (u∗, v∗) is local

asymptotically stable by (2.7). If ξ < ξ̂, then Qi(σi, 0) − ξχ(v∗)f2σi < 0
for some i. Thus (u∗, v∗) is unstable.

�

From Theorem 2.2, we can point out that the homogeneous steady state (u∗, v∗)
loses its stability when prey-taxis sensitivity coefficient is less than a certain thresh-

old value ξ̂. In what follows, we only discuss case (i): Qi(σi, 0) > 0 for all i ≥ 1 in
the proof of Theorem 2.2 while case (ii) can be treated by the same arguments.

3. Nonconstant positive steady states

In this section, we study the existence, stability and global bifurcation analysis
of the noncanstant steady states to model (1.3) in the 1D interval (0, l), l > 0

with prey-taxis sensitivity coefficient ξ < ξ̂. The eigenvalues of −∆ on (0, l) under
homogeneous Neumann boundary condition are σi = ( iπl )2, i = 0, 1, 2 . . . , and the

corresponding eigenfunctions are φi = cos( iπxl ), i = 0, 1, 2 . . . .

3.1. Existence of nonconstant steady states. In this part, we apply Crandall-
Rabinowitz bifurcation theory [4] with ξ as the bifurcation parameter to show that
the system (1.3) can generate a branch of nonconstant steady state solutions, which
bifurcate from the homogeneous steady state (u∗, v∗). To be precise, we consider
the corresponding strongly coupled elliptic system of model (1.3):

d1∆u+ u(1− u)− u2v

au2 + u+ g
= 0, x ∈ (0, l),

d2∆v − ξ∇(χ(v)∇u) + cv(1− v

δu
) = 0, x ∈ (0, l),

u′ = v′ = 0, x = 0, l.

(3.1)

Then (3.1) can be expressed as

F (ξ, u, v) = 0, (ξ, u, v) ∈ R×X ,

where

F (ξ, u, v) =

(
d1∆u+ u(1− u)− u2v

au2+u+g

d2∆v − ξ∇(χ(v)∇u) + cv(1− v
δu )

)
,

X = {(u, v) : u, v ∈ H2([0, l]), u′ = v′ = 0 at x = 0, l}.



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 7

It is obvious that the solution of F = 0 is the classical solution of (3.1), and F
is a continuously differentiable mapping from R× X to Y, Y = L2(0, l)× L2(0, l).
Moreover, F (ξ, u∗, v∗) = 0 for any ξ ∈ R.

If (α, u∗, v∗) is a bifurcation point of F = 0 with respect to the curve (ξ, u∗, v∗),
then the solution of F = 0 not lying on this curve in any neighbourhood of
(α, u∗, v∗). Base on (2.4) and (2.7), the derivative operator

F(u,v)(ξ, u
∗, v∗) =

(
d1∆ + f1 f2

g1 − ξχ(v∗)∆ d2∆ + g2

)
satisfies

(A+D∆)(ξi) = F(u,v)(ξi, u
∗, v∗) (3.2)

and Qi(σi, ξi) = 0 by Theorem 2.2. Therefore, F(u,v)(ξi, u
∗, v∗) has zero eigenvalue,

which implies that (ξi, u
∗, v∗) may be bifurcation points for all i ≥ 1 as follow. We

now show that the local bifurcation does occur at (ξi, u
∗, v∗) for all i ≥ 1.

Theorem 3.1. Assuming (H2) holds, and ξi 6= ξk for any k 6= i, then (ξi, u
∗, v∗) is

a bifurcation point with respect to the curve (ξ, u∗, v∗), where ξi is defined in (2.4).
That is there exists a constant δ > 0 such that for each positive integer i, system
(3.1) has non-constant positive solutions Γi(s) = (ξi(s), ui(s, x), vi(s, x)) for |s| < δ,
where ξi(s) ∈ R and (ui(s, x), vi(s, x)) ∈ X are smooth function of s. Furthermore,
Γi(s) around (ξi, u

∗, v∗) can be expressed as

(ui(s, x), vi(s, x)) = (u∗, v∗) + s(1, bi) cos
iπx

l
+ o(s),

ξi(s) = ξi + sK1 + s2K2 + . . . ,
(3.3)

where bi = d1σi−f1
f2

, K1, K2 are constants to be determined and

(ui(s, x), vi(s, x))− (u∗, v∗)− s(1, bi) cos
iπx

l

∈
{

(u, v) ∈ X :

∫ 1

0

(uûi + vv̂i)dx = 0
}
≡ Z,

where (ûi, v̂i) = (1, bi) cos iπxl .

Proof. For fixed i, we only need to prove the following conditions hold by [4, The-
orem 1.7].

(1) The partial derivatives Fξ, F(u,v) and F(ξ,u,v) exist and are continuous;
(2) dimN (F(u,v)(ξi, u

∗, v∗)) = dimY/R(F(u,v)(ξi, u
∗, v∗)) = 1,

(3) N (F(u,v)(ξi, u
∗, v∗)) = span {(ûi, v̂i)} and

F(ξ,u,v)(ξi, u
∗, v∗)(ûi, v̂i) /∈ R(F(u,v)(ξi, u

∗, v∗)).

From calculations, we obtain

Fξ(ξ, u, v) = −
(

0
∇(χ(v)∇u)

)
,

F(ξ,u,v)(ξ, u, v)(p, q) = −
(

0
∇(χ(v)∇p) +∇(χ′(v)q∇u)

) (3.4)

Clearly, Fξ, F(u,v) and F(ξ,u,v) are continuous. Let(
ϕ
ψ

)
=

∑
0≤i<∞

(
ai
bi

)
cos

iπx

l
∈ N (F(u,v)(ξi, u

∗, v∗)).



8 ZH. ZHAO, H. HU EJDE-2023/37

Then ∑
0≤i<∞

(
d1∆ + f1 f2

g1 − ξiχ(v∗)∆ d2∆ + g2

)(
ai
bi

)
cos

iπx

l
= 0, (3.5)

and therefore the coefficient matrix in (3.5) is singular, that is,

det

(
d1∆ + f1 f2

g1 − ξiχ(v∗)∆ d2∆ + g2

)
= det

(
−d1σi + f1 f2

g1 + ξiχ(v∗)σi −d2σi + g2

)
= 0, (3.6)

which implies

dimN (F(u,v)(ξi, u
∗, v∗)) = 1,

N (F(u,v)(ξi, u
∗, v∗)) = span{(ûi, v̂i)},

here (ûi, v̂i) = (1, bi) cos iπxl , with bi = d1σi−f1
f2

. The accompanying matrix of

F(u,v)(ξi, u
∗, v∗) is:

F ∗(u,v)(ξi, u
∗, v∗) =

(
d1∆ + f1 g1 − ξiχ(v∗)∆

f2 d2∆ + g2

)
.

It is easy to see that

N (F ∗(ξi,u,v)(u
∗, v∗)) = span {(û∗i , v̂∗i )} ,

where (û∗i , v̂
∗
i ) = (1, b∗i ) cos iπxl , with b∗i = − f2

g2−d2σi
. Thus

dimN
(
F ∗(u,v)(ξi, u

∗, v∗)
)

= 1,

and by the Fredholm alternative theorem (see Appendix D of [7]), we obtain

R(F(u,v)(ξi, u
∗, v∗)) =

[
N (F ∗(u,v)(ξi, u

∗, v∗))
]⊥

,

that is

dimY/R(F(u,v)(ξi, u
∗, v∗)) = dimN (F ∗(u,v)(ξi, u

∗, v∗)) = 1.

From the inner product,

(F(ξ,u,v)(ξi, u
∗, v∗)(ûi, v̂i), (û

∗
i , v̂
∗
i )) = −

∫ 1

0

σiχ(v∗)b∗i cos2 iπx

l
dx

=
χ(v∗)σilf2

2(g2 − d2σi)
> 0 .

Then F(ξ,u,v)(ξi, u
∗, v∗)(ûi, v̂i) /∈ R(F(u,v)(ξi, u

∗, v∗)). �

Remark 3.2. According to Theorem 3.1, each (ξi, u
∗, v∗) represents a bifurcation

point, that is, system (3.1) has infinite bifurcation points. And the emergence of
nontrivial steady-state bifurcation solutions is explained by the effect of prey-taxis
sensitivity coefficient ξ.

Remark 3.3. The bifurcation curve Γi(s), |s| < δ is part of a connected component
Γi of S where

S = {(ξ, u, v) ∈ R×X : F (ξ, u, v) = 0, (u, v) 6= (u∗, v∗)} .

Morever, Γi can be characterised by eigenfunction cos( iπxl ), when it is around
(ξi, u

∗, v∗).



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 9

3.2. Stability of nonconstant steady states. In this subsection, we investigate
the stability of the steady-state solution and provide criterion and explicit formulas
to determine the bifurcating direction by employing classical results from [5].

From assumption (H1), F is C3−smooth, and (ui(s, x), vi(s, x)) are C3− smooth
functions of s by Theorem 1.18 in [4], then we can expand them as follows:

ui(s, x) = u∗ + s cos
iπx

l
+ s2Φ1 + s3Φ2 + o(s3),

vi(s, x) = v∗ + sbi cos
iπx

l
+ s2Ψ1 + s3Ψ2 + o(s3),

(3.7)

where (Φk,Ψk) ∈ Z, k = 1, 2 as defined in Theorem 3.1. Moreover, from the
Taylor’s expansion, we have

χ(vi(s, x))

= χ(v∗) + sχ′(v∗)bi cos
iπx

l
+ s2(χ′(v∗)Ψ1 +

1

2
χ′′(v∗)b2i cos2 iπx

l
) + o(s3).

(3.8)

For the branch direction of Γi(s), we have the following Lemma.

Lemma 3.4. Suppose that the conditions of Theorem 3.1 are satisfied. Then for
each i ∈ N+, K1 = 0 implys the local bifurcation branch Γi(s), |s| < δ is of pitch-
fork type, where K1 is given by (3.3).

Proof. Let

M1 =
1

2

[
fuu(u∗, v∗) + 2fuv(u

∗, v∗)bi + fvv(u
∗, v∗)b2i

]
,

M2 =
1

2

[
guu(u∗, v∗) + 2guv(u

∗, v∗)bi + gvv(u
∗, v∗)b2i

]
,

To calculate the value of K1, we substitute (3.3), (3.7), and (3.8) into (3.1) and
collect the coefficients of the s2− terms, then we have

d1Φ′′1 + f1Φ1 + f2Ψ1 = −M1 cos2 iπx

l
,

d2Ψ′′1 − ξiχ(v∗)Φ′′1 + g1Φ1 + g2Ψ1 +K1(
iπ

l
)2χ(v∗) cos

iπx

l

= −ξi(
iπ

l
)2biχ(v∗) cos

2iπx

l
−M2 cos2 iπx

l
.

(3.9)

Multiplying both sides of the first equation of (3.9) by cos iπxl and then integrating
it over 0 to l, we obtain(

f1 − d1(
iπ

l
)2
)∫ 1

0

Φ1 cos
iπx

l
dx+ f2

∫ 1

0

Ψ1 cos
iπx

l
dx = 0. (3.10)

On the other hand, because (Φ1,Ψ1) ∈ Z, we have∫ 1

0

Φ1 cos
iπx

l
dx+ bi

∫ 1

0

Ψ1 cos
iπx

l
dx = 0. (3.11)

From (3.10) and (3.11), the determinant of the coefficient matrix is

det

(
f1 − d1( iπl )2 f2

1 bi

)
= −f2(1 + b2i ) 6= 0.



10 ZH. ZHAO, H. HU EJDE-2023/37

Therefore, the coefficient matrix is not singular, and we have∫ 1

0

Φ1 cos
iπx

l
dx =

∫ 1

0

Ψ1 cos
iπx

l
dx = 0, ∀i ∈ N+. (3.12)

Then K1 = 0 based on the second equation of (3.9) and (3.12). Thus, the local
bifurcation branch Γi(s), |s| < δ is of pitch-fork type. �

We next evaluate K2 to determine the bifurcation direction and the stability of
Γi(s), |s| < δ. To simplify the calculation, we introduce the following notation:

M3 = fuu(u∗, v∗) + fuv(u
∗, v∗)bi, M4 = fuv(u

∗, v∗) + fvv(u
∗, v∗)bi,

M5 = guu(u∗, v∗) + guv(u
∗, v∗)bi, M6 = guv(u

∗, v∗) + gvv(u
∗, v∗)bi,

M7 =
1

6
(fuuu(u∗, v∗) + 3fuuv(u

∗, v∗)bi + 3fuvv(u
∗, v∗)b2i + fvvv(u

∗, v∗)b3i ),

M8 =
1

6
(guuu(u∗, v∗) + 3guuv(u

∗, v∗)bi + 3guvv(u
∗, v∗)b2i + gvvv(u

∗, v∗)b3i ).

Integrating the equations of system (3.9) over 0 to l and combining K1 = 0, we
obtain ∫ 1

0

Φ1dx =
l(M2f2 −M1g2)

2(f1g2 − f2g1)
,

∫ 1

0

Ψ1dx =
l(M1g1 −M2f1)

2(f1g2 − f2g1)
. (3.13)

Multiplying the equations of system (3.9) by cos 2iπx
l and then integrating them

over 0 to l, once again combining K1 = 0 yields∫ 1

0

Φ1 cos
2iπx

l
dx =

E1

E0
,

∫ 1

0

Ψ1 cos
2iπx

l
dx =

E2

E0
, (3.14)

where

E0 = f1g2 − f2g1 −
4i2π2(ξiχ(v∗)f2 + g2d1 + f1d2)

l2
+

16i4π4d1d2

l4
,

E1 =
π2i2(ξibiχ(v∗)f2 + 2M1d2)

2l
+

(M2f2 −M1g2)l

4
,

E2 =
π2i2(2ξiχ(v∗)M1 − ξibiχ(v∗)f1 + 2M2d1)

2l
+

2π4i4ξibiχ(v∗)d1

l3

+
(M1g1 −M2f1)l

4
,

obviously, E0 is always nonzero by ξi 6= ξk for any k 6= i.
Substituting (3.7), (3.8) into (3.1), we equate s3-terms of (3.1),

d1Φ′′2 + f1Φ2 + f2Ψ2 = −M3Φ1 cos
iπx

l
−M4Ψ1 cos

iπx

l
−M7 cos3 iπx

l
,

d2Ψ′′2 − ξiχ(v∗)Φ′′2 + g1Φ2 + g2Ψ2 +K2χ(v∗)(
iπ

l
)2 cos

iπx

l

= ξiN −M5Φ1 cos
iπx

l
−M6Ψ1 cos

iπx

l
−M8 cos3 iπx

l
,

(3.15)

where

N = −χ′(v∗)( iπ
l

) sin
iπx

l
(biΦ

′
1 + Ψ′1) + χ′(v∗) cos

iπx

l
(biΦ

′′
1 − (

iπ

l
)2Ψ1)

+ χ′′(v∗)b2i (
iπ

l
)2 cos

iπx

l
(sin2 iπx

l
− 1

2
cos2 iπx

l
).



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 11

Following (Φ2,Ψ2) ∈ Z and the first equation of system (3.15), we have∫ 1

0

Φ2 cos
iπx

l
dx = −bi

∫ 1

0

Ψ2 cos
iπx

l
dx,∫ 1

0

Ψ2 cos
iπx

l
dx = − 1

(1 + b2i )f2

[M3

2
(

∫ 1

0

Φ1dx+

∫ 1

0

Φ1 cos
2iπx

l
dx)

+
M4

2
(

∫ 1

0

Ψ1dx+

∫ 1

0

Ψ1 cos
2iπx

l
dx) +

3M7l

8

]
.

By the second equation of system (3.15), we have

K2 = − l

(iπ)2χ(v∗)
(A1

∫ 1

0

Φ2 cos
iπx

l
dx+A2

∫ 1

0

Ψ2 cos
iπx

l
dx

+A3

∫ 1

0

Φ1 cos
2iπx

l
dx+A4

∫ 1

0

Ψ1 cos
2iπx

l
dx

+M5

∫ 1

0

Φ1dx+A6

∫ 1

0

Ψ1dx+A7),

(3.16)

where

A1 = 2(ξiχ(v∗)(
iπ

l
)2 + g1), A2 = 2(−d2(

iπ

l
)2 + g2),

A3 = M5 + 2ξiχ
′(v∗)bi(

iπ

l
)2, A4 = M6 − ξiχ′(v∗)(

iπ

l
)2,

A6 = M6 + ξiχ
′(v∗)(

iπ

l
)2, A7 =

3M8l

4
+ ξiχ

′′(v∗)
b2i i

2π2

8l
.

Therefore, we can evaluate K2 in terms of system parameters. For the stability of
the local bifurcating solution Γi(s), we have the following theorem.

Theorem 3.5. Assuming that the conditions of Theorem 3.1 are satisfied. Let K2

be given in (3.16), i0 be a non-negative integer such that ξi0 = ξ̂ defined by (2.4).
Then we have

(i) for each i 6= i0, Γi(s), |s| < δ is unstable;
(ii) Γi0(s), |s| < δ is stable if K2 < 0 and it is unstable if K2 > 0.

Proof. (i) Linearizing (3.1) around (ξi(s), ui(s, x), vi(s, x)), we obtain the eigenvalue
problem

F(u,v)(ξi(s), ui(s, x), vi(s, x))(u, v) = ρ(s)(u, v), (u, v) ∈ X . (3.17)

When s→ 0, (3.17) can be written as

F(u,v)(ξi, u
∗, v∗)(u, v) = ρ(0)(u, v), (3.18)

By (2.4) and (2.7), we have

Qi(σi, ξj) = f2χ(v∗)σi(ξi − ξj), 0 ≤ i, j <∞.

Indeed, if j 6= i0, then Qi0(σi0 , ξj) < 0, which implies (3.18) has a simple eigenvalue
ρ(0) > 0 by (3.2). Based on the standard eigenvalue perturbation theory in [14], for
smaller s, there exists a positive real part of the eigenvalue ρ(s) in the linearized
system (3.17). Thus the solution curve Γi(s), |s| < δ, i 6= i0 is unstable, which
means there must be i = i0 when the solution curve Γi(s), |s| < δ is stable.

The proof of (ii) is the same as [39, Theorem 4.2], we omit it here. �



12 ZH. ZHAO, H. HU EJDE-2023/37

It follows from the Remark 3.3, Lemma 3.4, and Theorem 3.5 that the bifurca-
tion curve Γi around (ξi, u

∗, v∗) is pitch-fork and the stability of the critical local
bifurcation branch Γi0(s), |s| < δ depends on K2. If Γi(s), |s| < δ is stable, then it
must be the right most one on the ξ-axis, and turns to the left at (ξi, u

∗, v∗).

3.3. Global bifurcation analysis. In this subsection, we apply the bifurcation
theory in [23] to illustrate the global structure of the bifurcation curve Γi, which
contains (ξi, u

∗, v∗), and present the following theorem to show that the local solu-
tion bifurcation can be extended to a global one.

Theorem 3.6. Assuming that the conditions of Theorem 3.1 are satisfied and σi 6=
c
d2

, the projection of bifurcation branch Γi on the ξ-axis is (−∞, ξi). Furthermore,

when ξ 6= ξi for any positive integer i and ξ < ξ̂, system (3.1) has at least one
nonconstant positive solution.

Proof. By the transformation (U, V ) = (u− u∗, v − v∗), system (3.1) becomes

d1U
′′ + (U + u∗)(1− (U + u∗))− (U + u∗)2(V + v∗)

a(U + u∗)2 + (U + u∗) + g
= 0 x ∈ (0, l),

d2V
′′ − ξ(χ(V + v∗)U ′)′ + c(V + v∗)(1− V + v∗

δ(U + u∗)
) = 0 x ∈ (0, l),

U ′(x) = V ′(x) = 0 x = 0, l.

We denote

H(ξ, U, V ) =

(
d1U

′′ + (U + u∗)(1− (U + u∗))− (U+u∗)2(V+v∗)
a(U+u∗)2+(U+u∗)+g

d2V
′′ − ξ(χ(V + v∗)U ′)′ + c(V + v∗)(1− V+v∗

δ(U+u∗) )

)
.

Let w = (p, q) ∈ X . Then we obtain

H(u,v)(ξ, U, V )w = H1w
′′ +H2(ξ, w,w′) = 0,

where χ′(V + v∗) denotes the derivative of χ(V + v∗) with respect to V ,

H1 =

(
d1 0

−ξχ(V + v∗) d2

)
,

H2(ξ, w,w′) =

(
h1

−ξ(χ(V + v∗)p′ + χ′(V + v∗)(U + u∗)′q)′ + h2

)
,

h1 =
[
1− 2(U + u∗)− (U + u∗)(U + u∗ + 2g)(V + v∗)

(a(U + u∗)2 + (U + u∗) + g)2

]
p

− (U + u∗)2

a(U + u∗)2 + (U + u∗) + g
q,

h2 =
c(V + v∗)2

δ(U + u∗)2
p+ (c− 2c(V + v∗)

δ(U + u∗)
)q.

It can be seen that, Trace(H1) > 0 and Det(H1) > 0, then the operator H1 is
elliptic by Definition 2.1 in [23], furthermore, which is strongly elliptic and satisfies
the Agmon’s condition for all θ ∈ [−π2 ,

π
2 ] by Case 2 of [23, Remark 2.5.5]. Thus the

operator H(u,v)(ξ, U, V ) is the Fredholm operator with zero index by [23, Theorem
3.3 and Remark 3.4].

Therefore, Γi must satisfy one of the following conditions:

(1) Γi on R×X is unbounded;



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 13

(2) Γi contains a point (ξ̄, u∗, v∗), where ξ̄ 6= ξi;
(3) Γi contains a point where (ξ;u, v) ∈ ∂W , where W = {(ξ;u, v) ∈ R × X :

u > 0, v > 0},
with regard to the global bifurcation result [23, Theorem 4.3]. As we know, the
positive solution of the system (3.1) bifurcates from (ξ, u∗, v∗) if and only if ξ = ξi.
Therefore, (2) is ruled out.

Evaluating F(u,v)(ξ, u, v) at (1, 0), we obtain

F(u,v)(ξ, 1, 0)(p, q) =

(
d1p
′′ − p− q

a+1+g

d2q
′′ + cq

)
.

It is easy to know that if the eigenvalues of the elliptic operator with the homo-
geneous Neumann boundary condition satisfy ( iπl )2 6= c

d2
, i = 0, 1, 2, . . . , then the

boundary (1, 0) is a non-degenerate equilibrium point. Thus (3) does not hold
either.

From the Theorem 1.1, it follows that any positive solution (u, v) of the system
(3.1) is bounded in L∞(Ω), and satisfies (u, v) ∈ C1,α(Ω) for α ∈ (0, 1). Then
the Sobolev embedding theorem implies that (u, v) possesses boundedness in the
norm of X . Thus, Γi extends to infinity in ξ, which projects on the ξ−axis as
(−∞, ξi). �

4. Numerical simulations

In this section, we present some numerical simulations to illustrate our theoreti-
cal results given in the previous sections. Meanwhile, we can find some interesting
patterns of striking structures. Based on Volume-filling mechanism, that is preda-
tors stop aggregating after their density attains certain threshold values vm and
χ(v) vanishes identically when v ≥ vm [1, 21], which satisfies the assumption (H).
We take

χ(v) :=

{
1− cos( 2πv

vm
) 0 ≤ v < vm,

0 v ≥ vm.
as the sensitivity function.

In Figure 1, the parameters of (1.3) are chosen to be d1 = 0.5, d2 = 1, a = 2
3 ,

g = 2.25, c = 3, δ = 1.5, vm = 2, which satisfy assumption (H2). By a simple
computation, we can obtain that the unique positive homogeneous steady state
solution is

(u∗, v∗) = (0.75, 1.125),

which is stable for the corresponding ODE system.
In Figure 2, parameters are the same as those for Figure 1, and we select the ini-

tial data as (u0, v0) = (u∗+0.01 cosx, v∗+0.01 cosx), which are small perturbations
of (u∗, v∗). By (2.3) and (2.4), it is easy to find that

ξ̂ = max
1≤i≤+∞

ξi = ξ5 ≈ −15.6163.

Then we observe the pattern formation for the density of predator v(x, t) when
ξ = −15.61 and −15.62 (see Figures 2(a) and 2(b)), and confirm that a steady

state bifurcation occurs at ξ̂ ≈ −15.6163. From Figures 2(a) and 2(c), we find

that the homogeneous steady state solution (u∗, v∗) is stable for any ξ > ξ̂. Turing

instability at (u∗, v∗) occurs when ξ < ξ̂, see Figures 2(b) and 2(d). These graphs
also support our stability analysis of steady state solution (u∗, v∗) in Theorem 2.2.



14 ZH. ZHAO, H. HU EJDE-2023/37

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
u

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

v

_u = 0
_v = 0
Stable Fixed Point

pd

Figure 1. Phase diagram.

Morever, the spatial pattern with mode i0 = 5 and eigenmode cos i0πxl = cos πx2 is
observed from Figures 2(b) and 2(d).

Figure 3 shows the effect of interval length on the stable wave number, when ξ
is near the critical bifurcation value. System parameters and initial data are taken
to be the same as those for Figure 2. Besides, we choose l = 13, 21 and 29, and cal-
culate the corresponding critical bifurcation values of each point to ξ7 ≈ −15.6361,
ξ11 ≈ −15.6078, and ξ15 ≈ −15.6031. We observe that the homogeneous steady
state solution (u∗, v∗) loses its stability to the stable wave mode cos 7πx

13 , cos 11πx
21

and cos 15πx
29 , respectively. Therefore, Theorem 3.5 is verified numerically in Fig-

ure 3. Moreover, as the length of interval is increases, the number of aggregates
in the spatial pattern become greater and the time for the system to arrive at a
non-constant steady state increases.

Figure 4 depicts how a small prey-taxis sensitivity coefficient ξ affects pattern
formation in system (1.3). Here the parameter values and initial values are still
chosen in the same way as in Figure 2. Choosing l = 15, and different prey-
taxis sensitivity coefficients: ξ = −20, ξ = −50 and ξ = −100, which are less
than the maximum critical bifurcation value ξ8 = −15.6235. we can get that v(x, t)
develop into spike functions when the prey-taxis sensitivity coefficient is small. This
validates the fact that a small prey-taxis sensitivity coefficient benefits heterogeneity
and aggregation of population species in (1.3).

We present a list of plots in Figure 5 to show that a variety of nonconstant steady
state solutions, spatio-temporal periodic solutions and spatio-temporal irregular so-
lutions are observed for system (1.3) as the interval length is increased when the
prey-taxis sensitivity coefficient ξ is much smaller than the critical bifurcation value.
The system parameters are taken as d1 = 0.5, d2 = 1, a = 1, g = 0.7125, c = 0.9,
δ = 0.9, vm = 2. We take initial values (u0, v0) = (u∗ + 0.05 cosx, v∗ + 0.05 cosx)
and fix ξ = −30, which is obviously far away from the critical bifurcation value.
It is clear that the solution can evolve through a series of emerging and merging.
Here merging refers to the combination of two aggregates into one via chemotactic
attraction, while emerging refers to the appearance of new aggregation peaks in the



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 15

(a) ξ = −15.61 (b) ξ = −15.62

(c) ξ = 0 (d) ξ = −17

Figure 2. Numerical simulation of v(x, t) of (1.3) over l = (0, 10).
Here d1 = 0.5, d2 = 1, a = 2/3, g = 2.25, c = 3, δ = 1.5 with the
initial data (u0, v0) = (u∗+0.01∗cos(x), v∗+0.01∗cos(x)), and ξ is
a variable coefficient. Solutions are plotted for (i) (ii): t ∈ [0, 5000],
(iii) (iv): t ∈ [0, 100].

space created by the merger. Except that, the spatio-temporal behavior of the in-
ternal spike transfer to the boundary can also be observed. These numerical results
reflect the coexistence of predator-prey species and the complexity of the spatio-
temporal structure of prey-taxis systems. For the relevant theoretical analysis we
can refer to [22].

3021
In Figure 6, system parameter are chosen to be the same as those for Figure

5, except that ξ ∈ [−30.2,−29.5], which are smaller than the critical bifurcation
value ξ8 = 6.0879. We find that the pattern formation is complex and extremely
sensitive to the prey-taxis sensitivity coefficient ξ. In particular, the dynamics are
same initially, while they vary significantly from t ≈ 200, transitioning between
regular and irregular patterning.

Finally, we find an interesting phenomenon, that is the intrinsic growth rate of
predators c can influence spatio-temporal periodicity when ξ is far away from the
critical bifurcation value. We take d1 = 0.5, d2 = 1, a = 2

3 , g = 2.25, δ = 1.5 and
l = 4, ξ = −40 in (1.3). The initial values are the same as Figure 5. It can be seen
directly from Figure 7 that there is a corresponding transition from nonconstant



16 ZH. ZHAO, H. HU EJDE-2023/37

(a) l = 13 (b) l = 21 (c) l = 29

Figure 3. Pattern formations of v(x, t) of (1.3) illustrate the ef-
fect of interval length on the stable wave number with ξ near the
critical bifurcation values. System parameters and initial data are
the same as in Figure 2.

Figure 4. Nonconstant steady state solutions of the system (1.3)
over l = (0, 15) show how a small prey-taxis sensitivity coefficient
ξ affects pattern formations. System parameters and initial data
are the same as in Figure 2.

steady state solutions to spatio-temporal periodic solutions as c increases from
2.823 to 2.824. As c increases from 2.824 to 2.925, we can observe the period of the
spatio-temporal periodic solutions becomes greater. And spatio-temporal periodic
solutions transform into nonconstant steady state solutions, when c increases from
3.005 to 3.006. We can suppose that the period becomes infinite in this process. Our
numerical simulations indicate that the intrinsic growth rate of predator c affects
the period of sustained coexistence of prey and predator species in the habitat,
which enriches the spatial-temporal patterns of the prey-taxis system.

5. Conclusions and discussion

In this article, we investigated the dynamics of predator-prey model of Leslie-
Gower type with Sigmoid functional response and prey-taxis under the homoge-
neous Neumann boundary condition. First the global existence of classical solutions
is established. Choosing the prey-taxis sensitivity coefficient ξ as the bifurcation
parameter, we establish conditions for Turing pattern formation using linear sta-
bility analysis. The existence of nonconstant positive steady state is derived from
the local bifurcation theory, and it is proved to be of pitch-fork type. Moreover, we



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 17

(a) l = 6 (b) l = 7 (c) l = 11 (d) l = 20

(e) l = 21 (f) l = 38 (g) l = 48 (h) l = 50

Figure 5. A variety of pattern formations are observed as the
interval length is increased. System parameters of (1.3) are given
as (d1, d2, a, g, c, δ) = (0.5, 1, 1, 0.7125, 0.9, 0.9) and ξ = −30.

(a) ξ = −30.2 (b) ξ = −30.1 (c) ξ = −30 (d) ξ = −29.9

(e) [ξ = −29.8 (f) ξ = −29.7 (g) ξ = −29.6 (h) ξ = −29.5

Figure 6. Transitions between regular and irregular patterning
in (1.3) when ξ ∈ [−30.2,−29.5] is far away from the critical bi-
furcation value. System parameter are chosen to be the same as
those for Figure 5.

also obtain formulas to determine the turning and stability of the local bifurcation
branch Γi(s). In addition, we show that the bifurcation curves can be extended
by applying the global bifurcation theory. Numerical simulations are provided to
illustrate and validate our theoretical results. Our findings also reveal that pattern
formation is affected by interval length, prey-taxis sensitivity coefficient ξ and the
intrinsic growth rate of predator c in (1.3).



18 ZH. ZHAO, H. HU EJDE-2023/37

c
=

2
.8

2
3

c
=

2
.8

2
4

c
=

2
.9

2
5

c
=

3
.0

0
5

c
=

3
.0

0
6

(a) u (b) v (c) phase portrait

Figure 7. Maps of solutions for the system (1.3) show how a
variable coefficient c affects pattern formations for the parameter
sets (d1, d2, a, g, δ, l, ξ) = (0.5, 1, 2/3, 2.25, 1.5, 4,−40) and initial
condition sets (u0, v0) = (u∗ + 0.05 cosx, v∗ + 0.05 cosx). Phase
portrait are plotted for t ∈ [100, 2000].



EJDE-2023/37 PATTERN FORMATION FOR A PREY-TAXIS MODEL 19

Particularly, when ξ is away from critical bifurcation value, we can observe
nonconstant steady state solutions, spatio-temporal periodic solutions and spatio-
temporal irregular solutions as the length interval increases. It is worth noting
that solutions evolve to a spatio-temporal pattern with clear temporal periodicity,
that is phase plane trajectory gives a closed orbit when c increases from 2.823 to
2.834. The period of spatio-temporal pattern becomes greater as c increases from
2.824 to 2.925. And the period tends to infinity as c tends to 3.006, that is the
spatio-temporal pattern backs to a nonconstant steady state solution. But the
mathematical behavior of the this process is not well understood yet. We speculate
that the effect of the intrinsic growth rate of predators in the pattern formation of
the predator-prey model is prevalent, which is worthwhile to explore.

Acknowledgments. The authors are grateful to the anonymous referees for their
helpful comments and valuable suggestions which have improved the presentation
of the manuscript.

References

[1] B. Ainseba, M. Bendahmane, A. Noussair; A reaction-diffusion system modeling predator-

prey with prey-taxis, Nonlinear Anal. Real World Appl., 9 (2008), no. 5, 2086–2105.
[2] N. Bellomo, A. Bellouquid, Y. Tao, M. Winkler; Toward a mathematical theory of keller-

segel models of pattern formation in biological tissues, Math. Models Methods Appl. Sci., 25

(2015), no. 9, 1663–1763.
[3] J. Cao, H. Sun, P. Hao, P. Wang; Bifurcation and Turing instability for a predator-prey

model with nonlinear reaction cross-diffusion, Appl.Math. Model., 89 (2021), 1663–1677.

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

[5] M. G. Crandall, P. H. Rabinowitz; Bifurcation, perturbation of simple eigenvalues and lin-
earized stability, Arch. Rational Mech. Anal., 52 (1973), 161–180.

[6] B. D. Dalziel, E. Thomann, J. Medlock, P. De Leenheer; Global analysis of a predator-prey

model with variable predator search rate, J. Math. Biol., 81 (2020), no. 1, 159–183.
[7] L. C. Evans; Partial differential equations, vol. 19, American Mathematical Soc., 2010.

[8] M. P. Hassell, J. H. Lawton, J. R. Beddington; Sigmoid functional responses by invertebrate

predators and parasitoids, J. Anim. Ecol., (1977), 249–262.
[9] X. He, S. Zheng; Global boundedness of solutions in a reaction-diffusion system of predator-

prey model with prey-taxis, Appl. Math. Lett., 49 (2015), 73–77.
[10] G. Hu, X. Li, Y. Wang; Pattern formation and spatiotemporal chaos in a reaction-diffusion

predator-prey system, Nonlinear Dynam., 81 (2015), no. 1-2, 265–275.

[11] J. Huang, S. Ruan, J. Song; Bifurcations in a predator-prey system of Leslie type with

generalized Holling type III functional response, J. Differential Equations, 257 (2014), no. 6,
1721–1752.

[12] J. Jorné; Negative ionic cross diffusion coefficients in electrolytic solutions, J. Theor. Biol.,
55 (1975), no. 2, 529–532.

[13] P. Kareiva, G. Odell; Swarms of predators exhibit “preytaxis” if individual predators use

area-restricted search, Am. Nat., 130 (1987), no. 2, 233–270.
[14] T. Kato; Perturbation theory for linear operators, Classics in Mathematics, Springer-Verlag,

Berlin, 1995, Reprint of the 1980 edition.

[15] R. Kubo, A. Ugajin, M. Ono; Molecular phylogenetic analysis of mermithid nematodes
(mermithida: Mermithidae) discovered from japanese bumblebee (hymenoptera: Bombinae)

and behavioral observation of an infected bumblebee, Appl. Entomol. Zool., 51 (2016), no. 4,

549–554.
[16] J. M. Lee, T. Hillen, M. A. Lewis; Pattern formation in prey-taxis systems, J. Biol. Dyn., 3

(2009), no. 6, 551–573.

[17] P. H. Leslie; Some further notes on the use of matrices in population mathematics,
Biometrika, 35 (1948), 213–245.



20 ZH. ZHAO, H. HU EJDE-2023/37

[18] P. H. Leslie, J. C. Gower; The properties of a stochastic model for the predator-prey type of

interaction between two species, Biometrika 47 (1960), 219–234.

[19] J. D. Murray; Mathematical biology. I, third ed., Interdisciplinary Applied Mathematics,
vol. 17, Springer-Verlag, New York, 2002, An introduction.

[20] W. Ni, M. Tang; Turing patterns in the Lengyel-Epstein system for the CIMA reaction,

Trans. Amer. Math. Soc., 357 (2005), no. 10, 3953–3969.
[21] K. J. Painter, T. Hillen; Volume-filling and quorum-sensing in models for chemosensitive

movement, Can. Appl. Math. Q., 10 (2002), no. 4, 501–543.

[22] K. J. Painter, T. Hillen; Spatio-temporal chaos in a chemotaxis model, Physica D, 240 (2011),
no. 4-5, 363–375.

[23] J. Shi, X. Wang; On global bifurcation for quasilinear elliptic systems on bounded domains,

J. Differential Equations, 246 (2009), no. 7, 2788–2812.
[24] Y. Song, X. Tang; Stability, steady-state bifurcations, and Turing patterns in a predator-prey

model with herd behavior and prey-taxis, Stud. Appl. Math., 139 (2017), no. 3, 371–404.
[25] X. Tang, Y. Song; Bifurcation analysis and turing instability in a diffusive predator-prey

model with herd behavior and hyperbolic mortality, Chaos Solitons Fractals, 81 (2015),

no. part A, 303–314.
[26] X. Tang, Y. Song; Stability, hopf bifurcations and spatial patterns in a delayed diffusive

predator-prey model with herd behavior, Appl. Math. Comput. 254 (2015), 375–391.

[27] Y. Tao; Global existence of classical solutions to a predator-prey model with nonlinear prey-
taxis, Nonlinear Anal. Real World Appl., 11 (2010), no. 3, 2056–2064.

[28] C. Tian, Z. Ling, Z. Lin; Turing pattern formation in a predator-prey-mutualist system,

Nonlinear Anal. Real World Appl., 12 (2011), no. 6, 3224–3237.
[29] A. M. Turing; The chemical basis of morphogenesis, Philos. Trans. Roy. Soc. London Ser.

B, 237 (1952), no. 641, 37–72.

[30] R. K. Upadhyay, A. Patra, B. Dubey, N. K. Thakur; A predator-prey interaction model with
self-and cross-diffusion in aquatic systems, J. Biol. Syst. 22 (2014), no. 04, 691–712.

[31] J. Wang, J. Wei, J. Shi; Global bifurcation analysis and pattern formation in homogeneous
diffusive predator-prey systems, J. Differential Equations, 260 (2016), no. 4, 3495–3523.

[32] M. Wang; Stability and Hopf bifurcation for a prey-predator model with prey-stage structure

and diffusion, Math. Biosci., 212 (2008), no. 2, 149–160.
[33] Q. Wang, Y. Song, L. Shao; Nonconstant positive steady states and pattern formation of 1D

prey-taxis systems, J. Nonlinear Sci., 27 (2017), no. 1, 71–97.

[34] X. Wang, W. Wang, G. Zhang; Global bifurcation of solutions for a predator-prey model with
prey-taxis, Math. Methods Appl. Sci., 38 (2015), no. 3, 431–443.

[35] H. F. Weinberger; Invariant sets for weakly coupled parabolic and elliptic systems, Rend.

Mat., (6) 8 (1975), 295–310.
[36] S. Wu, J. Shi, B. Wu; Global existence of solutions and uniform persistence of a diffusive

predator-prey model with prey-taxis, J. Differential Equations, 260 (2016), no. 7, 5847–5874.

[37] S. Wu, J. Wang, J. Shi; Dynamics and pattern formation of a diffusive predator-prey model
with predator-taxis, Math. Models Methods Appl. Sci. 28 (2018), no. 11, 2275–2312.

[38] L. Zhang, S. Fu; Global bifurcation for a holling-tanner predator-prey model with prey-taxis,
Nonlinear Anal. Real World Appl., 47 (2019), 460–472.

[39] T. Zhang, X. Liu, X. Meng, T. Zhang; Spatio-temporal dynamics near the steady state of a

planktonic system, Comput. Math. Appl., 75 (2018), no. 12, 4490–4504.

Zhihong Zhao
School of Mathematics and Physics, University of Science & Technology Beijing, Bei-

jing 100083, China

Email address: zzh@ustb.edu.cn

Huanqin Hu

School of Mathematics and Physics, University of Science & Technology Beijing, Bei-
jing 100083, China

Email address: bjkjdxhh@163.com


	1. Introduction
	1.1. Global solutions

	2. Stability analysis of homogeneous steady state
	3. Nonconstant positive steady states
	3.1. Existence of nonconstant steady states
	3.2. Stability of nonconstant steady states
	3.3. Global bifurcation analysis

	4. Numerical simulations
	5. Conclusions and discussion
	Acknowledgments

	References