EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 16, No. 2, 2023, 899-918
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Dynamical Behaviors of Four-Dimensional
Prey-Predator Model

A. G. Farhan1,∗, Nihad Sh. Khalaf2, Talat J. Aldhlki 1

1 Department of Mathematics, College of Basic Education, Mustansiriyah University,
Baghdad, Iraq
2 Department of Mathematics, College of Education for Women, Tikrit University, Tikrit,
Iraq

Abstract. In this paper, we study the dynamic behavior of a four-dimensional prey-predator
suggested model of four species. The four species are two prey and two predator species, each of
them grows logistically. The two prey live in diverse habitats and have the ability of group defense.
In the mentioned model, one predator feeds on the two prey, the top predator feeds on other three
species. The existents and, the boundedness of the positive solution, the existence and the local
stability of all possible equilibrium points, of the model are investigated. The model has seven
equilibrium points at most, four of them always exist and the others exist under certain conditions.
Three equilibrium points are not stable while the others are locally asymptotically stable, under
given conditions. For the coexistence point, a basin of attraction for it has been found. The
steady-state bifurcation relative to the mortality rate of the predators in the neighborhood of
three of the equilibrium points and the Hopf-bifurcation relative to the growth rate of the prey in
the neighborhood of two of the equilibrium points has been found. Finally, two numerical example
has been given to support the theoretical results.

2020 Mathematics Subject Classifications: 65P30, 34C23, 92D25

Key Words and Phrases: Basin of attraction, Bifurcation, Equilibrium point, Group defence,
stability

1. Introduction

Predators feed in a habitat that is relatively rich in food for some time, which means
to them that there are large numbers of prey, or that such prey is easy to catch. When
food is scarce, or prey is minor in the habitat, the predators look for another habitat with
enough food to live in for another period of time, this phenomenon is called switching, see
[8, 9, 13]. Group defense means that all prey animals attack predators collectively and do

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v16i2.4716

Email addresses: a.g.farhan.edbs@uomustansiriyah.edu.iq (A. G. Farhan),
nihad.shreef16@tu.edu.iq (N.Sh. Khalaf),
dr.talatjassim.edbs@uomustansiriyah.edu.iq (T.G. Aldhlki)

https://www.ejpam.com 899 © 2023 EJPAM All rights reserved.



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 900

not give them an opportunity to attack and prey on them. In 1920, Volterra introduced
a mathematical model that represents the interaction between prey and predator [1, 10].
Later many modified Volterra models were introduced by researchers interested in this
field. Various models, including prey with one predator, two predators with one prey, two
prey with two predators, and food chains of three or more species have been studied in
[2–4, 15, 16]. In the wild, lions are at the top of the food chain because of their strength
and ability to kill all animals and prey on them, but it is possible for a herd or a small
group of hyenas to attack a young lion or an old lion outside its kingdom. In general,
predators such as lions, hyenas, wolves, and others attack prey such as zebra, wild buffalo,
and other prey with the aim of killing and devouring them. Prey by nature has the ability
to live in groups that move from one habitat to another in search of food. Some prey
animals, such as buffaloes and zebras, have the ability to collectively defend and attack
predators in masse [3] and [8]. In [3], two prey-predatory models were studied, one of
the two models was with group prey defense, while the other model was without group
prey defense. The prey was supposed to live in two different habitats. Both predatory
species tend to move from one habitat to another in search of food. The first model has
been expanded into the following mathematical model that proposed and studied, in [5]
and [7] with switching index n = 1,2, respectively. The expanded model deals with, two
prey and two predators, all species grow logistically, the prey lives in diverse habitats with
prey group defense, one predator tends to switch habitats and feeds on the prey, while the
other predator feeds on only one prey.

ẋ1 = x1

(
g1

(
1− x1

k1

)
− α1x

n
2y1

xn1 + xn2
− βy2

)
,

ẋ2 = x2

(
g2

(
1− x2

k2

)
− α2x

n
1y1

xn1 + xn2

)
,

ẏ1 = y1

(
−µ1 +

δ1x1x
n
2

xn1 + xn2
+

δ2x
n
1x2

xn1 + xn2

)
,

ẏ2 = y2(−µ2 + γx1)

(1)

The aim of this work is to study the local qualitative behaviors in the neighborhood of the
equilibrium points of a proposed prey-predator mathematical model for four species with
an incomplete food chain. The four species, consist of two prey and two predators, and
each of them grows logistically. The two prey live in two different habitats and have the
ability to group defense. One predator or the top feeds on the other predator, in addition
to the two prey indicated in the model. We show that all positive solutions are bounded
under some conditions. In the third section, it was found that the model has, at most
seven equilibrium points. Local stability and local bifurcation of the equilibrium points
have been analyzed in section4 and section5 respectively. Two sets of parameters form
two systems of 4-dimensional differential equations, each of them is an example of the
model (2) presented in the sixth section of this paper to simulate the theoretical results
obtained in the previous sections. Finally a brave conclusion of this paper is given in the
last section.



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 901

2. THE MATHEMATICAL MODEL

In this work, the following model was proposed, which is a modification of model (1)

ẋ1 = x1

(
g1

(
1− x1

k1

)
− α1x2y1

x1 + x2
− α2x2y2

x1 + x2

)
,

ẋ2 = x2

(
g2

(
1− x2

k2

)
− β1x1y1

x1 + x2
− β2x1y2

x1 + x2

)
,

ẏ1 = y1

(
−µ1 +

ϵ1x1x2
x1 + x2

− ρ2y2

)
,

ẏ2 = y2

(
−µ2 +

ϵ2x1x2
x1 + x2

+ ρ1y1

)
(2)

where xi denote the density of the prey i = 1, 2; y1 denote the density of the predator
that fed on the two prey; y2 denote the density of the predator that fed on the other three
species (top predator). gi, i = 1, 2 is the growth rate of xi, i = 1, 2; ki, i = 1, 2, is the
carrying environmental capacity to xi, i = 1, 2; αi, i = 1, 2, is the rate of predation by the
predatoryi, i = 1, 2 with the prey x1;βi, i = 1, 2 is the rate of predation by the predator
yi, i = 1, 2 on prey x2; µi, i = 1, 2, is the mortality rate of predatorsyi, i = 1, 2; εi, i = 1, 2,
is the corresponding conversion rates to µi, i = 1, 2; ρ1, is the rate of change of y2 due
to the presence of y1,; and ρ2, is the rate of change of y1 due to the presence of y2. The
switching behavior of predators y1 and y2 is shown by, the functions α1x2y1(x1 + x2)

−1,
β1x1y1(x1 + x2)

−1 and, the functions α2x2y2(x1 + x2)
−1, β2x1y2(x1 + x2)

−1 respectively.
It is easy to show that all functions of the model (2) and their partial derivatives are
continuous, so that these functions are Lipschitizion functions on

R+4 =
{
(x1, x2, y1, y2) ∈ R4 : xi(0) ≥ 0, yi(0) ≥ 0, i = 1, 2, x1 + x2 > 0

}
.

Hence, the existence and the uniqueness of the solution of model (2) is guaranteed. Now,
we will show that the trajectories of all the positive solutions of the system (2) whose
initial conditions lie within the following region D are bounded

D :=
{
(x1, x2, y1, y2) ∈ R4, 0 < xi < ki, yi > 0, i = 1, 2

}
. (3)

Theorem 1. If εi ≤ αi + βi, i = 1, 2, ρ1 ≤ ρ2, then the trajectories of all the positive
solutions of the system (2) whose initial conditions belong to D are bounded.

Proof. The following real valued function:

z(t) = x1(t) + x2(t) + y1(t) + y2(t),

is a positive definitive on D. Therefore, we get

ż =
2∑

i=1

(
xigi

(
1− xi

ki

)
− µiyi −

(ρ2 − ρ1)

2
y
1
y2 −

x1x2yi
x1 + x2

(αi + βi−εi)

)



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 902

If 0 < ρ ≤ max {µ1, µ2}, then we obtain:

ż +ρu =
2∑

i=1

(
xigi

(
1− xi

ki
+

ρ

gi

)
+(ρ−µi)yi−

(ρ2 − ρ1)

2
y1y2 −

x1x2yi
x1 + x2

(αi + βi − εi)

)
.

It is obvious that

ż + ρz <

2∑
i=1

xigi

(
1− xi

ki
+

ρ

gi

)

=
2∑

i=1

xi
ki
(kigi − gixi + kiρ) <

2∑
i=1

xi
ki
(kigi + kiρ) <

2∑
i=1

ki(gi + ρ).

Put α =:
∑2

i=1 ki (gi + ρ). So it is clear that: 0 ≤ z (t) ≤ α
ρ + z (0) e−ρt, and when

t → ∞, 0 ≤ z(t) ≤ α

ρ
. (4)

So that, from (4), the trajectories of all the positive solutions of the system (2) with initial
conditions lie within the region D defined above and satisfy εi ≤ αi+βi, i = 1, 2, ρ1 ≤ ρ2
are bounded. and with this, we have completed the proof.

3. EXISTENCE OF EQUILIBRIUM POINTS

The system has four equilibrium points always exist, regardless of the values of the
system parameters plus three others, whose existence depends on the change in the values
of the system parameters. The equilibrium point are given as follows:

(i) The equilibrium points P0 = (0, 0, 0, 0), P1 = (k1, 0, 0, 0), P2 = (0, k2, 0, 0) and P3 =
(k1, k2, 0, 0) always exist.

(ii) The equilibrium point P4 = (x̃1, x̃2, ỹ1, 0), where:

x̃1 =
µ1(1 + x̃)

ϵ1
, x̃2 =

x̃1
x̌
, ỹ1 = g1

(
1− x̃1

k1

)
(1 + x̃)

α1
(5)

and x̃ is a positive root of the equation (7) below, provided that

x̃ <
ε1k1 − µ1

µ1
(6)

A1x̃
3 + A2x̃

2 + A3x̃+ A4 = 0, (7)

such that

A1 = g1β1k2µ1 > 0, A2 = g1β1k2(µ1 − k1ε1) < 0,



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 903

A3 = g2α1k1(k2ε1 − µ1), A4 = −g2α1k1µ1 < 0.

The number of sign changes between A1,A2,A3, and A4 is one if A3 < 0, and three
if A3 > 0. According to Descartes’ rule of sign, equation (7) has at most one positive
root if A3 < 0, and at most three positive root if A3 > 0. Therefore, the existence
of P4 depends on the existence of the positive root x̃ achieves the two equations (6)
and (7).

(iii) The equilibrium point P5 = (x̌1, x̌2, 0, y̌2) , where:

x̌1 =
µ2 (1 + x̌)

ε2
, x̌2 =

x̌1
x̌
, y̌2 = g1

(
1− µ2 (1 + x̌)

ε2k1

)
(1 + x̌)

α2
(8)

and x̌ is a positive root of equation (10) below, provided that:

x̌ <
ε2k1 − µ2

µ2
, (9)

B1x̌
3 + B2x̌

2 + B3x̌+ B4 = 0, (10)

such that:

B1 = g1β2k2µ2, B2 = g1β2k2 (µ2 − k1ε2) < 0,

B3 = g2α2k1 (k2ε2 − µ2) , B4 = −g2α2k1µ2.

The number of sign changes between B1,B2,B3, and B4 is one if B3 < 0, and, three
if B3 > 0. Based on Descartes’ rule of sign., eq.(10) has at most one positive root
if B3 < 0, and at most three positive root if B3 > 0. Therefore, the existence of P5

depends on the existence of the positive root x̌ achieves the two equations (9) and
(10).

(iv) The interior equilibrium point or the coexistence point P 6 = (x1, x2, y1, y2), where

x1 =
k1k2(1 + h) ((α1β2 − α2β1)hµ2 + ρ1(1 + h)(α2g2h− β2g1))

ρ1(1 + h)2 (k1 α2g2h
2 − k2β2g1) + k1k2ε2h2 (α1β2 − α2β1)

,

x2 = hx1, y1 =
(1 + h)µ2 − hε2

ρ1 (1 + h)
, y2 =

hε1 − (1 + h)µ1

ρ2 (1 + h)
,

(11)

and h = x2x
−1
1 , is a positive real root of the following equation:

Ah4 + Bh3 + Ch2 + Dh+ E = 0, (12)

where:

A = k1g2(α2µ1ρ1 − α1µ2ρ2 + ρ1ρ2g1),

B = 2k1g2(α2µ1ρ1 − α1µ2ρ2) + 3ρ1ρ2g1g2k1 − k2ρ1ρ2g1g2



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 904

+ k1k2g2(α1ϵ2ρ2 − α2ϵ1ρ1),

C = k1g2(α2µ1ρ1 − α1µ2ρ2) + k2g1(β1µ2ρ2 − β2µ1ρ1) − k1k2g2(α2ϵ1ρ1 − α1ϵ2ρ2)

+ k1k2g1(β2ϵ1ρ1 − β1ϵ2ρ2) + k1k2(α1β2 − α2β1)(µ1ϵ2 − µ2ϵ1) + 3ρ1ρ2g1g2(k1 − k2),

D = 2k2g1(β1µ2ρ2 − β2µ1ρ1)− 3k2ρ1ρ2g1g2 + ρ1ρ2g1g2k1 + k1k2g1(β2ϵ1ρ1 − β1ϵ2ρ2),

E = k2g1(β1µ2ρ2 − β2µ1ρ1 − ρ1ρ2g2).

If the signs of the coefficients A,B,C,D, and E are positive, a positive root cannot
be obtained, and this means that there is no interior equilibrium point. Otherwise,
there is a possibility of finding one positive root, two positive roots, three positive
roots, or four at most, according to Descartes’ rule of signs.

4. LOCAL STABILITY OF EQUILIBRIUM POINTS

In this section, the local stability of the equilibrium points of the system is studied, by
using the following Jacobian matrix J (x1, x2, y1, y2) of the system (2) at each equilibrium
point [6], [11], [12].

J(x1, x2, y1, y2) =


∂ẋ1
∂x1

∂ẋ1
∂x2

∂ẋ1
∂y1

∂ẋ1
∂y2

∂ẋ2
∂x1

∂ẋ2
∂x2

∂ẋ2
∂y1

∂ẋ2
∂y2

∂ẏ1
∂x1

∂ẏ1
∂x2

∂ẏ1
∂y1

∂ẏ1
∂y2

∂ẏ2
∂x1

∂ẏ2
∂x2

∂ẏ2
∂y1

∂ẏ2
∂y2


where,

∂ẋ1
∂x1

= g1

(
1− 2x1

k1

)
−

(α1y1+α2y2)x
2
2

(x1 + x2)
2 ,

∂ẋ1
∂x2

= −
(α1y1+α2y2)x

2
1

(x1 + x2)
2 ,

∂ẋ1
∂y1

= −α1x1x2
x1 + x2

,
∂ẋ1
∂y2

= −α2x1x2
x1 + x2

,

∂ẋ2
∂x1

= −
(β1y1+β2y2)x

2
2

(x1 + x2)
2 ,

∂ẋ2
∂x2

= g2

(
1− 2x2

k2

)
−

(β1y1+β2y2)x
2
1

(x1 + x2)
2 ,

∂ẋ2
∂y1

= − β1x1x2
x1 + x2

,
∂ẋ2
∂y2

= − β2x1x2
x1 + x2

,

∂ẏ1
∂x1

=
ε1x

2
2y1

(x1 + x2)
2 ,

∂ẏ1
∂x2

=
ε1x

2
1y1

(x1 + x2)
2 ,

∂ẏ1
∂y1

= −µ1 +
ε1x1x2
x1 + x2

− ρ2y2,
∂ẏ1
∂y2

= −ρ2y1,

∂ẏ2
∂x1

=
ε2x

2
2y2

(x1 + x2)
2 ,

∂ẏ2
∂x2

=
ε2x

2
1y2

(x1 + x2)
2 ,

∂ẏ2
∂y1

= ρ1y2,
∂ẏ2
∂y2

= −µ2 +
ε1x1x2
x1 + x2

+ ρ1y1,



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 905

In Theorems 2, 3 and 4, we will prove that the points P0, P1, and P2 are unstable point
while the remaining points are locally asymptoticaly stable points if they meet certain
conditions.

Theorem 2. Consider the system (2), then the equilibrium points

(i) P0 = (0, 0, 0, 0), P1 = (k1, 0, 0, 0) and P2 = (0, k2, 0, 0) are unstable equilibrium
points.

(ii) P3 = (k1, k2, 0, 0) locally asymptotically stable point if k1k2εi < (k1+k2) µi, i = 1, 2.

Proof.

(i) Suppose that P0 = (0, 0, 0, 0) is locally asymptotically stable. So that all the trajec-
tories (x1, x2, y1, y2) of the system (2) converge to (0, 0, 0, 0) as t → ∞. Then since
x1 > 0, we have that d

dt (lnx1) → g1, as t → ∞.

It is possible to find a small ball with center P0 and radius g1, such that inside it
we have d

dt (lnx1 ) ≥ g1
2 . if (x1, x2, y1, y2) converges to (0, 0, 0, 0), when t converges

to ∞, then there exists t0 > 0, such that;

x1(t0) > 0, x1(t) ≥ x1(t0)exp
(
g1(t−t0)

2

)
→ ∞, as t → ∞. So x1 → ∞. Similarly, if

x2(0) > 0, x1 → ∞, there is no a trajectory to the system (2) converges to (0, 0, 0, 0).
Hence P0 = (0, 0, 0, 0) is unstable point.

The Jacobian matrices of the system (2) at the equilibrium points P1 and P2 are

J(P1) =


−g1 0 0 0
0 g2 0 0
0 0 −µ1 0
0 0 0 −µ2

 and J(P2) =


g1 0 0 0
0 −g2 0 0
0 0 −µ1 0
0 0 0 −µ2

 ,

respectively. The eigenvalues of the two matrices J (P1) and J (P2) are:

λ1 = −g1 < 0, λ2 = g2 > 0, λ3 = −µ1 < 0, λ4 = −µ2 < 0

and
λ1 = g1 > 0, λ2 = −g2 < 0, λ3 = −µ1 < 0, λ4 = −µ2 < 0,

respectively. It is clear that, we have one positive eigenvalue in both cases. Therefore,
the two points P1 and P2 are saddle points and this means that they are not stable.

(ii) The Jacobian matrix of the system (2) at the equilibrium point P3 = (k1, k2, 0, 0) is

J(P3) =


−g1 0 −α1k1k2(k1 + k2)

−1 −α2k1k2(k1 + k2)
−1

0 −g2 −β2k1k2(k1 + k2)
−1 −β2k1k2(k1 + k2)

−1

0 0 ε1k1k2(k1 + k2)
−1 − µ1 0

0 0 0 ε2k1k2(k1 + k2)
−1 − µ2





A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 906

The solutions of the following characteristic equations |J (P3)− λI4×4| = 0, where
I4×4 is the identity matrix, are:

λ1 = −g1 < 0, λ2 = −g2 < 0, λ3 =
ε1k1k2
k1 + k2

− µ1, λ4 =
ε2k1k2
k1 + k2

− µ2.

So that, P3 = (k1, k2, 0, 0) is locally asymptotically stable if the following conditions
are met

k1k2εi < (k1 + k2)µi, i = 1, 2. (13)

Theorem 3. Consider the system (2), then

(i) The equilibrium point P4 = (x̃1, x̃2, ỹ1, 0) is locally asymptotically stable if

Γ̃ε2x̃1 + ρ1ỹ1 − µ2 < 0, trM̃ < 0,detM̃ < 0, and
3∑

i=1

trM̃M̃ii < detM̃. (14)

where M̃ =

Γ̃2α1x̃ỹ1 − k−1
1 x̃1g1 −Γ̃2α1x̃

2ỹ1 −Γ̃α1x̃1
−Γ̃2β1ỹ1 Γ̃2β1x̃ỹ1 − k−1

2 x̃2g2 −Γ̃β2x̃1
Γ̃ε1ỹ1 Γ̃ε1x̃

2ỹ1 Γ̃ε1x̃1 − µ1

 and Γ̃ = 1
1+x̃ .

(ii) The equilibrium point P5 = (x̌1, x̌2, 0, y̌2) is locally asymptotically stable point if

Γ̌ε1x̌1 − ρ2y̌2 − µ1 < 0, trM̌ < 0, detM̌ < 0, and

3∑
i=1

trM̌M̌ii < detM̌. (15)

where M̌ =

Γ̌2α1x̌y̌2 − k−1
1 x̌1g1 −Γ̌2α1x̌

2y̌2 − Γ̌α1x̌1
−Γ̌2β2y̌2 Γ̌2β2x̌y̌2 − k−1

2 x̌2g2 Γ̌β2x̌1
Γ̌2ε2y̌2 Γ̌2ε2x̌

2y̌2 Γ̌ε1x̌1 − µ2

 and Γ̌ = 1
1+x̌

Proof.

(i) The Jacobian matrix of the system (2) at equilibrium points P4 = (x̃1, x̃2, ỹ1, 0) , is

J (P4) =


Γ̃2α1x̃ỹ1 − k−1

1 x̃1g1 −Γ̃2α1x̃
2ỹ1 −Γ̃α1x̃1 −Γ̃α2x̃1

−Γ̃2β1ỹ1 Γ̃2β1x̃ỹ1 − k−1
2 x̃2g2 −Γ̃β1x̃1 −Γ̃β2x̃1

Γ̃2ε1ỹ1 Γ̃2ε1x̃
2ỹ1 Γ̃ε1x̃1 − µ1 −ρ2ỹ1

0 0 0 Γ̃ε2x̃1 + ρ1ỹ1 − µ2


such that Γ̃ = 1

1+x̃ .

The characteristic equation for the matrix J (P4) is as follows:

(
λ− Γ̃ε2x̃1 − ρ1ỹ1 + µ2

)(
λ3 − trM̃ λ2 +

3∑
i=1

M̃iiλ−detM̃

)
= 0,



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 907

where, M̃ =

Γ̃2α1x̃ỹ1 − k−1
1 x̃1g1 −Γ̃2α1x̃

2ỹ1 −Γ̃α1x̃1
−Γ̃2β1ỹ1 Γ̃2β1x̃ỹ1 − k−1

2 x̃2g2 −Γ̃β2x̃1
Γ̃2ε1ỹ1 Γ̃2ε1x̃

2ỹ1 Γ̃ε1x̃1 − µ1


and M̃ii, i = 1, 2, 3 is the diagonal minor’s of the matrix M̃ .

According to the criteria of Routh-Hurwitz, P4 = (x̃1, x̃2, ỹ1, 0) , is a locally asymp-
totically stable point provided that (14) are satisfied.

(ii) The Jacobi’s matrix of the system (2) at P5 = (x̌1, x̌2, 0, y̌2) is

J (P5) =


Γ̌2α1x̌y̌2 − k−1

1 x̌1g1 −Γ̌2α2x̌
2y̌2 −Γ̌α1x̌1 −Γ̌α2x̌1

−Γ̌2β2y2 Γ̌2β2x̌y̌2 − k−1
2 x̌2g2 −Γ̌β1x̌1 −Γ̌β2x̌1

0 0 Γ̌ε1x̌1 − ρ2y̌2 − µ1 0

Γ̌2ε2y̌2 Γ̌2ε2x̌
2y̌2 ρ1y̌2 Γ̌ε1x̌1 − µ2


such that Γ̌ = 1

1+x̌ .

Simple calculations yield that the characteristic equation of the matrixJ (P5) is:

(
Γ̌ε1x̌1 − ρ2y2 − µ1 − λ

)(
λ3 − trM̌ λ2 +

3∑
i=1

M̌ii λ−detM̌

)
= 0,

where

M̌ =

Γ̌2α1x̌y̌2 − k−1
1 x̌1g1 −Γ̌2α1x̌

2y̌2 −Γ̌α1x̌1
−Γ̌2β2y̌2 Γ̌2β2x̌y̌2 − k−1

2 x̌2g2 Γ̌β2x̌1
Γ̌2ε2y̌2 Γ̌2ε2x̌

2y̌2 Γ̌ε1x̌1 − µ2


and M̌ii, i = 1, 2, 3 is the diagonal minor’s matrix M̌.

According to the criteria of Routh-Hurwitz, any of the equilibrium points P5 =
(x̌1, x̌2, 0, y̌2) , is locally asymptotically provided that (15) are satisfied. And with
this, we have completed the proof.

Theorem 4. Consider the system (2), then P6 = (x1, x2, y1, y2) is a locally asymptotically
stable point if 

trM > 0,

trM ∆−
4∑

i=1

Mii > 0,

4∑
i=1

Mii

(
trM ∆−

4∑
i=1

Mii

)
−detM trM2

> 0,

(16)

where, for i = 1, 2, 3, 4, Mii is the diagonal minor’s of the matrix M = J(P6)



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 908

Proof. The Jacobian matrix of the system (2) at equilibrium points P6 = (x1, x2, y1, y2)
is as follows

J (P6) =


Q1 − g1x1

k1
−Q1

h α1E α2E

−h2Q2 hQ2 − g2x2

k2
β1E β2E

h2F 1 F1 0 −ρ2y1
h2F2 F2 ρ1y2 0

 =: M,

such that, for i = 1, 2, Qi = gi

(
1− xi

ki

)
1

1+h , E = −x2
1+h , and Fi =

εiyi
(1+h)2

, further, the

characteristic equation of the matrix M is

λ4 − trM λ3 +∆ λ2 −
4∑

i=1

Mii λ+ detM = 0,

where, Mii, i = 1, 2, 3, 4 is the diagonal minor’s of the matrix M, and

∆ =
g1x1g2x2
k1k2

+ ρ1ρ2y1y2 −
2∑

i=1

(
EFi

(
h2αi + βi

)
+ xi

giQ(i−(−1)i)

ki

)
.

According to the criteria of Routh-Hurwitz, P6 = (x1, x2, y1, y2) is locally asymptotically
stable point if the (16) is satisfied. And with this, we have completed the proof.

In the following theorem, we give conditions through which the stability region of the
coexistence point be asymptotically stable.

Theorem 5. Assume that P6 = (x1, x2, y1, y2) is locally asymptotically stable point, and
xi ≥ ki, i = 1, 2.Then, the set B, which is defined below represents an attraction basin for
P6.

B = {(x1, x2, y1, y2) : xi ≥ xi, y1 ≤ y1, y2 = y2} .

Proof. The function

V (x1, x2, y1, y2) =

2∑
i=1

(
xi − xi − xiln

xi
xi

+ yi − yi − yiln
yi
yi

)
is positive definite .

V̇ (x1, x2, y1, y2) =
2∑

i=1

(
ẋi

(
1− k1

xi

)
+ ẏi

(
1− yi

yi

))

=

2∑
i=1

(xi − xi)Gi (x1, x2) + (yi − yi)G3 (x1, x2),

such that

G1 (x1, x2) =

(
g1

(
1− x1

k1

)
− α1x2y1

x1 + x2
− α2x2y2

x1 + x2

)
,



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 909

G2 (x1, x2) =

(
g2

(
1− x2

k2

)
− α2x1y1

x1 + x2

)
,

G3 (x1, x2) =

(
−µ1 +

ε1x1x2
x1 + x2

− ρ2y2

)
.

It is obvious that (x1 − x1) > 0, (x2 − 2) > 0,G1 (x1, x2) < 0 and G2 (x1, x2) < 0, so that,

2∑
i=1

(xi − xi)Gi (x1, x2) < 0.

The function G3 (x1, x2) is increasing with respect to xi, i = 1, 2, because of the positivity
of it derivative with respect to xi for xi ≥ xi, i = 1, 2 as shown below:

∂G3 (x1, x2)

∂x1
=

ε1x
2
2

(x1 + x2)
2 > 0,

∂G3 (x1, x2)

∂x2
=

ε1x
2
1

(x1 + x2)
2 > 0, for xi ≥ xi, i = 1, 2.

Now, since G3 (x1, x2) = 0, then it is obtained that G3 (x1, x2) > 0, for xi ≥ xi, i = 1, 2.
Hence, (y1 − y1)G3 (x1, x2) < 0. So that,

2∑
i=1

(xi − xi)Gi (x1, x2) + (yi − yi)G3 (x1, x2) < 0,

V̇ (x1, x2, y1, y2) < 0, ∀ (x1, x2, y1, y2) ∈ B − {(x1, x2, y1, y2)} ,

and V̇ (x1, x2, y1, y2) = 0. So that, B is an attraction basin for P6 = (x1, x2, y1, y2) and
with this, we have completed the proof.

5. LOCAL BIFURCATION

This section is dedicated to the study of the local bifurcation of the three equilibrium
points. Local bifurcation occurs when a small change in the parameter value changes
the behavior of the equilibrium. If the corresponding real part of one eigenvalue of the
equilibrium passes through zero, a steady-state bifurcation occurs, and in the case that the
eigenvalues is not zeros but is purely imaginary, then this is a Hopf-bifurcation [12, 14].

(I) If the equilibrium point P (ϑ) is locally asymptotically stable and for ϑ = ϑ∗ we
have one eigenvalue λ (ϑ∗) = 0, the bifurcation that occurs is a steady state bifurcation
[6]. The following theorems show that steady-state bifurcation occurs at points P3 =
(k1, k2.0, 0) , P4 = (x̃1, x̃2, ỹ1, 0), and P5 = (x̌1, x̌2, 0, y̌2) for the parameters µi, i = 1 or 2,
µ1, and µ2, respectively.

Theorem 6. Suppose that, the equilibrium point P3 = (k1, k2.0, 0) of the system (2) is
locally asymptotically stable point, then

(i) A steady state bifurcation occurs as µ1 pass through µ∗
1 =

k1k2ε1
k1+k2

,

(ii) A steady state bifurcation occurs as µ1 pass through µ∗
2 =

k1k2ε2
k1+k2

.



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 910

Proof. (i) From Theorem 2 and the assumption or the theorem that we have that the
eigenvalues of the characteristic equation of J(P3) are

λ1 = −g1, λ2 = −g2, λ3 =
ε1k1k2
k1 + k2

− µ1 < 0 , λ4 =
ε2k1k2
k1 + k2

− µ2 < 0.

It obvious that P3 = (k1, k2.0, 0) do not depends on the parameter µ1, which means that
P3 do not change with the change of the value of µ1. Hence, if µ1 > µ∗

1 then P3 of the
system (2) is locally asymptotically stable point (nod point) and if µ1 < µ∗

1 then P3 of the
system (2) is unstable point (saddle point).
(ii) The proof is the same as the proof of the first part of this theorem, and with this, we
have completed the proof.

Theorem 7. (i) Suppose that, the equilibrium point P4 = (x̌1, x̌2, 0, y̌2) of the system
(2) exists and is locally asymptotically stable point then a steady state bifurcation
occurs as µ2 pass through µ∗

2 =
ε2x̌1
1+x̌ + ρ1y̌1.

(ii) Suppose that, the equilibrium point P5 of the system (2) exist and is locally asymp-
totically stable point then a steady state bifurcation occurs as µ1 pass through µ∗

1 =
ε1x̃1
1+x̃ − ρ2ỹ2.

Proof. (i) From Theorem 3, the characteristic equation of J(P4) is

(
λ− Γ̃ε2x̃1 − ρ1ỹ1 + µ2

)(
λ3 − trM̃ λ2 +

3∑
i=1

M̃ii λ−detM̃

)
= 0.

As it is clear from the equation that second algebraic expression
(
λ3 − trM̃λ2 +

∑3
i=1 M̃iiλ−detM̃

)
of the left side the characteristic equation does not depend on the parameter µ2. The equa-
tions (5), (6) and (8), show that the existence of P4 does not depends on µ2. Hence, if
µ2 > µ∗

2 the P3 of the system (2) is locally asymptotically stable point (nod point) and if
µ2 < µ∗

2 the P4 of the system (2) is unstable point (saddle point).
(ii) The proof is the same as the proof of first part of this theorem.

And with this, we have completed the proof.

(II) The local birth or death of a periodic solution from equilibrium as a parameter ϑ
passes through a critical value ϑ∗ is called Hopf bifurcation.
For the parameter g = g1 = g2, let P4 be a locally asymptotically stable equilibrium point
and

trM̃

3∑
i=1

M̃ii = detM̃, (17)

g∗ =
S2∓

√
S2
2−4S1S3

2S1
, is the positive root of the equation (17) where,

S1 = ÃÑ22, S2 = detÑ − Ã
(
Ñ11 + Ñ22

)
− B̃Ñ33, S3 = B̃

(
Ñ11 + Ñ22

)
,



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 911

à =
2∑

i=1

(
k−1
i x̃i

)
− Γ̃2x̃ỹ (α1 + β1) , B̃ = (µ1 − Γ̃ε1x̃1),

Ñ =

k−1
1 x̃1 − Γ̃2α1x̃ỹ Γ̃2α1x̃

2ỹ Γ̃α1x̃1
Γ̃2β1ỹ k−1

2 x̃2 − Γ̃2β1x̃ỹ Γ̃β2x̃1
−Γ̃ε1x̃

−1
2 ỹ −Γ̃ε1x̃x̃

−1
2 ỹ µ1 − Γ̃ε1x̃1

 ,

ỹ =

(
1− x̃1

k1

)
(1 + x̃)

α1
.

From (17), eigenvalues of (P4(g
∗)) are

λ1,2 (g
∗) = ±i

√√√√ 3∑
i=1

M̌ii, λ3 (g
∗) = trM̃(g), λ4(g

∗) = Γ̃ε2x̃1 + ρ1ỹ1 − µ2.

Theorem 8. Suppose that the equilibrium point P4 = (x̃1, x̃2, ỹ1, 0), of the system (2)
exists and is locally asymptotically stable for the parameter g = g1 = g2, and g∗ be the pos-
itive root of the equation (17), then a simple Hopf bifurcation occurs as g passes through g∗

provided that(
M̃11 + M̃22

)(
µ1 − Γ̃ε1x̃1

)
g−1 − M̃33

(
k−1
1 x̃1−Γ̃2α1x̃ỹ + k−1

2 x̃2 − Γ̃2β1x̃ỹ
)
̸=0.

Proof. Since ∂x̃1
∂g = ∂x̃2

∂g = ∂ỹ2
∂g = 0, ∂ỹ1∂g = ỹ1

g = ỹ, then the point P4 = (x̃1, x̃2, ỹ1, 0),
depends smoothly on the parameter g. If there is a simple pair of complex eigenvalues
λ1,2 (g) = u± iv of the Jacobian matrix (J (P4(g))) at the equilibrium point P4 (g) , such
that, it becomes a purely imaginary at g = g∗, while all other eigenvalues are real and
negative; and du(g∗)

dg ̸= 0, then at g∗we have a simple Hopf bifurcation, [5]. It is easy to
deduce,

du (g∗)

dg
=

(
d(detM̃)

dg −
∑3

i=1 M̃ii
dtrM̃
dg − trM̃ db

dg

)
2
∑3

i=1 M̃ii + 2
(
trM̃

)2 .

Because 2
∑3

i=1 M̃ii + 2
(
trM̃

)2
is positive, so it is sufficient to prove that

W :=
d(detM̃)

dg
− trM̃

d
∑3

i=1 M̃ii

dg
−

3∑
i=1

M̃ii
d trM̃
dg

̸= 0.

Simple calculation show that

W =
d(detM̃)

dg
−

3∑
i=1

M̃ii
d trM̃
dg

− trM̃
d (
∑3

i=1 M̃ii)

dg

= (2g) det Ñ − trM̃
(
Ñ11 + Ñ22 + 2gÑ33

)
−

3∑
i=1

M̃ii

(
k−1
1 x̃1 − Γ̃2α1x̃ỹ + k−1

2 x̃2 − Γ̃2β1x̃ỹ
)



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 912

= trM̃
(
Ñ11 + Ñ22

)
−
(
gÑ11 + gÑ22+g2Ñ33

)(
k−1
1 x̃1−Γ̃2α1x̃ỹ + k−1

2 x̃2 − Γ̃2β1x̃ỹ
)

=
(
M̃11 + M̃22

)
(µ1 − Γ̃ε1x̃1)g

−1 − M̃33

(
k−1
1 x̃1 − Γ̃2α1x̃ỹ + k−1

2 x̃2 − Γ̃2β1x̃ỹ.
)

And with this, we have completed the proof.

Now consider the following equation

trM̌
3∑

i=1

M̌ii = detM̌ (18)

g∗∗ =
r2∓

√
r22−4r1r3
2r1

is the positive root of the (18) where

r1 = ǍŇ22, r2 = detŇ − Ǎ
(
Ň11 + Ň22

)
− B̌Ň33, r3 = B̌

(
Ň11 + Ň22

)
,

Ǎ =

2∑
i=1

(
k−1
i x̃i

)
− Γ̃2x̃ỹ (α1 + β1) , B̌ = (µ1 − Γ̃ε1x̃1),

Ň =

Γ̌2α1x̌y̌ − k−1
1 −Γ̌2α1x̌

2y̌ −Γ̌α1x̌1
−Γ̌2β2y̌ Γ̌2β2x̌y̌ − k−1

2 x̌2 Γ̌β2x̌1
Γ̌2ε2y̌ Γ̌2ε2x̌

2y̌ Γ̌ε2x̌1 − µ2

 and Γ̌ =
1

1 + x̌

Theorem 9. Suppose that the equilibrium point P5 = (x̌1, x̌2, 0, y̌2) of the system (2) exists
and is locally asymptotically stable for the parameter g = g1 = g2, and g∗∗ be the positive
root of the equation (17), then a simple Hopf bifurcation occurs as g passes through g∗∗ pro-
vided that

(
M̌11 + M̌22

)
(µ2 − Γ̌ε1x̌1)g

−1 − M̌33

(
k−1
1 x̌1 − Γ̌2α1x̌y̌ + k−1

2 x̌2 − Γ̌2β2x̌y̌.
)
̸=0,

where y̌ = g−1y̌2, and g∗∗ is the positive root of the (17) given above.

Proof. The proof is the same as the proof of Theorem 8.

6. NUMERICAL EXAMPLE

In this section, two numerical examples are given to confirm the obtained theoretical
results in the above sections. Consider is the set of parameters in the following table

Table 1: model parameters1

I gi ki αi βi µi ϵi ρi
1 0.5 1.9 0.2 0.49 0.12 0.19 0.1

2 0.51 1.8 0.21 0.5 0.19 0.2 0.12

The point P3 = (k1, k1, 0, 0) is unstable, but with µ1 > µ∗
1 = 0.1756 is locally asymp-

totically stable as shown in Figure 1, this means P3 has a steady state bifurcation when µ1

pass throw µ∗
1 = 0.1756 . The point P4 = (1.6927, 1.0075, 0.731, 0), is unstable, see Figure



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 913

2, whether the equilibrium point P5 = (1.6934, 0.7095, 0, 0.87686), exist when µ2 = 0.1,
and is locally asymptotically stable as shown in Figure 3. The point P6 = (1.7098, 1.1715,
0.5096, 0.1007), exist and is locally asymptotically stable as shown in Figure 4.
Now consider the table bellow

Table 2: model parameters 2

I gi ki αi βi µi ϵi ρi
1 0.3 2 0.1 0.4 0.12 0.1 0.1

2 0.4 3.2 0.2 0.5 0.13 0.2 0.11

The point P3 = (k1, k1, 0, 0) is unstable, but with µ1 > µ∗
1 = 0.12308 and µ2 > µ∗

2 =
0.24615 is locally asymptotically stable as shown in Figure 6 This means P3 has a steady
bifurcation when µ1 pass throw µ∗

1 = 0.12308, and µ2 pass throw µ∗
2 = 0.24615.

The point P4 = (1.9630, 3.0872, 0.0907, 0), is unstable, see Figure 7, whether the equilib-
rium point P5 = (1.4549, 1.1749, 0, 0.9151), exist and is locally asymptotically stable as
shown in Figure 8. The point P6 = does not exist.

Figure 1: The trajectory of system (2), according to the parameters listed in Table 1 and µ1 > µ∗
1 = 0.1756,

starts in (1.7,1.6,0.05,0.08), which located close to P3 tend to P3.



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 914

Figure 2: The trajectory of system (2), according to the parameters listed in Table 1, starts in the initial point
(1.6, , 1, 0.7, 0.01), which located close to P4 and it is moving away from P4 and approaching P6.

Figure 3: The trajectory of system (2), according to the parameters listed in table 1, starts in the initial point
(1.6, 0.7, 0.02, 0.8), which located close to P5, is approaching P5.



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 915

Figure 4: The trajectory of system (2), according to the parameters listed in Table 1, starts in the initial point
(1.4, 1.3, 0.3, 0.05), which located close to P6 is approaching P6.

Figure 5: The trajectory of system (2), according to the parameters listed in Table 1, starts in the initial point
(1.4, 1.3, 0.3, 0.05), which located close to P6 is approaching P6.



A. G. Farhan, N. Sh. Khalaf, T. J. Aldhlki / Eur. J. Pure Appl. Math, 16 (2) (2023), 899-918 916

Figure 6: The trajectory of system (2), according to the parameters listed in Table 2 and µ1 > µ∗
1 = 0.12308,

and µ2 > µ∗
2 = 0.24615 starts in the initial point (1.8,2.9,0.3,0.1) that located close to P3 tends to P3.

Figure 7: The trajectory of system (2), according to the parameters listed in Table 2 starts in the initial point
(1.7, 2.9, 0.07, 0.03) that located close to P4 tends to P5



REFERENCES 917

Figure 8: The trajectory of system (2), according to the parameters listed in Table 2 starts in the initial point
(1.3, 1, 0.3, 0.8) is located close to P5 tends to P5

7. CONCLUSION

In this paper, the dynamical behavior of a four-dimensional prey-predator model for
four species was presented. The model has seven equilibrium points, four of them always
exist, whatever the values of the model parameters. The positive trajectory of the model
are bounded under some certain conditions. The system has three, unstable equilibrium
points and four equilibrium points that are asymptotically locally stable under some con-
ditions. The appearance of steady-state and, Hopf bifurcation near the equilibrium points
has been discussed. In the numerical example represented by Table 1, all the possible
equilibrium points exist, and in the example represented by Table 2, the coexistence point
disappeared. In both examples, the steady-state bifurcation was present in the neighbor-
hood of P3; P4 is not stable while P5 is locally asymptotically stable.

References

[1] A H Adamu. Mathematical analysis of predator-prey model with two preys and one
predator. International Journal of Engineering and Applied Sciences, 5:17–23, 2018.

[2] S Antipova. Mathematical modeling of competition two ideologies with internal con-
flicts. Bulletin of Sib GUTI, 4:27–42, 2022.

[3] R Bhattacharyya and B Mukhopadhyay. Spatial dynamics of nonlinear prey-predator
models with pry migration and predator switching. Ecological complexity, 3:160–169,
2006.



REFERENCES 918

[4] A Chatterjee and S Pal. Switching effects driven by predation on diffusive preda-
tor prey system. Applications and Applied Mathematics: An International journal,
16:682–704, 2021.

[5] A G Farhan. Study of a nonlinear model for two prey and two predators species.
Journal of the College of Basic Education, Mustansiriyah University, 100:73–86, 2018.

[6] A G Farhan. Lotka-volterra model with prey-predators food chain. Iraqi Journal of
Science, 0:56–63, 2020.

[7] A G Farhan. On the mathematical model of two-prey and two-predator species. Iraqi
Journal of Science, 3:608–619, 2020.

[8] A G Farhan and M A Rasheed. On dynamics of nonlinear prey–predator model
involving group defense. J. of College of Education, Mustansiriyah University, 5:249–
264, 2014.

[9] P. Fisher. Sur legilibre de faunas: Interactions des moules, des poupreset des
cripedaed. C.R. Soc. Biologeogr, 92:47–48, 1934.

[10] V Hadžiabdić, MMehuljić, and J Bektešević. Lotka-volterra model with two predators
and their prey. TEM Journal, 6:132–136, 2017.

[11] C S Holling. Principles of insect predation. Annual Review of Entomology, 6:163–182,
1961.

[12] S L Jianhua and W H Ni. Steady-state bifurcation and hopf bifurcation for a diffusive
leslie–gower predator–prey model. Computers and Mathematics with Applications,
70:3043–3056, 2015.

[13] J H Lawton, J R Beddington, and R Bonser. Switching in invertebrate predators.
Springer US, Boston, MA, 1974.

[14] J L Moiola and G Chen. Hopf bifurcation analysis: a frequency domain approach.
World Scientific, 1996.

[15] S Rahman, S Islam, and S Sarwardi. Effects of prey refuge with holling type iv func-
tional response dependent prey predator model. International Journal of Modelling
and Simulation, 0:1–19, 2023.

[16] M Tansky. Switching effects in prey–predator system. J. Theor. Biol, 70:263–271,
1978.