316 

© 2025 The Author(s). Published by College of Education for Pure Science (Ibn Al-Haitham), 

University of Baghdad. This is an open-access article distributed under the terms of the Creative Commons 

Attribution 4.0 International License 

 
 

 
The Local Bifurcation of Food Chain Prey-Predator Model with Crowley-

Martin-Type of Functional Response 

 

Asmaa Abdulhussein Aziz1*   and Azhar Abbas Majeed 2  
1,2Department of Mathematics, College of Science, University of Baghdad, Baghdad, Iraq. 

*Corresponding Author. 

 

Received: 4 June 2023                     Accepted: 26 February 2024                  Published: 20 April 2025 

doi.org/10.30526/38.2.3553 

 
 Abstract                                                         

     The conditions under which the local bifurcation including (saddle-node, transcritical, and 

pitchfork) of all stable points in a model of a food chain occur are examined in this paper. The 

Growly-Martin model and the emotion of fear have both been important factors in the 

development of a functional reaction in a food chain model, It has been shown that a 

transcritical bifurcation and a pitchfork bifurcation can be discovered near to each of the sites 

𝐴1, A2, 𝐴3𝑎𝑛𝑑 𝐴4, and that a saddle-node bifurcation can be found close to the point where 

positive equilibrium is located. In conclusion, a numerical simulation was run in order to 

illustrate how the proposed model may display bifurcation in its behavior. This was done in 

order to show how the impact of parameters on the dynamics of the proposed model.  

Keywords: Transcritical pitchfork, Sotomayor's theorem, Local bifurcation, Global 

bifurcation. 

 

1. Introduction  

     Bifurcation theory is the scientific study of how a moving system's structure changes over 

time. Changes that can be qualitative or visual in integral family curves, vector fields, 

differential equation solutions. Bifurcation happens when a small, smooth change in the values 

of a system's parameters (bifurcation parameters) causes a sudden shift in its "specific" or 

"topological" behavior. This term is used a lot when mathematicians study systems that change 

over time (1-5). There are two major categories of bifurcations: local and global. Local 

bifurcations like the saddle node, transcritical pitchfork, period-doubling flip, Hopf, and 

Neimark secondary Hopf bifurcation, which can all be studied purely by analyzing how their 

local stability properties change, happen as parameters cross critical thresholds. Global 

bifurcation happens when stability and bigger invariant sets, like periodic orbits, collide. 

Contrary to local bifurcations, this alters the structure of the pathways in phase space in a non-

local manner. Topological shifts are widespread because they can have an impact on areas over 

incredibly large distances; for instance, the homoclinic (6-8).Pamuk and Cay (9) looked at the 

steady state with respect to Hopf bifurcation of a feedback diffusion method, which is used to 

control how vascular cells and inhibitors talk to each other. Mukherjee and Maji (10) 

https://creativecommons.org/licenses/by/4.0/
https://creativecommons.org/licenses/by/4.0/
https://orcid.org/0009-0007-7885-1840
mailto:asmaaaziz1998@gmail.com
https://orcid.org/0000-0002-9652-6132
mailto:azhar.majeed@sc.uobaghdad.edu.iq


IHJPAS. 2025, 38(2) 

317 
 

established the criteria for local bifurcations for all equilibrium points and Hopf for positive 

equilibrium points in a prey-predator model with prey haven. While Majeed and Ali (11) 

presented a model of bifurcation with refuge that included a predator-structured stage food 

chain. In recent years, many researchers have studied the importance of bifurcation. Majeed 

and Alabacy (12) found the conditions for local bifurcation at all equilibrium points a prey-

predator-refuge-harvesting model. Since then, lots of researchers have considered the fear 

effect of predators in the dynamical study on prey-predator models (13-14). Many other 

academics have studied the local bifurcation in recent years, including (15-19). 

Finally, in this work, the occurrence of the local bifurcation (LB) of the proposed system have 

been discussed. 

 

2. Materials and Methods 

2.1. Model formulation  

Consider the following system that given in (20). 
 𝑑𝑅1

𝑑𝑡
=

𝑟1𝑅1

1+𝐾1𝑅2
− 𝑚1𝑅1

2 −
𝛽1𝑅1𝑅2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
= 𝑓1(𝑅1 . 𝑅2 . 𝑅3).  

𝑑𝑅2

𝑑𝑡
=

𝑟2𝑅2

1+𝐾2𝑅3
− 𝑚2𝑅2

2 +
𝑙1𝑅1𝑅2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

𝛽2𝑅2𝑅3

(1+𝑎3𝑅2)(1+𝑎4𝑅3)
− 𝑑1𝑅2 = 𝑓2 (𝑅1 . 𝑅2 . 𝑅3)      (1)                                                                                                            

 
𝑑𝑅3

𝑑𝑡
=

𝑙2𝑅2𝑅3

(1+𝑎3𝑅2)(1+𝑎4𝑅3)
− 𝑑2𝑅3 = 𝑓3 (𝑅1 . 𝑅2 . 𝑅3).  

Following is a table describing the positive parameters of system (1). 

Table 1. The system's parameters (1): 

Parameters Description Parameters 

THE DENSITY OF PREY . INTERMEDIATE PREDATOR AND TOP − PREDATOR AT TIME T.  RESPECTIVELY RII. = 1. 2. 3 

THE INTRINSIC GROWTH  RATE OFTHE PREY AND THE INTERMEDIATE PREDATOR RESPECTIVELY RII. = 1. 2 

THE RATE OF INTERNAL COMPETITION FOR THE PREY AND THE INTERMEDIATE PREDATOR RESPECTIVELY MI.I = 1. 2 

THE RATE OF FEAR OF PREY AND INTERMEDIATE PREDATOR RESPECTIVELY KI. I = 1. 2 

ATTACK RATE FOR PREY AND INTERMEDIATE PREDATOR RESPECTIVELY ΒI. I = 1. 2 

FOOD TRANSFER RATE FOR INTERMEDIATE PREDATOR AND TOP PREDATOR RESPECTIVELY LI. I = 1. 2 

HANDLING TIME OF PREY AND INTERMEDIATE PREDATOR RESPECTIVELY AI. I = 1. 3 

MAGNITUDE OF DISTRUBANCE AMONG INTERMEDIATE PREDATOR AND TOP PREDATOR  RESPECTIVELY AI. I = 2. 4 

THE NATURAL DEATH RATE OF INTERMEDIATE PREDATOR AND TOP PREDATOR RESPECTIVELY DI. I = 1. 2 

 

2.2. Local bifurcation analysis 

In this subsection, the local bifurcation of the model (1) has been examined, with a particular 

emphasis on the changes that occur around each equilibrium point as a result of changes in the 

parameter values governing the dynamic behavior. With the assistance of Sotomayor's 

theorem, our mission is to develop higher order conditions with the intention of ensuring the 

appearance of the most typical of the area's local bifurcations. 

Now, according to Jacobean matrix  𝐽(𝑅1. 𝑅2. 𝑅3) of the system (1) which is given in (16) as 

follows: 

𝐽 =  [𝑢𝑖𝑗]3×3
                                                                                                                                            (2)     

where: 

u11 =
r1

1+K1R2
− 2m1R1 −

β1R2

(1+a1R1)2(1+a2R2)
 ، 𝑢12 − (

r1R1K1

(1+K1R2)2
+

β1R1

(1+a1R1)(1+a2R2)2
) ، 𝑢13 = 0، 

u21 =
𝑙1R2

(1+a1R1)2(1+a2R2)
 ، 𝑢22 =

r2

1+K2R3
− 2m2R2 +

𝑙1R1

(1+a1R1)(1+a2R2)2
−

β2R3

(1+a3R2)2(1+a4R3)
d1، 

u23 = −
r2R2K2

(1 + K2R3)2
−

β2R2

(1 + a3R2)(1 + a4R3)2
، u31 = 0،u32 =

𝑙2R3

(1 + a3R2)2(1 + a4R3)
 ، 



IHJPAS. 2025, 38(2) 

318 
 

u33 =
𝑙2R2

(1+a3R2)(1+a4R3)2
− d2. 

For any non- zero vector T =(𝑡1. 𝑡2. 𝑡3)
𝑇: 

𝐷2𝐹𝜇(𝑋. 𝜇)(𝑇. 𝑇) = [𝑌𝑖1]3×1.                                                                                                                 (3) 

𝑌11 = 2((
𝛽1𝑅2𝑎1

(1+𝑎1𝑅1)3(1+𝑎2𝑅2)
− 𝑚1) 𝑡1

2 − (
𝑟1𝐾1

(1+𝐾1𝑅2)2
+

𝛽1

(1+𝑎1𝑅1)2(1+𝑎2𝑅2)2
) 𝑡1𝑡2 + (

𝑟1𝐾1
2

(1+𝐾1𝑅2)3
+

𝛽1𝑎2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)3
)  𝑡2

2𝑅1 , 

𝑌21 = −2(
𝑙1R2𝑎1

(1+a1R1)3(1+a2R2)
𝑡1

2 −
𝑙1

(1+a1R1)2(1+a2R2)2
𝑡1𝑡2 − (

r2K2
2R2

(1+K2R3)3
+

β2 a4R2

(1+a3R2)(1+a4R3)3
)𝑡3

2 + (𝑚2 +
𝑙1R1𝑎2

(1+a1R1)(1+a2R2)3
−

β2R3𝑎3

(1+a3R2)3(1+a4R3)
) 𝑡2

2 + (
r2K2

(1+K2R3)2
+

β2

(1+a3R2)2(1+a4R3)2
) 𝑡2𝑡3) , 

𝑌31 = −2(
𝑙2𝑅3𝑎3

(1+a3R2)3(1+a4R3)
𝑡2

2 −
𝑙2

(1+a3R2)2(1+a4R3)2
𝑡2𝑡3 +

𝑙2𝑅2𝑎4

(1+a3R2)(1+a4R3)3
𝑡3

2) , 

and   𝐷3𝐹𝜇(𝑋. 𝜇)(𝑇. 𝑇. 𝑇) = [𝑊𝑖1]3×1 ,                                                                                      (4) 

𝑊11 = 6 [(
𝑟1𝐾1

2

(1+𝐾1𝑅2)3
+

𝛽1𝑎2

(1+𝑎1𝑅1)2(1+𝑎2𝑅2)3
) 𝑡1𝑡2

2 + (
𝑡2

1+𝑎2𝑅2
−

R2𝑎1𝑡1

1+𝑎1𝑅1
)

𝛽1𝑎1𝑡1
2

(1+𝑎1𝑅1)3(1+𝑎2𝑅2)
−

(
𝑟1𝐾1

3𝑅1

(1+𝐾1𝑅2)4
+

𝛽1𝑎2
2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)4
) 𝑡2

3] . 

𝑊21 = 6[(
r2K2

2

(1+K2R3)3
+

β2 a4

(1+a3R2)2(1+a4R3)3
) 𝑡2𝑡3

2 + (
R2𝑎1𝑡1

1+a1R1
−

𝑡2

1+𝑎2𝑅2
)

𝑙1𝑎1𝑡1
2

(1+a1R1)3(1+𝑎2𝑅2)
+

(
𝑙1R1𝑎2

2

(1+a1R1)(1+a2R2)4
−

β2R3𝑎3
2

(1+a3R2)4(1+a4R3)
) 𝑡2

3 + (
β2 a4𝑡3

(1+a3R2)3(1+a4R3)2
−

𝑙1𝑎2𝑡1

(1+a1R1)2(1+a2R2)3
)𝑡1

2 −

(
r2K2

3𝑅2

(1+K2R3)4
+

𝛽1𝑎2
2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)4
) 𝑡3

3], 

𝑊31 = 6 [(
𝑎3𝑅3𝑡2

(1+a3R2)
−

𝑡3

(1+a4R3)
)

𝑙2𝑎3𝑡2
2

(1+a3R2)3(1+a4R3)
+ (

𝑎4𝑅2𝑡3

(1+a4R3)
−

𝑡2

(1+a3R2)
)

𝑙2𝑎4𝑡3
2

(1+a3R2)(1+a4R3)3
] , 

where 𝑋 = (𝑅1. 𝑅2. 𝑅3) and (𝜇 ) any parameter. 

Theorem (1): System (1) with the parameter value 𝑑`1 = 𝑑1 = 𝑟2 +
𝑙1𝑟1

𝑚1+𝑎1𝑟1
 , has a 

transcritical and pitchfork bifurcation at 𝐴1 = (
𝑟1

𝑚1
. 0. 0), if the next conditions hold: 

𝑍`1 ≠ 𝑍`2,                                                                                                                                   (5) 

𝑍`3 ≠ 𝑍`4.                                                                                                                               (6) 

Where: 

𝑍`1 = (
𝑙1𝑁

(1+a1R1)2(1+a2R2)2
+

β2R3𝑎3

(1+a3R2)3(1+a4R3)
),  

𝑍`2 = (
𝑙1R2𝑎1𝑁2

(1+a1R1)3(1+a2R2)
+ 𝑚2 +

𝑙1R1𝑎2

(1+a1R1)(1+a2R2)3
), 

𝑍`3 = (
𝑙1R2𝑎1

2𝑁3

(1+a1R1)4(1+a2R2)
+

𝑙1R1𝑎2
2

(1+a1R1)(1+a2R2)4
),  

𝑍`4 = (
𝑙1𝑎1𝑁2

(1+a1R1)3(1+a2R2)2
+

𝑙1𝑎2𝑁

(1+a1R1)2(1+a2R2)3
+

β2R3𝑎3
2

(1+a3R2)4(1+a4R3)
). 

Proof: By using the Jacobian matrix in equation  (1.7) in (20) 𝐽`1 = 𝐽1(𝐴1. 𝑑`1) = [𝑞`𝑖𝑗]3×3
, 

where 𝑞`𝑖𝑗 = 𝑞𝑖𝑗 , except 𝑞`22 = 0 .Then, the characterizing equation for𝐽`1 has an eigenvalue 

of zero, which is λ1R2
 at  𝑑`1 = 𝑑1 ,  

Now, let 𝑇[1] = (𝑡1
[1]

. 𝑡2
[1]

. 𝑡3
[1]

)𝑇 be the eigenvector associated with an eigenvalue λ1R2
= 0. 

Thus,(𝐽`1
 −  λ1R2

𝐼)𝑇[1] = 0.  this gives: 

𝑡1
[1]

= 𝑁𝑡2
[1]

. 𝑡3
[1]

= 0.𝑤ℎ𝑒𝑟𝑒 𝑁 = −(
𝑟1𝐾1

𝑚1
+

𝛽1

𝑚1+𝑎1𝑟1
)and 𝑡2

[1]
any real number that is not zero. 



IHJPAS. 2025, 38(2) 

319 
 

Let 𝑀[1] = (𝑚1
[1]

. 𝑚2
[1]

. 𝑚3
[1]

)𝑇 be the eigenvector associated with an eigenvalue λ1R2
= 0 of 

the matrix 𝐽`1
 𝑇.Then.  (𝐽`1

 𝑇 −  λ1R2 
𝐼) 𝑀[1] = 0 . 

By solving this equation for. 𝑀[1] = (0.𝑚2
[1]

. 0)
𝑇

.where 𝑚2
[1]

 any real number that is not zero. 

Now consider this: 
𝜕𝑓

𝜕𝑑1
= 𝑓𝑑1

(𝑋. 𝑑1) = (
𝜕𝑓1

𝜕𝑑1
.
𝜕𝑓2

𝜕𝑑1
.
𝜕𝑓3

𝜕𝑑1
)
𝑇

= (0.−𝑅2. 0)𝑇 . 

So, 𝑓𝑑1
(𝐴1. 𝑑`1) = (0. 0. 0)𝑇 and hence (𝑀[1])

𝑇
𝑓𝑑1

(𝐴1. 𝑑`1) = 0 

Using Sotomayor's theorem, it is impossible to satisfy the saddle-node bifurcation condition. 

The first condition for transcritical bifurcation is therefore satisfied. Now 

𝐷𝑓𝑑1
(𝑋. 𝑑1) = [

0     0     0
0 − 1   0
0     0     0

], 

where, 𝐷𝑓𝑑1
(𝑋. 𝑑1) represents the derivative of 𝑓𝑑1

(𝑋. 𝑑1) with respect to 𝑋 = (𝑅1. 𝑅2. 𝑅3)
𝑇 . 

Furthermore, it is observed that: 

𝐷𝑓𝑑1
(𝐴1. 𝑑1)𝑇

[1] = [
0     0     0
0 − 1   0
0     0     0

] [
𝑁𝑡2

[1]

𝑡2
[1]

0

] = [

0

−𝑡2
[1]

0

] 

(𝑀[1])
𝑇
[𝐷𝑓𝑑1

(𝐴1. 𝑑`1)𝑇
[1]] = (0.𝑚2

[1]
. 0)

𝑇

(0.−𝑡2
[1]

. 0) = −𝑡2
[1]

𝑚2
[1]

≠ 0. 

By substituting 𝑇[1] in equation (3) we get: 

𝐷2𝐹𝜇(𝐴1. 𝑑`1)(𝑇
[1]. 𝑇[1]) = (𝑌`𝑖𝑗).                                                                                                   (7) 

𝑌`11 = 2𝑡2
2[(

R2𝑎1𝑁2

(1+𝑎1𝑅1)2
−

𝑁

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+

R1𝑎2

(1+a2R2)2
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (

𝐾1𝑅1

(1+𝐾1𝑅2)
−

𝑁)
𝑟1𝐾1

(1+𝐾1𝑅2)2
− 𝑚1𝑁

2], 

𝑌`21 = 2𝑡2
2[(

𝑁

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

R1𝑎2

(1+a2R2)2
−

𝑙1R2𝑎1𝑁2

(1+a1R1)2
)

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
− 𝑚2 +

β2R3𝑎3

(1+a3R2)3(1+a4R3)
], 

𝑌`31= 
−2𝑡2

2𝑙2𝑅3𝑎3

(1+a3R2)3(1+a4R3)
 , 

Hence, it was obtained 

(𝑀[1])
𝑇
[𝐷2𝐹𝜇(𝐴1. 𝑑`1)(𝑇

[1]. 𝑇[1])] = 2 (𝑡2
[1]

)
2

𝑚2
[1]

(Z`1 − 𝑍`2) ≠ 0. 

This indicates that the system (1) exhibits a transcritical bifurcation at 𝐴1 with a parameter 

 𝑑`1 = 𝑑1, If condition (5) not satisfied then. 

By substituting 𝑇[1] in equation (4) we get: 

𝐷3𝐹𝜇(𝐴2. 𝑑`1)(𝑇
[1]. 𝑇[1]. 𝑇[1]) = (𝑊`𝑖𝑗),                                                                                         (8) 

𝑊`11 = 6𝑡2
3[(𝑁 −

𝐾1𝑅1

(1+𝐾1𝑅2)
)

𝑟1𝐾1
2

(1+𝐾1𝑅2)3
+ (

𝑎1𝑁2

(1+𝑎1𝑅1)2(1+𝑎1𝑅1)
+

𝑎2𝑁

(1+𝑎1𝑅1)(1+𝑎2𝑅2)2
−

R2𝑎1
2𝑁3

(1+𝑎1𝑅1)3
−

𝑎2
2𝑅1

(1+𝑎2𝑅2)2
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
], 

𝑊`21 = 6𝑡2
3[(

R2𝑎1
2𝑁3

(1+a1R1)3
−

𝑎1𝑁2

(1+a1R1)2(1+𝑎2𝑅2)
−

𝑎2𝑁

(1+𝑎1𝑅1)(1+a2R2)2
+

R1𝑎2
2

(1+a2R2)3
)

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

β2R3𝑎3
2

(1+a3R2)4(1+a4R3)
], 

𝑊`31=
6𝑡2 3𝑙2𝑅3𝑎3

2

(1+a3R2)4(1+a4R3)
, 

(𝑀[1])
𝑇
[𝐷3𝐹𝜇(𝐴2. 𝑑`1)(𝑇

[1]. 𝑇[1]. 𝑇[1])] = 6 (𝑡2
[1]

)
3

𝑚2
[1](𝑍`3 − 𝑍`4) ≠ 0 

Which grantee that there is pitch fork bifurcation at 𝐴1 where 𝑑`1 = 𝑑1. 



IHJPAS. 2025, 38(2) 

320 
 

Theorem (2): System (1) with the parameter value  �̅�2 = 𝑑2 =
𝑙2�̅�2

1+𝑎3�̅�2
 ,has a transcritical and 

pitchfork bifurcation at 𝐴2 = (0. �̅�2. 0). if the next conditions hold: 

    Z̅1 ≠ �̅�2.                                                                                                                                                (9) 

 Z̅3 ≠ �̅�4 .                                                                                                                                (10)  

Where: 

Z̅1 =
𝑙2∆

(1+a3R2)2(1+a4R3)2
 , �̅�2 = (

𝑙2𝑅3𝑎3∆2

(1+a3R2)3(1+a4R3)
+

𝑙2𝑅2𝑎4

(1+a3R2)(1+a4R3)3
) , 

�̅�3 =
𝑙2

(1+a3R2)(1+a4R3)
(

𝑅3𝑎3
2∆2

(1+a3R2)3
+

𝑅2𝑎4
2

(1+a4R3)3
) , �̅�4 =

𝑙2∆

(1+a3R2)2(1+a4R3)2
(

𝑎4

(1+𝑎4𝑅3)
+

𝑎3∆

(1+a3R2)
). 

Proof: By using the Jacobian matrix in equation  (1.8) in (20) 𝐽2̅ = 𝐽2(𝐴2. �̅�2) = [�̅�𝑖𝑗]3×3
, 

where �̅�𝑖𝑗 = 𝑎𝑖𝑗 , except �̅�33 = 0.Then, the characterizing equation for𝐽2̅ has an eigenvalue of 

zero, which is λ2R3
 at  �̅�2 = 𝑑2 . 

Now, let 𝑇[2] = (𝑡1
[2]

. 𝑡2
[2]

. 𝑡3
[2]

)𝑇  be the eigenvector associated with an eigenvalue  λ2R2
= 0. 

Thus,(𝐽2̅ −  λ2R2
𝐼)𝑇[2] = 0.  this gives:𝑡1

[2]
= 0.  𝑡2

[2]
= ∆𝑡3

[2]
. 

 where ∶ 

∆= −(
𝑟2𝑅2𝐾2

𝑟2−𝑑1
+

𝛽2𝑅2

(1+𝑎3𝑅2)(𝑟2−𝑑1)
), 𝑎𝑛𝑑 𝑡3

[2]
 any real number that is not zero. 

Let 𝑀[2] = (𝑚1
[2]

. 𝑚2
[2]

. 𝑚3
[2]

)𝑇 be the eigenvector associated with an eigenvalue  λ2R2
= 0 of 

the matrix 𝐽2̅
 𝑇.Then.  (𝐽2̅

 𝑇 −  λ2R2
𝐼) 𝑀[2] = 0 . 

By solving this equation for . 𝑀[2] = (0. 0.𝑚3
[2]

)
𝑇

.  and 𝑚3
[2]

 any real number that is not zero. 

Now consider this: 
𝜕𝑓

𝜕𝑑2
= 𝑓𝑑2

(𝑋. 𝑑2) = (
𝜕𝑓1

𝜕𝑑2
.
𝜕𝑓2

𝜕𝑑2
.
𝜕𝑓3

𝜕𝑑2
)
𝑇

= (0. 0. −𝑅3)
𝑇. 

So, 𝑓𝑑2
(𝐴2. �̅�2) = (0. 0. 0)𝑇 and hence (𝑀[2])

𝑇
𝑓𝑑2

(𝐴2. �̅�2) = 0. 

Using Sotomayor's theorem, it is impossible to satisfy the saddle-node bifurcation condition. 

The first condition for transcritical bifurcation is therefore satisfied. Now 

𝐷𝑓𝑑(𝑋. 𝑑2) = [
0     0     0
0     0     0
0     0 − 1

].  

where, 𝐷𝑓𝑑2
(𝑋. 𝑑2) represents the derivative of 𝑓𝑑2

(𝑋. 𝑑2) with respect to 𝑋 = (𝑅1. 𝑅2. 𝑅3)
𝑇 . 

Furthermore, it is observed that: 

𝐷𝑓𝑑2
(𝐴2. 𝑑2)𝑇

[2] = [
0     0     0
0     0     0
0     0 − 1

] [

0

∆𝑡3
[2]

𝑡3
[2]

] = [

0
0

−𝑡3
[2]

] 

(𝑀[2])
𝑇
[𝐷𝑓𝑑2

(𝐴2. �̅�2)𝑇
[2]] = (0. 0.𝑚3

[2]
) (0. 0. −𝑡3

[2]
)
𝑇

= −𝑚3
[2]

𝑡3
[2]

≠ 0 

By substituting 𝑇[2]𝑖𝑛  equation (3) we get: 

𝐷2𝐹𝜇(𝐴2. �̅�2)(𝑇
[2]. 𝑇[2]) = (�̅�𝑖𝑗) .                                                                                                   (11) 

�̅�11 = 2∆2𝑅1𝑡3
[2]

(
𝑟1𝐾1

2

(1+𝐾1𝑅2)3
+

𝛽1𝑎2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)3
), 

�̅�21 = 2 𝑡3
[2]

[(
R3𝑎3∆2

(1+a3R2)2
+

 a4𝑅2

(1+a4R3)2
−

∆

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
)

𝛽2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

𝑙1R1𝑎2∆2

(1+a1R1)(1+a2R2)3
+

(∆ −
𝐾2𝑅2

(1+K2R3)
)

r2K2

(1+K2R3)2
− 𝑚2∆

2], 

�̅�31 = 2𝑡3
[2] 𝑙2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
(

1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

𝑅3𝑎3∆2

(1+a3R2)2
−

𝑅2𝑎4

(1+a4R3)2
), 

Hence, it was obtained. 



IHJPAS. 2025, 38(2) 

321 
 

(𝑀[2])
𝑇
[𝐷2𝐹𝜇(𝐴2. �̅�2)(𝑇

[2]. 𝑇[2])] = 2 (𝑡3
[2]

)
2

𝑚3
[2]

(Z̅1 − �̅�2) ≠ 0. 

This indicates that the system (1) exhibits a transcritical bifurcation at  𝐴2 with a parameter 

 �̅�2 = 𝑑2. If condition (9) not satisfied, then. 

By substituting 𝑇[2] in equation (4) we obtain: 

𝐷3𝐹𝜇(𝐴2. �̅�2)(𝑇
[2]. 𝑇[2]. 𝑇[2]) = (�̅�𝑖𝑗).                                                                                         (12) 

�̅�11 = −6∆2 (𝑡3
[2]

)
3

(
𝛽1𝑎2

2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)4
+

𝑟1𝐾1
3𝑅1

(1+𝐾1𝑅2)4
), 

�̅�21 = 6(𝑡3
[2]

)
3

[(∆ −
K2𝑅2

(1+K2R3)
)

r2K2
2

(1+K2R3)3
+ (

𝑙1∆3

(1+a2R2)3
−

𝛽1

(1+𝑎2𝑅2)3
)

𝑎2
2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+

(
𝑎4∆

(1+𝑎3𝑅2)(1+𝑎4𝑅3)
+

 a4

(1+a4R3)2
−

R3𝑎3
2∆2

(1+a3R2)2
)

𝛽2∆

(1+𝑎3𝑅2)2(1+𝑎4𝑅3)
], 

�̅�31 =
6(𝑡3

[2]
)
3
𝑙2

(1+a3R2)(1+a4R3)
(

𝑅3𝑎3
2∆2

(1+a3R2)3
+

𝑅2𝑎4
2

(1+a4R3)3
−

𝑎4∆

(1+a3R2)(1+a4R3)2
−

𝑎3∆2

(1+a3R2)2(1+a4R3)
). 

(𝑀[2])
𝑇
[𝐷3𝐹𝜇(𝐴2. �̅�2)(𝑇

[2]. 𝑇[2]. 𝑇[2])] = 6 (𝑡3
[2]

)
3

𝑚3
[2]

(�̅�3 − �̅�4) ≠ 0 

Which grantee that there is pitch fork bifurcation at 𝐴2 where �̅�2 = 𝑑2. 

Theorem (3):  System (1) with the parameter value  �̂�2 = 𝑑2 =
𝑙2�̂�2

1+𝑎3�̂�2
 , has a transcritical and 

pitchfork bifurcation at 𝐴3 = (�̂�1. �̂�2. 0), if the next conditions hold: 

�̂�1 ≠ �̂�2.                                                                                                                                                 (13) 

�̂�3 ≠ �̂�4.                                                                                                                                                 (14) 

Where: 

�̂�1 =
𝑙2∅2

(1+a3R2)2(1+a4R3)2
 ، �̂�2 =

𝑙2𝑅3𝑎3∅2
2

(1+a3R2)3(1+a4R3)
+

𝑙2𝑅2𝑎4

(1+a3R2)(1+a4R3)3
)  ،  

�̂�3 =
𝑙2

(1+a3R2)(1+a4R3)
(

𝑅3𝑎3
2∅2

3

(1+a3R2)3
+

𝑅2𝑎4
2

(1+a4R3)3
) ،  �̂�4 =

𝑙2∅2

(1+a3R2)2(1+a4R3)2
(

𝑎4

(1+𝑎4𝑅3)
+

𝑎3∅2

(1+a3R2)
). 

Proof: By using the Jacobian matrix inequation (1.9) in (20)  𝐽3 = 𝐽3(𝐴3. �̂�2) = [�̂�𝑖𝑗]3×3
, 

where �̂�𝑖𝑗 = 𝑑𝑖𝑗 , except �̂�33 = 0 .Then, the characterizing equation for 𝐽3 has an eigenvalue 

of zero, which is λ3R3
 at  �̂�2 = 𝑑2 . 

Now, let 𝑇[3] = (𝑡1
[3]

. 𝑡2
[3]

. 𝑡3
[3]

)𝑇 be the eigenvector associated with an eigenvalue  λ3R3
= 0. 

Thus,(𝐽3 −  λ3R3
𝐼)𝑇[3] = 0.  this gives: 𝑡1

[3]
= ∅1𝑡3

[3]
. 𝑡2

[3]
= ∅2𝑡3

[3]
. 

where: ∅1 =
𝑑23𝑑12

𝑑11𝑑22−𝑑12𝑑21
 ، ∅2 =

𝑑23𝑑11

𝑑12𝑑21−𝑑11𝑑22
. 

and 𝑡3
[3]

 any real number that is not zero. 

Let 𝑀[3] = (𝑚1
[3]

. 𝑚2
[3]

. 𝑚3
[3]

)𝑇 be the eigenvector associated with an eigenvalue  λ3R3
= 0 of 

the matrix 𝐽3
 𝑇.Then.  (𝐽3

 𝑇 −  λ3R3
𝐼)𝑀[3]  = 0 . 

By solving this equation for . 𝑀[3] = (0. 0.𝑚3
[3]

)
𝑇

,and 𝑚3
[3]

 any real number that is not zero. 

Now consider this:
𝜕𝑓

𝜕𝑑2
= 𝑓𝑑2

(𝑋. 𝑑2) = (
𝜕𝑓1

𝜕𝑑2
.
𝜕𝑓2

𝜕𝑑2
.
𝜕𝑓3

𝜕𝑑2
)
𝑇

= (0. 0. −𝑅3)
𝑇 . 

So, 𝑓𝑑2
(𝐴3. �̂�2) = (0. 0. 0)𝑇and hence (𝑀[3])

𝑇
𝑓𝑑2

(𝐴3. �̂�2) = 0. 

Using Sotomayor's theorem, it is impossible to satisfy the saddle-node bifurcation condition. 

The first condition for transcritical bifurcation is therefore satisfied. Now 

𝐷𝑓𝑑2
(𝑋. 𝑑2) = [

0     0     0
0     0     0
0     0 − 1

].  



IHJPAS. 2025, 38(2) 

322 
 

where, 𝐷𝑓𝑑2
(𝑋. 𝑑2) represents the derivative of 𝑓𝑑2

(𝑋. 𝑑2) with respect to 𝑋 = (𝑅1. 𝑅2. 𝑅3)
𝑇 . 

Furthermore, it is observed that: 

𝐷𝑓𝑑2
(𝐴3. 𝑑2)𝑇

[3] = [
0     0     0
0     0     0
0     0 − 1

] [

∅1𝑡3
[3]

∅2𝑡3
[3]

𝑡3
[3]

] = [

0
0

−𝑡3
[3]

] 

(𝑀[3])
𝑇
[𝐷𝑓𝑑2

(𝐴3. �̂�2)𝑇
[3]] = (0. 0.𝑚3

[3]
) (0. 0. −𝑡3

[3]
)
𝑇

= −𝑚3
[3]

𝑡3
[3]

≠ 0 

By substituting 𝑇[3] in equation (3) we get: 

𝐷2𝐹𝜇(𝐴3. �̂�2)(𝑇
[3]. 𝑇[3]) = (�̂�𝑖𝑗).                                                                                                   (15) 

�̂�11 = 2(𝑡3
[3]

)
2

[(
𝑅2𝑎1∅1

2

(1+𝑎2𝑅2)2
−

∅1∅2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+

𝑎2𝑅2∅2
2

(1+𝑎2𝑅2)2
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+

(
𝐾1𝑅1∅1

2

(1+𝐾1𝑅2)
−∅1∅2)

𝑟1𝐾1

(1+𝐾1𝑅2)2
− 𝑚1∅1

2], 

�̂�21 = 2(𝑡3
[3]

)
2

[(
∅1∅2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

R2𝑎1∅1
2

(1+a1R1)2
−

R1𝑎2∅2
2

(1+a2R2)2
)

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (

a3R3∅2
2

(1+a3R2)2
+

∅2

(1+a3R2)(1+a4R3)
+

 a4𝑅2

(1+a4R3)2
)

𝛽2

(1+a3R2)(1+a4R3)
−𝑚2∅2

2 + (
𝐾2𝑅2

(1+K2R3)
− ∅2)

r2K2

(1+K2R3)2
], 

�̂�31 =
2(𝑡3

[3]
)
2
𝑙2

(1+a3R2)(1+a4R3)
(

∅2

(1+a3R2)(1+a4R3)
−

𝑅3𝑎3∅2
2

(1+a3R2)2
 −

𝑅2𝑎4

(1+a4R3)2
). 

Hence, it was obtained 

(𝑀[3])
𝑇
[𝐷2𝐹𝜇(𝐴3. �̂�2)(𝑇

[3]. 𝑇[3])] = 2 (𝑡3
[3]

)
2

𝑚3
[3]

(�̂�1 − �̂�2) ≠ 0. 

This indicates that the system (1) exhibits a transcritical bifurcation at 𝐴2 with a 

parameter �̂�2 = 𝑑2. If condition (13) not satisfied then. 

By substituting 𝑇[3]  in equation (4) we get: 

𝐷3𝐹𝜇(𝐴3. �̂�2)(𝑇
[3]. 𝑇[3]. 𝑇[3]) = (�̂�𝑖𝑗),                                                                                        (16) 

�̂�11 = 6(𝑡3
[3]

)
3

[(
𝑎2∅1∅2

2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)2
−

∅2
3𝑎2

2𝑅1

(1+𝑎2𝑅2)3
−

R2𝑎1
2∅1

3

(1+𝑎1𝑅1)3
𝑎1∅1

2∅2

(1+𝑎1𝑅1)2(1+𝑎2𝑅2)
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (∅1 −

𝐾1𝑅1∅1

(1+𝐾1𝑅2)
)

𝑟1𝐾1
2∅2

2

(1+𝐾1𝑅2)3
] , 

�̂�21 = 6(𝑡3
[3]

)
3

[(∅2 −
K2𝑅2

(1+K2R3)
)

r2K2
2

(1+K2R3)3
−

𝛽1𝑎2
2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)4
+ (

R2𝑎1
2∅1

3

(1+a1R1)3
+

𝑎1∅1
2∅2

(1+a1R1)2(1+a2R2)
−

𝑎2∅1∅2
2

(1+a1R1)(1+a2R2)2
+

R1𝑎2
2∅2

3

(1+a2R2)3
) 

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (

𝑎4∅2

(1+𝑎3𝑅2)(1+𝑎4𝑅3)2
+

 a4∅2
2

(1+a3R2)(1+a4R3)2
−

R3𝑎3
2∅2

3

(1+a3R2)3
)

𝛽2𝐼

(1+𝑎3𝑅2)2(1+𝑎4𝑅3)
], 

�̂�31 =
6(𝑡3

[3]
)
3
𝑙2

(1+a3R2)(1+a4R3)
(

𝑅3𝑎3
2∅2

3

(1+a3R2)3
+

𝑅2𝑎4
2

(1+a4R3)3
−

𝑎4∅2

(1+a3R2)(1+a4R3)2
−

𝑎3∅2
2

(1+a3R2)2(1+a4R3)
). 

Hence,  

(𝑀[3])
𝑇
[𝐷3𝐹𝜇(𝐴3. �̂�2)(𝑇

[3]. 𝑇[3]. 𝑇[3])] = 6 (𝑡3
[3]

)
3

𝑚3
[3]

(�̂�3 − �̂�4) ≠ 0. 

 By condition (14). Which grantee that there is pitch fork bifurcation at 𝐴3 where �̅�2 = 𝑑2. 

Theorem (4):   System (1) with the parameter value �̿�1 = 𝑟1 =
𝛽1�̿�2(1+𝐾1�̿�2)

(1+𝑎2�̿�2)
, has a transcritical 

and pitchfork bifurcation at 𝐴4 = (0. �̿�2. �̿�3), if the next conditions hold: 

�̿�1 ≠ �̿�2.                                                                                                                                                (17) 

�̿�3 ≠ �̿�4.                                                                                                                                                (18) 

Where: 

�̿�1 =
𝑟1𝐾1

2𝐼1
2

(1+𝐾1𝑅2)3
+

𝛽1𝑎2𝐼1
2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)3
+

𝛽1𝑅2𝑎1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)3
, 



IHJPAS. 2025, 38(2) 

323 
 

�̿�2 = (𝑚1 +
𝑟1𝐾1𝐼1

(1+𝐾1𝑅2)2
+

𝛽1𝐼1

(1+𝑎1𝑅1)2(1+𝑎2𝑅2)2
), 

�̿�3 = (
𝑟1𝐾1

2

(1+𝐾1𝑅2)3
𝐼1

2 +
𝛽1𝑎1

(1+𝑎1𝑅1)3(1+𝑎2𝑅2)2
𝐼1+

𝛽1𝑎2

(1+𝑎1𝑅1)2(1+𝑎2𝑅2)3
𝐼1

2 −
β1R2𝑎1

2

(1+𝑎1𝑅1)4(1+a2R2)
), 

�̿�4 = (
β1R2𝑎1

2

(1+𝑎1𝑅1)4(1+a2R2)
+

𝑟1𝐾1
3𝑅1

(1+𝐾1𝑅2)4
𝐼1

3 +
𝛽1𝑎2

2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)4
𝐼1

3). 

Proof: By using the Jacobian matrix from equation (1.10) in (20)  𝐽4̿ = 𝐽4(𝐴4. �̿�1) = [�̿�𝑖𝑗]3×3
, 

where �̿�𝑖𝑗 = 𝑒𝑖𝑗 . except �̿�11 = 0 .Then, the characterising equation for 𝐽4̿ has an eigenvalue of 

zero, which is λ4R1
= 0 at �̿�1 = 𝑟1. 

Now, let 𝑇[4] = (𝑡1
[4]

. 𝑡2
[4]

. 𝑡3
[4]

)𝑇 be the eigenvector associated with an eigenvalue λ4R1
= 0. 

Thus,(𝐽4̿ −  λ4R1
𝐼)𝑇[4] = 0.  this gives: 

𝑡2
[4]

= 𝐼1𝑡1
[4]

. 𝑡3
[4]

= 𝐼2𝑡1
[4]

. 

where: 

 𝐼1 =
𝑒21𝑒33

𝑒23𝑒32−𝑒22𝑒33
 ,  𝐼2 =

𝑒32𝑒21

𝑒33𝑒22−𝑒23𝑒32
  . 

and 𝑡1
[4]

 any real number that is not zero. 

Let 𝑀[4] = (𝑚1
[4]

. 𝑚2
[4]

. 𝑚3
[4]

)𝑇 be the eigenvector associated with an eigenvalue  λ4R1
= 0 of 

the matrix 𝐽4̿
 𝑇.Then.  (𝐽4̿

 𝑇 −  λ4R1
𝐼) 𝑀[4] = 0 . 

By solving this equation for . 𝑀[4] = (𝑚1
[4]

. 0. 0)
𝑇

and 𝑚1
[4]

 any real number that is not zero. 

Now consider this: 
𝜕𝑓

𝜕𝑟1
= 𝑓𝑟1(𝑋. 𝑟1) = (

𝜕𝑓1

𝜕𝑟1
.
𝜕𝑓2

𝜕𝑟1
.
𝜕𝑓3

𝜕𝑟1
)
𝑇

= (
𝑅1

1+𝐾1𝑅2
. 0. 0)

𝑇

. 

So.  𝑓𝑟1(𝐴4. �̿�1) = (0. 0. 0)𝑇and hence  (𝑀[4])𝑇𝑓𝑟1(𝐴4. �̿�1) = 0. 

Using Sotomayor's theorem, it is impossible to satisfy the saddle-node bifurcation condition. 

The first condition for transcritical bifurcation is therefore satisfied. Now 

𝐷𝑓𝑟1(𝑋. 𝑟1) = [ 

1

1 + 𝐾1�̿�2

  0     0

      0           0      0 
     0           0      0

].  

where, 𝐷𝑓𝑟1(𝑋. 𝑟1) represents the derivative of 𝑓𝑟1(𝑋. 𝑟1) with respect to 𝑋 = (𝑅1. 𝑅2. 𝑅3)
𝑇 . 

Furthermore, it is observed that: 

𝐷𝑓𝑟1(𝐴4. 𝑟1)𝑇
[4] = [

1

1 + 𝐾1𝑅2
  0     0

      0           0      0 
     0           0      0

] [

𝑡1
[4]

𝐼1𝑡1
[4]

𝐼2𝑡1
[4]

] =

[
 
 
 𝑡1

[4]

1 + 𝐾1�̿�2

0
0 ]

 
 
 

 

(𝑀[4])
𝑇
[𝐷𝑓𝑟1(𝐴4. �̿�1)𝑇

[4]] = (𝑚1
[4]

. 0. 0) (
𝑡1
[4]

1+𝐾1�̿�2
. 0. 0)

𝑇

=
𝑚1

[4]
𝑡1
[4]

1+𝐾1�̿�2
≠ 0, 

By substituting 𝑇[4] in equation (3) we get: 

𝐷2𝐹𝜇(𝐴4. �̿�1)(𝑇
[4]. 𝑇[4]) = (�̿�𝑖𝑗),                                                                                                    (19) 

�̿�11 = 2(𝑡1
[4]

)
2

[(
𝑅2𝑎1𝐼1

2

(1+𝑎2𝑅2)2
−

𝐼1𝐼2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+

𝑎2𝑅2𝐼2
2

(1+𝑎2𝑅2)2
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
− 𝑚1𝐼1

2 +

(
𝐾1𝑅1𝐼1

2

(1+𝐾1𝑅2)
−𝐼1𝐼2)

𝑟1𝐾1

(1+𝐾1𝑅2)2
], 

�̿�21 = 2(𝑡1
[4]

)
2

[(
𝐼1𝐼2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

R2𝑎1𝐼1
2

(1+a1R1)2
−

R1𝑎2𝐼2
2

(1+a2R2)2
)

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (

𝐾2𝑅2

(1+K2R3)
−

𝐼2)
r2K2

(1+K2R3)2
+ (

a3R3𝐼2
2

(1+a3R2)2
+

𝐼2

(1+a3R2)(1+a4R3)
+

 a4𝑅2

(1+a4R3)2
)

𝛽2

(1+a3R2)(1+a4R3)
−𝑚2𝐼2

2], 



IHJPAS. 2025, 38(2) 

324 
 

�̿�31 =
2(𝑡1

[4]
)2 𝑙2

(1+a3R2)(1+a4R3)
(

𝐼2

(1+a3R2)(1+a4R3)
−

𝑅3𝑎3𝐼2
2

(1+a3R2)2
 −

𝑅2𝑎4

(1+a4R3)2
). 

Hence, it was obtained  

(𝑀[4])
𝑇
[𝐷2𝐹𝜇(𝐴4. �̿�1)(𝑇

[4]. 𝑇[4])] = 2 (𝑡1
[4]

)
2

𝑚1
[4]

(𝑍1 − 𝑍2) ≠ 0 . 

This indicates that the system (1) exhibits a transcritical bifurcation at 𝐴4 with a parameter 

�̿�1 = 𝑟1, and no pitch fork bifurcation at 𝐴4 where �̿�1 = 𝑟1. 

By substituting 𝑇[4] in equation (4) we get: 

𝐷3𝐹𝜇(𝐴4. �̿�1)(𝑇
[4]. 𝑇[4]. 𝑇[4]) = (�̿�𝑖𝑗),                                                                                          (20) 

�̿�11 = 6(𝑡1
[4]

)
3

[(
𝑎2𝐼1𝐼2

2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)2
−

𝐼2
3𝑎2

2𝑅1

(1+𝑎2𝑅2)3
−

R2𝑎1
2𝐼1

3

(1+𝑎1𝑅1)3
+

𝑎1𝐼1
2𝐼2

(1+𝑎1𝑅1)2(1+𝑎2𝑅2)
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (𝐼1 −

𝐾1𝑅1𝐼1

(1+𝐾1𝑅2)
)

𝑟1𝐾1
2𝐼2

2

(1+𝐾1𝑅2)3
], 

�̿�21 = 6(𝑡1
[4]

)
3

[(𝐼2 −
K2𝑅2

(1+K2R3)
)

r2K2
2

(1+K2R3)3
−

𝛽1𝑎2
2𝑅1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)4
+ (

R2𝑎1
2𝐼1

3

(1+a1R1)3
+

𝑎1𝐼1
2𝐼2

(1+a1R1)2(1+a2R2)
−

𝑎2𝐼1𝐼2
2

(1+a1R1)(1+a2R2)2
+

R1𝑎2
2𝐼2

3

(1+a2R2)3
) 

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (

𝑎4𝐼2

(1+𝑎3𝑅2)(1+𝑎4𝑅3)2
+

 a4𝐼2
2

(1+a3R2)(1+a4R3)2
−

R3𝑎3
2𝐼2

3

(1+a3R2)3
)

𝛽2𝐼

(1+𝑎3𝑅2)2(1+𝑎4𝑅3)
], 

�̿�31 = 6(𝑡1
[4]

)
3 𝑙2

(1+a3R2)(1+a4R3)
(

𝑅3𝑎3
2𝐼2

3

(1+a3R2)3
+

𝑅2𝑎4
2

(1+a4R3)3
−

𝑎4𝐼2

(1+a3R2)(1+a4R3)2
−

𝑎3𝐼2
2

(1+a3R2)2(1+a4R3)
). 

Hence,  

(𝑀[4])
𝑇
[𝐷3𝐹𝜇(𝐴4. �̿�1)(𝑇

[4]. 𝑇[4]. 𝑇[4])] = 6 (𝑡1
[4]

)
3

𝑚1
[4]

(�̿�3 − �̿�4) ≠ 0 

If condition (17). Which grantee that there is pitch fork bifurcation at 𝐴4 where �̿�1 = 𝑟1. 

Theorem (5): Suppose that conditions (1.11b), (1.11e) and, reverse condition (1.11c) in (17) 

with the following conditions are satisfied: 

𝑟11𝑟22 < 𝑟12𝑟21.                                                                                                                                    (21) 

�̃�1 ≠ �̃�2.                                                                                                                                                 (22) 

Where, 

�̃�1 =
𝑙2𝛾2

(1+a3R2)2(1+a4R3)2
،�̃�2 = (

𝑙2𝑅3𝑎3𝛾2
2

(1+a3R2)3(1+a4R3)
+

𝑙2𝑅2𝑎4

(1+a3R2)(1+a4R3)3
)  ، 𝛾2 =

−�̃�33

�̃�32
< 0. 

Then the system near 𝐴5 = (�̃�1. �̃�2. �̃�3). has a saddle–node bifurcation at  

�̃�2 = 𝑑2 =
�̃�11�̃�32�̃�23

�̃�11�̃�22−�̃�12�̃�21
  . 

Proof: By using the Jacobian matrix in equation (1.11) in (20)   𝐽5 = 𝐽5(𝐴5. �̃�2) = [�̃�𝑖𝑗]3×3
, 

where �̃�𝑖𝑗 = 𝑟𝑖𝑗 .  except �̃�33 =
𝑙2�̃�2

(1+𝑎3�̃�2)(1+𝑎4�̃�3)2
− �̃�2 .Then, the characterising equation for 𝐽5 

has an eigenvalue of zero, which is (say λ5R3
= 0) if and if 𝜌3 = 0. 

Now, let 𝑇[5] = (𝑡1
[5]

. 𝑡2
[5]

. 𝑡3
[5]

)𝑇 be the eigenvector associated with an eigenvalue  λ5R3
= 0. 

Thus,(𝐽5 −  λ5R3
𝐼)𝑇[5] = 0.  this gives: 𝑡1

[5]
= 𝛾1𝑡3

[5]
.  𝑡2

[5]
= 𝛾2𝑡3

[5]
 .  and 𝑡3

[5]
 any real number 

that is not zero.  𝐴𝑛𝑑 𝛾1 =
�̃�12�̃�33

�̃�11�̃�32
. 

Let 𝑀[5] = (𝑚1
[5]

. 𝑚2
[5]

. 𝑚3
[5]

)𝑇 be the eigenvector associated with an eigenvalue λ5R3
= 0 of 

the matrix 𝐽5
 𝑇.Then .  (𝐽5

 𝑇 −  λ5R3
𝐼) 𝑀[5] = 0 . 

By solving this equation for . 𝑀[5] = (𝛾3𝑚3
[5]

. 𝛾4𝑚3
[5]

. 𝑚3
[5]

)
𝑇

. where 𝑚3
[5]

 any real number that 

is not zero, and 𝛾3 =
�̃�21�̃�33

�̃�11�̃�23
> 0.       𝛾4 =

−�̃�33

�̃�23
< 0.  



IHJPAS. 2025, 38(2) 

325 
 

Now consider this: 
𝜕𝑓

𝜕𝑑2
= 𝑓𝑑2

(𝑋. 𝑑2) = (
𝜕𝑓1

𝜕𝑑2
.
𝜕𝑓2

𝜕𝑑2
.
𝜕𝑓3

𝜕𝑑2
)
𝑇

= (0. 0. −𝑅3)
𝑇.  

So.  𝑓𝑑2
(𝐴5. �̃�2) = (0. 0. −�̃�3)

𝑇
and hence  (𝑀[5])

𝑇
𝑓𝑑2

(𝐴5. �̃�2) = −�̃�3𝑚3
[5]

≠ 0. 

By substituting 𝑇[5] in equation  (3) we get: 

𝐷2𝐹𝜇(𝐴5. �̃�2)(𝑇
[5]. 𝑇[5]) = (�̃�𝑖𝑗).                                                                                                   (23) 

�̃�11 = 2(𝑡3
[5]

)
2

[(
𝐾1𝑅1𝛾1

2

(1+𝐾1𝑅2)
−𝛾1𝛾2)

𝑟1𝐾1

(1+𝐾1𝑅2)2
− 𝑚1𝛾1

2 + (
𝑅2𝑎1𝛾1

2

(1+𝑎2𝑅2)2
−

𝛾1𝛾2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+

𝑎2𝑅2𝛾2
2

(1+𝑎2𝑅2)2
)

𝛽1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
], 

�̃�21 = 2(𝑡3
[5]

)
2

[(
𝛾1𝛾2

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
−

R2𝑎1𝛾1
2

(1+a1R1)2
−

R1𝑎2𝛾2
2

(1+a2R2)2
)

𝑙1

(1+𝑎1𝑅1)(1+𝑎2𝑅2)
+ (

𝐾2𝑅2

(1+K2R3)
−

𝛾2)
r2K2

(1+K2R3)2
+ (

a3R3𝛾2
2

(1+a3R2)2
+

𝛾2

(1+a3R2)(1+a4R3)
+

 a4𝑅2

(1+a4R3)2
)

𝛽2

(1+a3R2)(1+a4R3)
−𝑚2𝛾2

2], 

�̃�31 =
2(𝑡3

[5]
)
2
𝑙2

(1+a3R2)(1+a4R3)
(

𝛾2

(1+a3R2)(1+a4R3)
−

𝑅3𝑎3𝛾2
2

(1+a3R2)2
 −

𝑅2𝑎4

(1+a4R3)2
). 

Hence, it was obtained. 

(𝑀[5])
𝑇
[𝐷2𝐹𝜇(𝐴5. �̃�2)(𝑇

[5]. 𝑇[5]. 𝑇[5])] = 2 (𝑡3
[5]

)
2

𝑚3
[5]

(�̃�1 − �̃�2) ≠ 0. 

So, by condition (21) system (1) exhibits a saddle-node bifurcation at 𝐴5 for any value of the 

parameter �̃�2 = 𝑑2 but no pitch fork bifurcation at 𝐴5 where �̃�2 = 𝑑2. 

 

3. Numerical simulation 

     In this section, the dynamical behavior of the system (1) was investigated. Calculations may 

be performed for a one set of parameters with a different starting point to examine the analytical results 

and understand how the parameters impact the dynamic model. Figure 1. (a-d), shows that system (1) 

has positive solution. 

 𝑟1 = 0.5.𝑚1 = 0.5. 𝐾1 = 0.003. 𝛽1 = 0.3. 𝑙1 = 0.3. 𝑎1 = 0.1. 𝑎2 = 0.02. 𝑑1 = 0.1. 𝑟2 = 0.5.        

𝑚2 = 0.4. 𝐾2 = 0.07. 𝛽2 = 0.7. 𝑙2 = 0.7. 𝑎3 = 0.01. 𝑎4 = 0.02. 𝑑2 = 0.4.                          (24) 

 

 

Figure 𝟏. (𝐚 − 𝐝). Time series of system solution (1) begin with different starting points (0.5. 0.6. 0.7).  

(0.4. 0.5. 0.8).  and  (0.2. 0.9. 0.1) . (𝐛) Path of  R1  depending on time, (𝐜) Path of   R2  depending on time,(𝐝) 

Path of  R3 depending on time. 

0.2
0.3

0.4
0.5

0.6
0.7

0.8
0.9

0.4
0.5

0.6
0.7

0.8
0.9

1
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

 

prey

a

intermediate predator
 

t
o

p
 p

r
e
d

a
t
o

r

prey 

intermediate predator

top predator

(0.6748,0.5691,0.4943)
(0.2,0.9,0.1)

(0.5,0.6,0.7)

(0.4,0.5,0.8)

0 100 200 300 400 500 600 700 800 900 1000
0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

b

Time

p
r
e
y

 

 

R
1

R
2

R
3

0 100 200 300 400 500 600 700 800 900 1000

0.4

0.5

0.6

0.7

0.8

0.9

1

Time

in
te

r
m

e
d

ia
te

 p
r
e
d

a
to

r

c

 

 
R

1

R
2

R
3

0 100 200 300 400 500 600 700 800 900 1000
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8
d

Time

to
p

 p
r
e

d
a

to
r

 

 

R
1
 

R
2
 

R
3



IHJPAS. 2025, 38(2) 

326 
 

Now, to study the effect of parameters on the dynamic behavior of the system (1), the system 

(1) was numerically solved to the data in (24) by changing a single parameter each time the 

results are given in Table (2). 

The result of changing the value of the parameter  𝑟1 in the range  

0.15 ≤ 𝑟1 < 3.7 solution approaches  A4, as described in Figure 2 ,(a)for typical value  𝑟1 = 2 

. in the range 0.1 ≤ 𝑟1 < 0.15 the solution approaches to  A5, as described in Figure 2,(b) for 

typical value 𝑟1 = 0.14 . 

 
 

Figure 2. (a).(T.S.) of system solution (1) for the given values in the equation (24) with 𝑟1 = 2   ,which 

approaches to 𝐴4 = (0. 0.5716. 0.2356). (b). (T.S.) of system solution (1) for the given values in the 

equation (24) with 𝑟1 = 0.14, which approaches to  A5 = (2.25226. 0.5627. 1.0574) . 

The impact of varying the parameter  𝑑1  in the range 0.001 ≤ 𝑑1 < 0.475 the solution 

approaches to  A5, as shown in Figure 3, (a) for typical value 𝑑1 = 0.2 . in the range  0.475 ≤

𝑑1 < 0.799   solution approaches   A3 , as described in figure 3, (b) for typical value 𝑑1 =

0.4799 . in the range 0.799 ≤ 𝑑1 < 1  solution approaches A1, as shown in Figure 3, (c)   for 

typical value 𝑑1 = 0.18. 

 

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

Figure 3. (a),(T.S.) of system solution (1) for the given values in the equation  (24) with 

𝑑1 = 0.2 ,which approaches to  A5 = (0.6981. 0.53. 0.19). (b) (T.S.). of system solution (1) for the given values 

in the equation (24) with 𝑑1 = 0.4799 ,which approaches to   A3 = (0.6932. 0.5383. 0). (c) (T.S.) of system 

solution (1) for the given values in the equation (24) with 𝑑1 = 0.81   ,which approaches to A1 = (0.999. 0. 0). 

 

0 100 200 300 400 500 600 700 800 900 1000
0

0.5

1

1.5

2

2.5

3

Time

p
o

p
u

la
ti

o
n

 

 

prey

intermediate predator

top predator

0 100 200 300 400 500 600 700 800 900 1000
0

0.5

1

1.5

2

2.5

3

Time

p
o

p
u

la
ti

o
n

 

 

prey

intermediate predator

top predator

0 100 200 300 400 500 600 700 800 900 1000
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Time

p
o

p
u

la
ti

o
n

 

 

prey

intermediate predator

top  predator

0 100 200 300 400 500 600 700 800 900 1000
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time

p
o

p
u

la
ti

o
n

 

 

prey

intermediate predator

top predator

0 100 200 300 400 500 600 700 800 900 1000
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time

p
o

p
u

la
ti

o
n

 

 

prey

intermediate predator

top predator



IHJPAS. 2025, 38(2) 

327 
 

For the parameter  𝑑2  in the range 0.1 ≤ 𝑑2 < 0.869 solution approaches  A5, as shown in 

Figure 4,(a) for typical value   𝑑2 = 0.5 ,the range growing further  0.869 ≤ 𝑑2 < 1 solution 

approaches  A3 as described in Figure 4,(b) for typical value   𝑑2 = 0.921 , in the range 

 

  
Figure 4. (a) (T.S.) of system solution (1) for the given values in the equation (24) with  d2 = 0.5,which 

approaches to A5 = (0.5875. 0.7139. 0.3857). (b) (T.S.) of system solution (1) for the given values in the 

equation (24)  with  d2 = 0.921 ,which approaches to   A3 = (0.2753. 1.2058. 0).   

 

Table 2. The numerical result. 

Range of parameter Stable The bifurcation point 

0.1 ≤ 𝑟1 < 0.15 

0.15 ≤ 𝑟1 < 3.7 
𝐴4 

𝐴5 

𝑟1 = 0.15 

0.0001 ≤ 𝑟2 < 0.1 

0.1 ≤ 𝑟2 < 0.5 
𝐴2 

𝐴5 

𝑟2 = 0.1 

0.003 ≤ 𝐾1 ≤ 0.999 

0.999 ≤ 𝐾1 < 1 
𝐴5 

𝐴4 

𝐾1 = 0.999 

0.00007 ≤ 𝐾2 < 0.07 𝐴5  

0.0009 ≤ 𝑚1 < 0.001 

0.001 ≤ 𝑚1 < 0.5 
𝐴1 

𝐴5 

𝑚1 = 0.001 

0.4 ≤ 𝑚2 < 1.2 𝐴5  

0.3 ≤ 𝛽1 < 0.99 

0.99 ≤ 𝛽1 < 2 
𝐴5 

𝐴1 

𝛽
1

= 0.99 

0.7 ≤ 𝛽2 < 2.5 𝐴5  

0.0001 ≤ 𝑙1 < 0.3 𝐴5  

0.1 ≤ 𝑙2 < 0.7 𝐴5  

0.0001 ≤ 𝑎1 < 1.6 𝐴5  

0.008 ≤ 𝑎2 < 2 𝐴5  

0.0006 ≤ 𝑎3 < 1.2 𝐴5  

0.0002 ≤ 𝑎4 < 3 𝐴5  

0.001 ≤ 𝑑1 < 0.475 

0.475 ≤ 𝑑1 < 0.799 

0.799 ≤ 𝑑1 < 1 

𝐴5 

𝐴3 

𝐴1 

𝑑1 = 0.475 

𝑑1 = 0.799 

0.01 ≤ 𝑑2 < 0.869 

0.869 ≤ 𝑑2 < 1 
𝐴5 

𝐴3 

𝑑2 = 0.869 

 

5. Conclusion  

     In this work, the occurrence of local bifurcation have been discussed with an appropriate 

conditions of food chain which contains a prey–intermediate predator-top predator model with 

fear and Crowley-Martin-type of functional response have been studied, transcritical and pitch 

fork bifurcation occurrence near 𝐴1. 𝐴2. 𝐴3 𝑎𝑛𝑑 𝐴4, while a saddle-node bifurcation occurs 

0 100 200 300 400 500 600 700 800 900 1000

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

0.8

Time

P
o

p
u

la
ti

o
n

a

 

 

prey

intermediate predator

top predator

0 100 200 300 400 500 600 700 800 900 1000
0

0.2

0.4

0.6

0.8

1

1.2

1.4
b

Time

P
o

p
u

la
ti

o
n

 

 

prey

intermediate predator

top predator



IHJPAS. 2025, 38(2) 

328 
 

near 𝐴5.Finally, in order to demonstrate how local bifurcation occurs in this system, numerical 

simulations are employed. The results of these simulations are as follows: 

*System (1) has no periodic dynamics 

*The parameters 𝑟𝑖. 𝑑𝑖 . 𝑖 = 1. 2 , 𝐾1. 𝑚1.  𝛽1 have a significant impact in the system dynamics 

(1), as opposed to other parameters  𝐾2. 𝑚2.  𝛽2. 𝑙𝑖. 𝑖 = 1. 2 𝑎𝑛𝑑 𝑎𝑖 . 𝑖 = 1. 2. 3. 4 the solution 

continues to approach the equilibrium point. 

 

Acknowledgment 

     Our researcher extends his Sincere thanks to the editor and members of the preparatory 

committee of the Ibn AL-Haitham Journal of Pure and Applied Science. 

 

Conflict of Interest  

     There are no conflicts of interest. 

 

Funding  

     There is no funding for the article. 

 

References 

1. Meng XY, Huo HF, Zhang XB. Stability and global Hopf bifurcation in a delayed food web 

consisting of a prey and two predators. Commun Nonlinear Sci Numer Simul. 2011;16:4335–4348. 

https://doi.org/10.1016/j.cnsns.2011.03.009 

2. Arancibia-Ibarra C, Aguirre P, Flores J, van Heijster P. Bifurcation analysis of a predator-prey model 

with predator intraspecific interactions and ratio-dependent functional response. Appl Math Comput. 

2021;402:126152. https://doi.org/10.1016/j.amc.2021.126152 

3. Magudeeswaran S, Vinoth S, Sathiyanathan K, Sivabalan M. Impact of fear on delayed three species 

food-web model with Holling Type-II functional response. Int J Biomath. 2022;15:2250014. 

https://doi.org/10.1142/S1793524522500140 

4. Majeed AA, Ismaeeb MH. The bifurcation analysis of prey-predator model in the presence of stage 

structured with harvesting and toxicity. In: Journal of Physics: Conference Series; IOP Publishing; 

2019; Vol. 1362, p. 12155. https://doi.org/10.1088/1742-6596/1362/1/012155 

5. Din Q. Stability, bifurcation analysis and chaos control for a predator-prey system. J Vib Control. 

2019;25:612–626. https://doi.org/10.1177/1077546318790871 

6. Prasad KD, Sasmal SK. Dynamics of anti-predator behavior and effect of fear on prey–predator 

model. J Biol Syst. 2022;30:887–912. https://doi.org/10.1142/S0218339022500322 

7. Shuai Z, Peng Y, Liu X, Li Z, Guerrero JM, Shen ZJ. Parameter stability region analysis of islanded 

microgrid based on bifurcation theory. IEEE Trans Smart Grid. 2019;10:6580–6591. 

8. Huang C, Li H, Cao J. A novel strategy of bifurcation control for a delayed fractional predator–prey 

model. Appl Math Comput. 2019;347:808–838. https://doi.org/10.1016/j.amc.2018.11.031 

9.    Pamuk S, İrem ÇAY. Stability and Hopf bifurcation analysis of a mathematical model in tumor 

angiogenesis. Anadolu Univ J Sci Technol A-Applied Sci Eng. 2018;19:50–57. 

10. Mukherjee D, Maji C. Bifurcation analysis of a Holling Type II predator-prey model with refuge. 

Chin J Phys. 2020;65:153–162. https://doi.org/10.1016/j.cjph.2020.02.012 

11. Majeed AA, Alabacy ZK. The persistence and bifurcation analysis of an ecological model with fear 

effect involving prey refuge and harvesting. In: AIP Conference Proceedings; AIP Publishing LLC; 

2022; Vol. 2394, p. 70001. https://doi.org/10.1063/5.0121877 

12. Zhang H, Cai Y, Fu S, Wang W. Impact of the fear effect in a prey-predator model incorporating a 

prey refuge. Appl Math Comput. 2019;356:328–337. https://doi.org/10.1016/j.amc.2019.03.034 

13. Wang Y, Cao J, Huang C. Dynamics of a fractional three-species food chain system with mixed 

functional responses and fear effect. 2023. https://doi.org/10.21203/rs.3.rs-2466067/v1 

https://doi.org/10.1016/j.cnsns.2011.03.009
https://doi.org/10.1016/j.amc.2021.126152
https://doi.org/10.1142/S1793524522500140
https://doi.org/10.1088/1742-6596/1362/1/012155
https://doi.org/10.1177/1077546318790871
https://doi.org/10.1142/S0218339022500322
https://doi.org/10.1016/j.amc.2018.11.031
https://doi.org/10.1016/j.cjph.2020.02.012
https://doi.org/10.1063/5.0121877
https://doi.org/10.1016/j.amc.2019.03.034


IHJPAS. 2025, 38(2) 

329 
 

14. Othman KB, Amen AI. Periodic solutions of the forest pest system via Hopf bifurcation and 

averaging theory. Iraqi J Sci. 2022;5496–5509. https://doi.org/10.24996/ijs.2022.63.12.35 

15. Aziz MM, Mohammed SAU. Analysis of stability and chaos of discrete time system with local 

bifurcation. In: 2022 8th International Conference on Contemporary Information Technology and 

Mathematics (ICCITM); IEEE; 2022; p. 425–429. 

16. Majeed NS. Local bifurcation analysis for a special type of SIR epidemic model. Int J Sci Res IJSR. 

2017;6:1143–1146. 

17. Majeed SN. Dynamical study of an SIR epidemic model with nonlinear incidence rate and regress 

of treatment. Ibn AL-Haitham Journal for Pure and Applied Sciences. 2018;510–522. 

https://doi.org/10.30526/2017.IHSCICONF.1810 

18. Jalal AA, Amen AI, Sulaiman NA. Darboux integrability of a generalized 3D chaotic Sprott ET9 

system. Baghdad Sci J. 2022;19:542. https://doi.org/10.21123/bsj.2022.19.3.0542 

19. Akhtar S, Ahmed R, Batool M, Shah NA, Chung JD. Stability, bifurcation and chaos control of a 

discretized Leslie prey-predator model. Chaos Solitons Fractals.2021;152:111345. 

https://doi.org/10.1016/j.chaos.2021.111345 

20. Aziz AA, Majeed AA. The fear effects on the dynamical study of an ecological model with Crowley-

Martin-Type of functional-response. In: AIP Conference Proceedings; AIP Publishing; 2024; Vol. 

3097. https://doi.org/10.1063/5.0209494 

 

https://doi.org/10.24996/ijs.2022.63.12.35
https://doi.org/10.30526/2017.IHSCICONF.1810
https://doi.org/10.21123/bsj.2022.19.3.0542
https://doi.org/10.1016/j.chaos.2021.111345
https://doi.org/10.1063/5.0209494