428 

© 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 web Prey-Predator Model involving fear and  

anti-Predator behavior 

 

Hanna Rasool Hadi1*  and Azhar Abbas Majeed2  
1,2Department of Mathematics, College of Science, University of Baghdad, City Baghdad, Country Iraq. 

*Corresponding Author. 

  

 Received:12 June 2023                             Accepted:27 August 2023                 Published:20 January 2025  

 doi.org/10.30526/38.1.3598 

 
 Abstract                                                            

     In this paper, the conditions under which the occurrence of the local bifurcation (such as saddle-

node (SN), transcritical (TC), and pitchfork (PT)) of all stable points of a food web model have been 

investigated. Fear and anti-predator responses involving the Holling-type IV and Growly-Martin 

functional responses have been found. It has been shown that there are transcritical and pitchfork 

bifurcations near 𝐻3, 𝐻4𝑎𝑛𝑑 𝐻6, as well as a saddle-node bifurcation close to the positive equilibrium 

point. In addition, there is a saddle-node bifurcation in close to the positive equilibrium point. These 

divergences have materialised into existence. In conclusion, to prove that the analytical results are 

correct, a numerical simulation of a set of parameters and starting conditions has been used. 

Keywords: Prey- Predator, Local bifurcation, Global bifurcation, Sotomayor’s theorem. 

 

1. Introduction  

       In the past few decades, a prominent topic in population dynamics has been the study of 

dynamical behaviors for massive predator-prey models. Numerous outcomes have been reported.  

Any biological or environmental characteristics are liable to change over time in the actual world. 

Since the selection pressures on systems in a fluctuating environment  differ from those in a stable 

environment, the impacts of a periodically altering environment are crucial for evolutionary 

theory(1,2). In the meantime, the gestational time delay is a frequent instance since, typically, the 

predator's present birth rate is determined by the amount of prey it has consumed in the past (3–5). 

The bifurcation analysis has been proposed to predict the parameter stability boundary of a nonlinear 

system. Bifurcation occurs when a small continuous change of parameter values causes a sudden 

change in system behavior  Compared with the small-signal analysis that analyzes the perturbation 

fixed equilibrium (6,7), the bifurcation analysis can perform parametric stability analysis that traces 

equilibrium solutions as the parameters change(8–10). 

Bifurcation theory analyses quantitative changes in the phase portrait, such as the emergence and 

disappearance of equilibria, periodic orbits, or more complex phenomena like weird attractors. 

Understanding of nonlinear dynamical systems depends critically on the techniques and findings of 

bifurcation theory(11–13). 

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


IHJPAS. 2025, 38 (1) 
 

429 
 

Bifurcation can take many distinct forms. (saddle-node, transcritical, pitchfork, period-doubling, 

Neimark-sacker, and hopf bifurcations) but the most significant ones are (SN), (TC), (PT), and hopf 

bifurcations. If the differential equation x=f(x,𝜇) had taken into account the point at which this 

dynamic system's behavior qualitatively changed, or what is described mathematically as a 

bifurcation(14–16). 

Many researchers have been studying bifurcation in dynamical systems accurately in the ordinary 

differential equations in both linear and nonlinear, autonomous and non-autonomous(17,18). 

 Orrell and Smith(19), studied visualizing bifurcation in high dimensional systems. Yuan and Zhao 

(20) studied  bifurcation for a fractional order  predator–prey system involving two nonidentical . 

While Majeed and Ismaeeb (21) studied the bifurcation analysis of a prey-predator model in the 

presence of a stage structured with harvesting and toxicity 

In this paper, the conditions of the local bifurcation (LB) of the food web model with fear and anti-

predator model have been found. 

 

2. Model Formulation (22): 

     The following section includes an ecological model with three species that had been proposed: the 

first and second preys with just one predator, which are indicated by the magnitude of their 

populations at the time. 𝐿1(t), L2(t)and L3(t) respectively . 

 𝑑𝐿1

𝑑𝑡
=  𝑆1𝐿1 (1 −

𝐿1

𝐾1
) −

𝐶1𝐿1𝐿3

𝑎1 + 𝐿1
2  

𝑑𝐿2

𝑑𝑡
=

𝑆2  𝐿2

1 + 𝑓1  𝐿3
− 𝐶𝐿2

2 −
𝐶2𝐿2𝐿3

(1 + 𝑎2𝐿2)(1 + 𝑎3𝐿3)
                                                                                    (1) 

𝑑𝐿3

𝑑𝑡
=

𝑔1𝐿1𝐿3

𝑎1 + 𝐿1
2 +

𝑔2𝐿2𝐿3

(1 + 𝑎2𝐿2)(1 + 𝑎3𝐿3)
− 𝜇𝐿1𝐿3 − 𝑑𝐿3 

The following table can be used to describe the positive system (1) parameters. 

Table 1: The system's (1) parameters are as follows: 

parameters Biological meaning 

𝑆𝑖 . 𝑖 = 1.2 First and second preys  growth rates, respectively 

𝐾1 > 0 Carrying capacity of first prey 

𝐶𝑖 > 0, 𝑖 = 1,2 Attack rate of the first prey and second prey by predator respecively 

𝑎1 > 0 Protection rate to the prey and predator by the environment 

C The rate of internal competition for the second prey 

𝑓1 The second prey's fear rate from predator 

𝑔𝑖 , 𝑖 = 1,2 Conversion rate of food from first prey and second prey to predator 

respecively 

𝑎2 Handling time of the second prey  

𝑎3 The magnitude of interference among predators 

𝜇 The rate of anti-predator behavior of first prey to predator 

d The natural death rate of predator 

 

 

 
 



IHJPAS. 2025, 38 (1) 
 

430 
 

3. Local bifurcation analysis: 

        The analysis of the (LB) of model (1) has been explored in this section, with a particular 

emphasis on the changes that occur around each equilibrium point when the parameter values in the 

dynamic behavior change. With the support of Sotomayor's theorem, our goal is to give higher order 

conditions that guarantee the appearance of the most frequent local bifurcations. 

Now, according to Jacobean matrix  𝐽(𝐿1, 𝐿2, 𝐿3) of the system (1) which is given in (22) as follows: 

𝐽 = [𝑎𝑖𝑗]
3×3

                                                                                                                                                         (2) 

where: 

𝑎11 = 𝑆1 −
2𝑆1𝐿1

𝐾1
−

𝐶1𝐿3(𝑎1−𝐿1
2)

(𝑎1+𝐿1
2)

2 , 𝑎12 = 0 ,  𝑎13 = −
𝐶1𝐿1

𝑎1+𝐿1
2 ,  

𝑎21 = 0 , 𝑎22 =
𝑆2

1 + 𝑓1𝐿3
− 2𝐶𝐿2 −

𝐶2𝐿3

(1 + 𝑎2𝐿2)(1 + 𝑎3𝐿3)
  

𝑎23 = −
𝑆2𝑓1𝐿2

(1+𝑓1𝐿3)2 −
𝐶2𝐿3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2  

𝑎31 =
𝑔1𝐿3(𝑎1−𝐿1

2)

(𝑎1+𝐿1
2)

2 − 𝜇𝐿3. 𝑎32 =
𝑔2𝐿3

(1+𝑎2𝐿2)2  (1+𝑎3𝐿3)  
 .  

 𝑎33 =
𝑔1𝐿1

𝑎1+𝐿1
2 +

𝑔2𝐿2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2 − 𝜇𝐿1 − 𝑑. 

For any non- zero vector V =(𝑣1, 𝑣2, 𝑣3)𝑇: 

𝐷2𝐹𝜇(𝑋, 𝜇)(𝑉, 𝑉) = [𝐴𝑖1]3×1,                                                                                                              (3) 

𝐴11 = −
2𝑆1

𝐾1
+

4𝑎1𝐶1𝐿1𝐿3

(𝑎1 + 𝐿1
2)2

𝑣1
2 −

𝐶1

(𝑎1 + 𝐿1
2)

((𝑎1 + 𝐿1
2) − 2𝐿1

2 + (𝑎1 + 2𝐿1))𝑣1𝑣3, 

𝐴21 = −2𝐶 +
𝐶2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)
𝑣1

2 − (
2𝑆2𝑓1

(1+𝑓1𝐿3)2 −
𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +
𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2) 𝑣2𝑣3 +

(
2𝑆2𝑓1

2𝐿3

(1+𝑓1𝐿3)3 −
𝐶2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3) 𝑣3
2, 

𝐴31 =
2𝑔1𝐿1𝐿3

𝑎1+𝐿1
2 𝑣1

2 + (
𝑔1(1−2𝐿1

2)

𝑎1+𝐿1
2 + 𝑔1 − 2𝜇) 𝑣1𝑣3 −

2𝑎2𝑔2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
𝑣2

2 −
𝑔2

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

1)𝑣2𝑣3 −
2𝑔2𝑎3𝐿3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3 𝑣3
2. 

and 

𝐷3𝐹𝜇(𝑋, 𝜇)(𝑉, 𝑉, 𝑉) = [𝐵𝑖1]3×1 ,                                                                                                         (4) 

𝐵11 =
4𝐶1𝑎1𝐿3(1−4𝐿1

2)

(𝑎1+𝐿1
2)

2 𝑣1
3 −

1

𝑎1+𝐿1
2 (2𝐿1(𝐶1 − 2) − 2𝐶1𝐿1((𝑎1 + 𝐿1

2) − 2𝐿1
2) + 2𝐶1(1 − 𝐿1(𝑎1 +

2𝐿1)) +
4𝐶1𝑎1𝐿1

𝑎1+𝐿1
2 ) 𝑣1

2𝑣3, 



IHJPAS. 2025, 38 (1) 
 

431 
 

𝐵21 = −
2𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
𝑣1

2𝑣2 + (
2𝐶2𝑎2

2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2
+

𝐶2𝑎2

(1+𝑎2𝐿2)
) 𝑣2

2𝑣3 +
𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3
𝑣2𝑣3

2 +

𝐶2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 𝑣1
2𝑣3 + (

4𝑆2𝑓1𝐿3

(1+𝑓1𝐿3)3 −
𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +
2𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)
) 𝑣2𝑣3

2 +

(
2𝑆2𝑓1

2(1−2𝑓1𝐿3)

(1+𝑓1𝐿3)6 +
4𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3) 𝑣3
3, 

𝐵31 =
2𝑔1𝐿3(𝑎1−𝐿1

2)

(𝑎1+𝐿1
2)

2 𝑣1
3 −

2𝑔1𝐿1(2𝑎1+1)

(𝑎1+𝐿1
2)

2 𝑣1
2𝑣3 +

6𝑎2
2𝑔2𝐿3

(1+𝑎2𝐿2)4(1+𝑎3𝐿3)
𝑣2

3 +
2𝑎2𝑔2

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

2)𝑣2
2𝑣3 +

2𝑔1𝐿1

𝑎1+𝐿1
2 𝑣1

2𝑣3 +
𝑎3𝑔2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3
(2𝑎2𝐿3 + 𝑎2𝐿2 − 2)𝑣2𝑣3

2 −
2𝑔2𝑎3(1−2𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)4 𝑣3
3. 

where 𝑋 = (𝐿1, 𝐿2, 𝐿3) and 𝜇 be any parameter. 

Theorem (1):System (1), with the value of the parameter  

 �́� = 𝑑 = 𝜇𝐾1 −
𝑔1𝑘1

𝑎1+𝐾1
2 −

𝑔2𝑆2

𝐶+𝑎2𝑆2
 , has a transcritical and pitchfork bifurcation at 𝐻3 = (𝐾1,

𝑆2

𝐶
, 0, ) if 

the following conditions: 

𝜇𝐾1 >
𝑔1𝑘1

𝑎1 + 𝐾1
2 +

𝑔2𝑆2

𝐶 + 𝑎2𝑆2
 ,                                                                                                                          (5) 

1

2
> 𝐿1

2                                                                                                                                                                   (6) 

�́�1 ≠ �́�2,                                                                                                                                           (7) 

�́�3 ≠ �́�4.                                                                                                                                           (8) 

where: 

�́�1 =
2𝑔1𝑀1

2𝐿1𝐿3

𝑎1+𝐿1
2 +

𝑔1𝑀1(1−2𝐿1
2)

𝑎1+𝐿1
2 + 𝑔1𝑀1 +

𝑔2𝑀2

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

�́�2 =
2𝑎2𝑔2𝑀2

2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
+

𝑎2𝐿2𝑔2𝑀2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 + 2𝜇𝑀1 +
2𝑔2𝑎3𝐿3

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

�́�3 =
2𝑔1𝑀1

3𝐿3(𝑎1−𝐿1
2)

(𝑎1+𝐿1
2)

2 +
6𝑎2

2𝑔2𝑀1
3𝐿3

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

2𝑔1𝑀1
2𝐿1

𝑎1+𝐿1
2 +

2𝑎2
2𝑔2𝑀2

2𝐿2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +

𝑎3𝑔2𝑀2

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

�́�4 =
2𝑔1𝑀1

2𝐿1(2𝑎1+1)

(𝑎1+𝐿1
2)

2 +
4𝑎2𝑔2𝑀2

2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
2𝑎3𝑔2𝑀2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +
2𝑔2𝑎3(1−2𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)4. 

Proof: Using the Jacobian matrix in eq. (1.8) in (22) 

  𝐽3 = 𝐽3(𝐻3, 𝑑) = [�́�𝑖𝑗]
3×3

, where �́�𝑖𝑗 = 𝑟𝑖𝑗 , except �́�33 = 0 .  

The characterizing eq. of �́�3 then has an eigenvalue of zero(say  λ3L3
). 

Now, let 𝑉1 = (𝑣1
[1]

, 𝑣2
[1]

, 𝑣3
[1]

)𝑇be the eigenvector corresponding to the eigenvalue  λ3L3
= 0. 

Thus,(𝐽3 −  λ3L3
𝐼)𝑉1 = 0, this gives: 

𝑣1
[1]

= 𝑀1𝑣3
[1]

, 𝑣2
[1]

= 𝑀2 𝑣3
[1]

 , where 𝑀1 = −
𝐶1𝐾1

(𝑎1+𝐾1
2)𝑆1

 , 𝑀2 = −
𝑆2𝑓1

2𝐶
. 

and 𝑣3
[1]

 any non-zero real number . 



IHJPAS. 2025, 38 (1) 
 

432 
 

 Consider 𝑍[1] = (𝑧1
[1]

, 𝑧2
[1]

, 𝑧3
[1]

)𝑇 be the eigenvector connected with an eigenvalue  λ3L3
= 0 of the 

matrix 𝐽3
 𝑇.Then, (𝐽3

 𝑇 −   λ3L3
𝐼) 𝑍[1] = 0. 

By solving this equation for , 𝑍[1] = (0,0, 𝑧3
[1]

)
𝑇

, where 𝑍3
[1]

 any non-zero real number. 

Now, consider that: 

𝜕𝑓

𝜕𝑑
= 𝑓𝑑(𝑋, 𝑑) = (

𝜕𝑓1

𝜕𝑑
,
𝜕𝑓2

𝜕𝑑
,
𝜕𝑓3

𝜕𝑑
)

𝑇

= (0, 0, −𝐿3)𝑇 . 

So, 𝑓𝑑(𝐻3, �́�) = (0,0,0)𝑇 and hence (𝑍[1])
𝑇

𝑓𝑑(𝐻3, �́�) = 0 

The (SD) bifurcation requirement cannot be satisfied according to Sotomayor's theorem. As a result, 

the first requirement for (TC) bifurcation is realized. Now 

𝐷𝑓𝑑(𝑋, 𝑑) = [
0 0 0
0 0 0
0 0 −1

] 

where, 𝐷𝑓𝑑(𝑋, 𝑑) represents the derivative of 𝑓𝑑(𝑋, 𝑑) with respect to 𝑋 = (𝐿1, 𝐿2, 𝐿3)𝑇 . 
Furthermore, it is observed that: 

𝐷𝑓𝑑(𝐻3, 𝑑)𝑉[1] = [
0 0 0
0 0 0
0 0 −1

] [

𝑀1𝑣3
[1]

𝑀2 𝑣3
[1]

𝑣3
[1]

] = [

0
0

−𝑣3
[1]

] 

(𝑍[1])
𝑇

[𝐷𝑓𝑑(𝐻3, 𝑑)𝑉[1]] = (0,0, 𝑧3
[1]

)
𝑇

(0,0, −𝑣3
[1]

) 

                      = −𝑣3
[1]

𝑧3
[1]

≠ 0. 

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

𝐷2𝐹𝜇(𝐻3, 𝑑)(𝑉[1], 𝑉[1]) = (𝐴𝑖1
́ )

3×1
, 

�́�11 = −
2𝑆1

𝐾1
+ (𝑣3

[1]
)

2

[
4𝑎1𝐶1𝐿1𝐿3

(𝑎1+𝐿1
2)

2 𝑀1
2 −

𝐶1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) − 2𝐿1

2 + (𝑎1 + 2𝐿1))𝑀1], 

�́�21 = −2𝐶 + (𝑣3
[1]

)
2

[
𝐶2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)
𝑀1

2 − (
2𝑆2𝑓1

(1+𝑓1𝐿3)2 −
𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2) 𝑀2 + (
2𝑆2𝑓1

2𝐿3

(1+𝑓1𝐿3)3 −
𝐶2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3)], 

�́�31 = (𝑣3
[1]

)
2

[
2𝑔1𝑀1

2𝐿1𝐿3

𝑎1+𝐿1
2 + (

𝑔1(1−2𝐿1
2)

𝑎1+𝐿1
2 + 𝑔1 − 2𝜇) 𝑀1 −

2𝑎2𝑔2𝑀2
2𝐿3

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

𝑎2𝐿2𝑔2𝑀2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

𝑔2𝑀2

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

2𝑔2𝑎3𝐿3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3
].  

Hence, it was obtained 

(𝑍[1])
𝑇

[𝐷2𝐹𝜇(𝐻3, 𝑑)(𝑉[1], 𝑉[1])] = 𝑧3
[1]

(𝑣3
[1]

)
2

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

So, by condition (7) that system (1) has a (TC) bifurcation at 𝐻3 with a parameter �́� = 𝑑 ,If the 

requirement (7) is not met. 



IHJPAS. 2025, 38 (1) 
 

433 
 

By substituting 𝑉[1] into (4) we get: 

𝐷3𝐹𝜇(𝐻3, 𝑑)(𝑉[1], 𝑉[1], 𝑉[1]) = (�́�𝑖1)
3×1

 

�́�11 = (𝑣3
[1]

)
3

[
4𝐶1𝑎1𝐿3(1−4𝐿1

2)

(𝑎1+𝐿1
2)

2 𝑀1
3 −

𝑀1
2

𝑎1+𝐿1
2 (2𝐿1(𝐶1 − 2) − 2𝐶1𝐿1((𝑎1 + 𝐿1

2) − 2𝐿1
2) + 2𝐶1(1 −

𝐿1(𝑎1 + 2𝐿1)) +
4𝐶1𝑎1𝐿1

𝑎1+𝐿1
2 )], 

�́�21 = (𝑣3
[1]

)
3

[−
2𝑀1

2𝑀2𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
+

2𝐶2𝑀2
2𝑎2

2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
𝑀2

2𝐶2𝑎2

(1+𝑎2𝐿2)
+

𝑀2𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +

𝐶2𝑀1
2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2
+ (

4𝑆2𝑀2𝑓1𝐿3

(1+𝑓1𝐿3)3
−

𝑀2𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3
+

2𝑀2𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)
) + (

2𝑆2𝑓1
2(1−2𝑓1𝐿3)

(1+𝑓1𝐿3)6
+

4𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3
)], 

�́�31 = (𝑣3
[1]

)
3

[
2𝑔1𝑀1

3𝐿3(𝑎1−𝐿1
2)

(𝑎1+𝐿1
2)

2 −
2𝑔1𝑀1

2𝐿1(2𝑎1+1)

(𝑎1+𝐿1
2)

2 +
6𝑎2

2𝑔2𝑀1
3𝐿3

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

2𝑎2
2𝑔2𝑀2

2𝐿2

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

4𝑎2𝑔2𝑀2
2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
2𝑔1𝑀1

2𝐿1

𝑎1+𝐿1
2 +

𝑎3𝑔2𝑀2

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

2𝑎3𝑔2𝑀2

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

2𝑔2𝑎3(1−2𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)4. 

(𝑍[1])
𝑇

[𝐷3𝐹𝜇(𝐻3, 𝑑)(𝑉[1], 𝑉[1], 𝑉[1])] = 𝑧3
[1]

(𝑣3
[1]

)
3

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

So, there is pitchfork bifurcation at 𝐻3 where�́� = 𝑑. 

Theorem (2): Suppose that condition (1.9b), revers conditions (1.9c) and (1.9d) in (22) with the 

following conditions are satisfied: 

𝐿1
2 <

1

4 
,                                                                                                                                                               (8) 

𝐿3 < 𝑚𝑖𝑛 {
1

𝑎3
,

1

2𝑓1
},                                                                                                                                         (9) 

�̇�31�̇�13 = �̇�32�̇�23                                                                                                                                               (10) 

�̇�1 ≠ �̇�2,                                                                                                                                                            (11) 

�̇�3 ≠ �̇�4,                                                                                                                                            (12)  

Where: 

�̇�1 = (𝑣3
[2]

)
2

(
4𝑎1𝐶1�̇�3𝛼1̇

2𝐿1𝐿3

(𝑎1+𝐿1
2)

2 +
2𝐶1𝛼1�̇�3̇ 𝐿1

2

(𝑎1+𝐿1
2)

+
�̇�1

2�̇�4𝐶2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)
+

�̇�2�̇�4𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +
2�̇�4𝑆2𝑓1

2𝐿3

(1+𝑓1𝐿3)3 +

2𝑔1�̇�1
2�̇�4𝐿1𝐿3

𝑎1+𝐿1
2 + (

𝑔1(1−2𝐿1
2)

𝑎1+𝐿1
2 + 𝑔1) �̇�1 +

𝑔2�̇�2

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

�̇�2 =
2𝑆1�̇�3

𝐾1
+ 2𝐶�̇�4 + (𝑣3

[2]
)

2
(

𝐶1�̇�1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) + (𝑎1 + 2𝐿1)) + (

2𝑆2𝑓1

(1+𝑓1𝐿3)2 +

𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2
) �̇�2�̇�4 +

𝐶2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3
+ 2𝜇�̇�1 +

2𝑎2𝑔2�̇�2
2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
+

𝑎2𝐿2𝑔2�̇�2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2
+

2𝑔2𝑎3𝐿3

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



IHJPAS. 2025, 38 (1) 
 

434 
 

�̇�3 =
2𝐶1�̇�1

3�̇�3𝑎1𝐿3(1−4𝐿1
2)

(𝑎1+𝐿1
2)

2 + �̇�1
2�̇�3𝐶1𝐿1 +

�̇�1
2�̇�3𝐿1

𝑎1+𝐿1
2 +

�̇�1
2�̇�3𝐶1𝐿1(𝑎1+2𝐿1)

𝑎1+𝐿1
2 +

4�̇�1
2�̇�3𝐶1𝑎1𝐿1

(𝑎1+𝐿1
2)

2 +

2𝐶2�̇�2
2�̇�4𝑎2

2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2
+

�̇�2
2�̇�4𝐶2𝑎2

(1+𝑎2𝐿2)
+

�̇�2�̇�4𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3
+

𝐶2�̇�1
2�̇�4

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2
+

4𝑆2�̇�2�̇�4𝑓1𝐿3

(1+𝑓1𝐿3)3
+

2�̇�2�̇�4𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)
+

2𝑆2𝑓1
2(1−2𝑓1𝐿3)

(1+𝑓1𝐿3)6 +
4𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3+
2𝑔1�̇�1

3𝐿3(𝑎1−𝐿1
2)

(𝑎1+𝐿1
2)

2 +
6𝑎2

2𝑔2�̇�1
3𝐿3

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

2𝑎2
2𝑔2�̇�2

2𝐿2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2
+

2𝑔1�̇�1
2𝐿1

𝑎1+𝐿1
2 +

𝑎3𝑔2�̇�2

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

�̇�4 =
�̇�1

2�̇�3𝐶1𝐿1

𝑎1+𝐿1
2 +

2�̇�1
2�̇�3𝐶1𝐿1

3

𝑎1+𝐿1
2 +

�̇�1
2�̇�3𝐶1

𝑎1+𝐿1
2 +

2�̇�1
2�̇�3�̇�2𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
+

�̇�2𝐶2𝑎2�̇�4(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +
2𝑔1�̇�1

2𝐿1(2𝑎1+1)

(𝑎1+𝐿1
2)

2 +

4𝑎2𝑔2�̇�2
2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
2𝑎3𝑔2�̇�2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +
2𝑔2𝑎3(1−2𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)4. 

Then, system (1) with parameter value : 

 �̇� = 𝑑 = −
(𝑎1+�̇�1

2)(𝑒11𝑒23𝑒32+𝑒13𝑒31𝑒22)+𝑔1�̇�1𝑒11𝑒22−𝜇�̇�1𝑒11𝑒22(𝑎1+�̇�1
2)

𝑒11𝑒22(𝑎1+�̇�1
2)

, has a transcritical and pitch fork 

bifurcation at 𝐻4, similary for 𝐻5 

Proof: Using the Jacobian matrix in eq. (1.9) in (22) 

 𝐽4̇ = 𝐽4(𝐻4, �̇�) = [�̇�𝑖𝑗]
3×3

, where �̇�𝑖𝑗 = 𝑒𝑖𝑗 , except 

 �̇�33 =
𝑔1�̇�1

(𝑎1+�̇�1
2)

− 𝜇�̇�1 − �̇� .  

the characterizing eq. of  𝐽4̇ then has an eigenvalue of zero (say λ4L3
 ). 

Now, let 𝑉[2] = (𝑣1
[2]

, 𝑣2
[2]

, 𝑣3
[2]

)𝑇be the eigenvector corresponding to the eigenvalue  λ4L3
= 0. 

Thus,(𝐽4̇ −  λ4L3
𝐼)𝑉[2] = 0, this gives: 

𝑣1
[2]

= �̇�1𝑣3
[2]

, 𝑣2
[2]

= �̇�2𝑣3
[2]

    . 

where: �̇�1 = −
�̇�13

�̇�11
 , �̇�2 = −

�̇�23

�̇�22
 , and 𝑣3

[2]
 any non-zero real number . 

Consider  𝑍[2] = (𝑧1
[2]

, 𝑧2
[2]

, 𝑧3
[2]

)𝑇 be the eigenvector connected with an eigenvalue  λ4L3
= 0 of the 

matrix 𝐽4̇
 𝑇.Then, (𝐽4̇

 𝑇 −   λ4L3
𝐼) 𝑍[2] = 0. 

By solving this equation for , 𝑍[2] = (�̇�3𝑧3
[2]

, �̇�4𝑧3
[2]

, 𝑧3
[2]

)
𝑇

, where: 

�̇�3 = −
�̇�31

�̇�11
 , �̇�4 = −

𝑒32̇

𝑒22̇
 , and 𝑧3

[2]
 any non-zero real number. 

Now, consider that: 

𝜕𝑓

𝜕𝑑
= 𝑓𝑑(𝑋, 𝑑) = (

𝜕𝑓1

𝜕𝑑
,
𝜕𝑓2

𝜕𝑑
,
𝜕𝑓3

𝜕𝑑
)

𝑇

= (0,0, −𝐿3)𝑇 . 

So, 𝑓𝑑(𝐻4, 𝑑) = (0,0,0)𝑇 and hence (𝑍[2])
𝑇

𝑓𝑑(𝐻4, 𝑑) = 0. 

The (SD) bifurcation requirement cannot be satisfied according to Sotomayor's theorem. As a result, 

the first requirement for (TC) bifurcation is realized. Now 



IHJPAS. 2025, 38 (1) 
 

435 
 

𝐷𝑓𝑑(𝑋, 𝑑) = [
0 0 0
0 0 0
0 0 −1

], 

where, 𝐷𝑓𝑑(𝑋, 𝑑) represents the derivative of 𝑓𝑑(𝑋, 𝑑) with respect to 𝑋 = (𝐿1, 𝐿2, 𝐿3)𝑇 . 
Furthermore, it is observed that: 

𝐷𝑓𝑑(𝐻4, 𝑑)𝑉[2] = [
0 0 0
0 0 0
0 0 −1

] [

�̇�1𝑣3
[2]

�̇�2𝑣3
[2]

𝑣3
[2]

] = [

0
0

−𝑣3
[2]

] 

(𝑍[2])
𝑇

[𝐷𝑓𝑑(𝐻4, 𝑑)𝑉[2]] = (�̇�3𝑧3
[2]

, �̇�4𝑧3
[2]

, 𝑧3
[2]

)
𝑇

(0,0, −𝑣3
[2]

) 

                                              = −𝑣3
[2]

𝑧3
[2]

≠ 0 

𝐷2𝐹𝜇(𝐻4, 𝑑)(𝑉[2], 𝑉[2]) =  [𝐴𝑖1
̇ ]

3×1
 

�̇�11 = −
2𝑆1

𝐾1
+ (𝑣3

[2]
)

2

[
4𝑎1𝐶1𝛼1̇

2𝐿1𝐿3

(𝑎1+𝐿1
2)

2 −
𝐶1�̇�1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) + (𝑎1 + 2𝐿1)) +

2𝐶1𝛼1̇𝐿1
2

(𝑎1+𝐿1
2)

], 

�̇�21 = −2𝐶 + (𝑣3
[2]

)
2

[
�̇�1

2𝐶2𝐿3

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

2𝑆2𝑓1

(1+𝑓1𝐿3)2 +
𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2) �̇�2 +
�̇�2𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

(
2𝑆2𝑓1

2𝐿3

(1+𝑓1𝐿3)3 −
𝐶2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3)], 

�̇�31 = (𝑣3
[2]

)
2

[
2𝑔1�̇�1

2𝐿1𝐿3

𝑎1+𝐿1
2 + (

𝑔1(1−2𝐿1
2)

𝑎1+𝐿1
2 + 𝑔1) �̇�1 − 2𝜇�̇�1 −

2𝑎2𝑔2�̇�2
2𝐿3

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

𝑎2𝐿2𝑔2�̇�2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

𝑔2�̇�2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 −
2𝑔2𝑎3𝐿3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3]. 

Hence, it was obtained . 

(𝑍[2])
𝑇

[𝐷2𝐹𝜇(𝐻4, 𝑑)(𝑉[2], 𝑉[2])] = 𝑧3
[2]

(�̇�1 − �̇�2) ≠ 0, 

This implies that system (1) has a (TC) bifurcation at 𝐻4 with a parameter �̇� = 𝑑, If the requirement 

(11) unsatisfied then. 

By substituting 𝑉[2] in (4) we get: 

𝐷3𝐹𝜇(𝐻4, 𝑑)(𝑉[2], 𝑉[2], 𝑉[2]) = (𝐵𝑖1
̇ )

3×1
 

�̇�11 = 2 (𝑣3
[1]

)
3

[
2𝐶1�̇�1

3𝑎1𝐿3(1−4𝐿1
2)

(𝑎1+𝐿1
2)

2 −
�̇�1

2𝐶1𝐿1

𝑎1+𝐿1
2 + �̇�1

2𝐶1𝐿1 +
�̇�1

2𝐿1

𝑎1+𝐿1
2 −

2�̇�1
2𝐶1𝐿1

3

𝑎1+𝐿1
2 −

�̇�1
2𝐶1

𝑎1+𝐿1
2 +

�̇�1
2𝐶1𝐿1(𝑎1+2𝐿1)

𝑎1+𝐿1
2 +

4�̇�1
2𝐶1𝑎1𝐿1

(𝑎1+𝐿1
2)

2 ], 

�̇�21 = (𝑣3
[1]

)
3

[−
2�̇�1

2�̇�2𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
+

2𝐶2�̇�2
2𝑎2

2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
�̇�2

2𝐶2𝑎2

(1+𝑎2𝐿2)
+

�̇�2𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +

𝐶2�̇�1
2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 + (
4𝑆2�̇�2𝑓1𝐿3

(1+𝑓1𝐿3)3 −
�̇�2𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +
2�̇�2𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)
) + (

2𝑆2𝑓1
2(1−2𝑓1𝐿3)

(1+𝑓1𝐿3)6 +

4𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3)], 



IHJPAS. 2025, 38 (1) 
 

436 
 

�̇�31 = (𝑣3
[1]

)
3

[
2𝑔1�̇�1

3𝐿3(𝑎1−𝐿1
2)

(𝑎1+𝐿1
2)

2 −
2𝑔1�̇�1

2𝐿1(2𝑎1+1)

(𝑎1+𝐿1
2)

2 +
6𝑎2

2𝑔2�̇�1
3𝐿3

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

2𝑎2
2𝑔2�̇�2

2𝐿2

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

4𝑎2𝑔2�̇�2
2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2
+

2𝑔1�̇�1
2𝐿1

𝑎1+𝐿1
2 +

𝑎3𝑔2�̇�2

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

2𝑎3𝑔2�̇�2

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

2𝑔2𝑎3(1−2𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)4. 

(𝑍[2])
𝑇

[𝐷3𝐹𝜇(𝐻4, 𝑑)(𝑉[2], 𝑉[2], 𝑉[2])] = (𝑣3
[2]

)
3

𝑧3
[4]

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

So, Pitch fork bifurcation can be found at 𝐻4 where �̇� = 𝑑. 

Theorem (3):  system (1) with parameter value : 

 𝑆1̅ = 𝑆1 =
𝐶1�̅�3

𝑎1
,  has a transcritical and pitchfork  bifurcation at 𝐻6.if the next conditions hold: 

�̅�1 ≠ �̅�2                                                                                                                                                            (13) 

     �̅�3 ≠ �̅�4                                                                                                                                                                                                               (14)                                       

Where: 

�̅�1 = (𝑣1
[3]

)
2

[
4𝑎1𝐶1𝐿1𝐿3

(𝑎1+𝐿1
2)

2 +
2𝐶1𝐿1

2

(𝑎1+𝐿1
2)

], 

�̅�1 =
2𝑆1

𝐾1
+ (𝑣1

[3]
)

2 𝐶1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) + (𝑎1 + 2𝐿1)), 

�̅�3 =
2𝐶1𝜎1

3𝑎1𝐿3(1−4𝐿1
2)

(𝑎1+𝐿1
2)

2 + �̇�1
2𝐶1𝐿1 +

𝜎1
2𝐿1

𝑎1+𝐿1
2 +

𝜎1
2𝐶1𝐿1(𝑎1+2𝐿1)

𝑎1+𝐿1
2 +

4𝜎1
2𝐶1𝑎1𝐿1

(𝑎1+𝐿1
2)

2 , 

�̅�4=
𝜎1

2𝐶1𝐿1

𝑎1+𝐿1
2 +

2𝜎1
2𝐶1𝐿1

3

𝑎1+𝐿1
2 −

𝜎1
2𝐶1

𝑎1+𝐿1
2. 

Proof: Using  the Jacobian matrix is presented in equation (1.10) in (22), 

 𝐽6̅ = 𝐽6(𝐻6, 𝑆1̅) = [𝑡�̅�𝑗]
3×3

, where 𝑡�̅�𝑗 = 𝑡𝑖𝑗 , except 𝑡1̅1 = 0 .  

Then the characterizing equation of  𝐽6̅  had a zero eigenvalue (say λ6L1
 ).   

Now, let 𝑉[3] = (𝑣1
[3]

, 𝑣2
[3]

, 𝑣3
[3]

)𝑇be the eigenvector corresponding to the eigenvalue  λ6L1
= 0. 

Thus,(𝐽6̅ −  λ6L1
𝐼)𝑉[3] = 0, this gives: 

 𝑣2
[3]

= 𝜎1𝑣1
[3]

 , 𝑣3
[3]

= 𝜎2𝑣1
[3]

, 

Where: 

 𝜎1 = −
𝑡31𝑡23

𝑡32𝑡23+𝑡22𝑡33
, 𝜎2 =

𝑡31𝑡22

𝑡32𝑡23+𝑡22𝑡33
 and 𝑣1

[6]
 any non-zero real number.  

Let 𝑍[3] = (𝑧1
[3]

, 𝑧2
[3]

, 𝑧3
[3]

)𝑇 be the eigenvector connected with an eigenvalue  λ6L1
= 0 of the matrix 

𝐽6̅
 𝑇.Then, (𝐽6̅

 𝑇 −   λ6L1
𝐼) 𝑍[3] = 0. 

By solving this equation for , 𝑍[3] = (𝑧1
[3]

, 0,0)
𝑇

,  

Where 𝑧1
[3]

 any non-zero real number. 



IHJPAS. 2025, 38 (1) 
 

437 
 

Now, consider that: 

𝜕𝑓

𝜕𝑆1
= 𝑓𝑆1

(𝑋, 𝑆1) = (
𝜕𝑓1

𝜕𝑆1
,
𝜕𝑓2

𝜕𝑆1
,
𝜕𝑓3

𝜕𝑆1
)

𝑇

= (𝐿1(1 −
1

𝐾1
),0,0)

𝑇

. 

So, 𝑓𝑆1
(𝐻6, 𝑆1̅) = (0,0,0)𝑇 , and hence  (𝑍[3])

𝑇
𝑓𝑆1

(𝐻6, 𝑆1̅) = 0. 

The (SD) bifurcation requirement cannot be satisfied according to Sotomayor's theorem. As a result, 

the first requirement for (TC) bifurcation is realized. Now 

 

𝐷𝑓�̅�1
(𝑋, 𝑆1̅) = [

1 −
1

𝐾1
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
(𝐻6, 𝑆1̅)𝑉[3] = [

1 −
1

𝐾1
0 0

0 0 0
0 0 0

] [

𝑣1
[3]

𝜎1𝑣1
[3]

𝜎2𝑣1
[3]

] = [
(1 −

1

𝐾1
) 𝑣1

[3]

0
0

] 

(𝑍[3])
𝑇

[𝐷𝑓�̅�1
(𝐻6, 𝑆1̅)𝑉[3]] = (𝑧1

[3]
, 0,0)

𝑇

((1 −
1

𝐾1
) 𝑣1

[3]
, 0,0) 

                      = (1 −
1

𝐾1
) 𝑣1

[3]
𝑧1

[3]
≠ 0 

 

By substituting 𝑉[6] in (3) we get: 

𝐷2𝐹𝜇(𝐻6, 𝑆1)(𝑉[3], 𝑉[3]) = [�̅�𝑖1]3×1 , 

�̅�11 = −
2𝑆1

𝐾1
+ (𝑣1

[3]
)

2

[
4𝑎1𝐶1𝐿1𝐿3

(𝑎1+𝐿1
2)

2 −
𝐶1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) + (𝑎1 + 2𝐿1)) +

2𝐶1𝐿1
2

(𝑎1+𝐿1
2)

], 

�̅�21 = −2𝐶 + (𝑣3
[1]

)
2

[
𝐶2𝐿3

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

2𝑆2𝑓1

(1+𝑓1𝐿3)2 +
𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2) 𝜎1 +
𝜎1𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

(
2𝑆2𝑓1

2𝜎2
2𝐿3

(1+𝑓1𝐿3)3 −
𝜎2

2𝐶2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3)], 

�̅�31 = (𝑣1
[3]

)
2

[
2𝑔1𝐿1𝐿3

𝑎1+𝐿1
2 + (

𝑔1(1−2𝐿1
2)

𝑎1+𝐿1
2 + 𝑔1) 𝜎2 − 2𝜇𝜎1 −

2𝑎2𝑔2𝜎1𝐿3

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

𝑎2𝐿2𝑔2𝜎1𝜎2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

𝑔2𝜎1𝜎2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 −
2𝑔2𝑎3𝜎1𝐿3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3]. 

Hence, it was obtained 

(𝑍[3])
𝑇

[𝐷2𝐹𝜇(𝐻6, 𝑆1)(𝑉[3], 𝑉[3])] = 𝑧1
[3](�̅�1 − �̅�2) ≠ 0, 

So, system (1) has  (TC) bifurcation at 𝐻6 with a parameter  𝑆1̅ = 𝑆1. If the requirement (13) not 

satisfied then . 



IHJPAS. 2025, 38 (1) 
 

438 
 

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

𝐷3𝐹𝜇(𝐻6, 𝑆1)(𝑉[3], 𝑉[3], 𝑉[3]) = (�̅�𝑖1)3×1 

�̅�11 = 2 (𝑣1
[3]

)
3

[
2𝐶1𝜎1

3𝑎1𝐿3(1−4𝐿1
2)

(𝑎1+𝐿1
2)

2 −
𝜎1

2𝐶1𝐿1

𝑎1+𝐿1
2 + �̇�1

2𝐶1𝐿1 +
𝜎1

2𝐿1

𝑎1+𝐿1
2 −

2𝜎1
2𝐶1𝐿1

3

𝑎1+𝐿1
2 −

𝜎1
2𝐶1

𝑎1+𝐿1
2 +

𝜎1
2𝐶1𝐿1(𝑎1+2𝐿1)

𝑎1+𝐿1
2 +

4𝜎1
2𝐶1𝑎1𝐿1

(𝑎1+𝐿1
2)

2 ], 

�̅�21 = (𝑣1
[3]

)
3

[−
2�̇�1

2�̇�2𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
+

2𝐶2�̇�2
2𝑎2

2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
�̇�2

2𝐶2𝑎2

(1+𝑎2𝐿2)
+

�̇�2𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3 +

𝐶2�̇�1
2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2
+ (

4𝑆2�̇�2𝑓1𝐿3

(1+𝑓1𝐿3)3
−

�̇�2𝐶2𝑎2(1−𝑎3𝐿3)

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)3
+

2�̇�2𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)
) + (

2𝑆2𝑓1
2(1−2𝑓1𝐿3)

(1+𝑓1𝐿3)6
+

4𝐶2𝑎3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3)], 

�̅�31 = (𝑣1
[3]

)
3

[
2𝑔1�̇�1

3𝐿3(𝑎1−𝐿1
2)

(𝑎1+𝐿1
2)

2 −
2𝑔1�̇�1

2𝐿1(2𝑎1+1)

(𝑎1+𝐿1
2)

2 +
6𝑎2

2𝑔2�̇�1
3𝐿3

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

2𝑎2
2𝑔2�̇�2

2𝐿2

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

4𝑎2𝑔2�̇�2
2

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)2 +
2𝑔1�̇�1

2𝐿1

𝑎1+𝐿1
2 +

𝑎3𝑔2�̇�2

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

2𝑎3𝑔2�̇�2

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

2𝑔2𝑎3(1−2𝑎3𝐿3)

(1+𝑎2𝐿2)(1+𝑎3𝐿3)4. 

(𝑍[3])
𝑇

[𝐷3𝐹𝜇(𝐻6, 𝑆1)(𝑉[3], 𝑉[3], 𝑉[3])] = 2 (𝑣1
[3]

)
3

𝑧1
[3]

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

So, Pitchfork bifurcation can be found at 𝐻6 where 𝑆1̅ = 𝑆1. 

Theorem (4):   If conditions (1.11c) and (1.11d) in (22) are reversing if the following requirements 

are fulfilled: 

ℎ̃31ℎ̃31 = ℎ̃32ℎ̃23 .                                                                                                                                           (15)  

ℎ̃13 = ℎ̃22                                                                                                                                                          (16) 

�̃�1 > �̃�2                                                                                                                                                              (17) 

�̃�1 ≠ �̃�2,                                                                                                                                                            (18) 

where: 

�̃�1 =
4𝑎1𝐶1𝐿1𝐿3

(𝑎1+𝐿1
2)

2 �̃�1 +
𝐶2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)
�̃�1

2 +
2𝑆2𝑓1

2𝐿3

(1+𝑓1𝐿3)3
+

2𝑔1𝐿1𝐿3

𝑎1+𝐿1
2 �̃�1

2 + (
𝑔1(1−2𝐿1

2)

𝑎1+𝐿1
2 + 𝑔1 − 2𝜇) �̃�1, 

�̃�2 =
𝐶1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) − 2𝐿1

2 + (𝑎1 + 2𝐿1))�̃�1 + (
2𝑆2𝑓1

(1+𝑓1𝐿3)2 +
𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2) �̃�1 +
2𝑎2𝑔2𝐿3

(1+𝑎2𝐿2)3(1+𝑎3𝐿3)
�̃�2

2 +
𝑔2

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2
(𝑎2𝐿2 − 1)�̃�2 +

2𝑔2𝑎3𝐿3

(1+𝑎2𝐿2)(1+𝑎3𝐿3)3. 

Then, system (1) with parameter value:  

   �̃� = 𝑑 =
�̃�1−�̃�2

(ℎ11ℎ22)(𝑎1+�̃�1
2)(1+𝑎2�̃�2)(1+𝑎3�̃�3)2  , where: 

�̃�1 = −(𝑎1 + �̃�1
2)(1 + 𝑎2�̃�2)(1 + 𝑎3�̃�3)

2
(ℎ11ℎ23ℎ32 + ℎ13ℎ31ℎ22) + (ℎ11ℎ22)(1 + 𝑎2�̃�2)(1 +

𝑎3�̃�3)
2

𝑔1�̃�1 + (ℎ11ℎ22)(𝑎1 + �̃�1
2)𝑔2�̃�2, 

�̃�2 = −(ℎ11ℎ22)(𝑎1 + �̃�1
2)(1 + 𝑎2�̃�2)(1 + 𝑎3�̃�3)2𝜇�̃�1. 



IHJPAS. 2025, 38 (1) 
 

439 
 

has a saddle-node bifurcation at 𝐻7 = (�̃�1, �̃�2, �̃�3). 

Proof: Using the Jacobian matrix given by eq. (1.11) in (22)  

 𝐽7 = 𝐽7(𝐻7, �̃�) = [ℎ̃𝑖𝑗]
3×3

, where ℎ̃𝑖𝑗 = ℎ𝑖𝑗  ,except 

 ℎ̃33 =
𝑔1�̃�1

𝑎1 + �̃�1
2

+
𝑔2�̃�2

(1 + 𝑎2�̃�2)(1 + 𝑎3�̃�3)2
− 𝜇�̃�1 − 𝑑 

 the characterizing eq. 𝐽7 had a zero eigenvalue (say λ7L3
= 0) if and only if �̿�3 = 0 . 

Now, let 𝑉[7] = (𝑣1
[7]

, 𝑣2
[7]

, 𝑣3
[7]

)𝑇be the eigenvector corresponding to the eigenvalue  λ7L3
= 0. 

Thus,(𝐽7 −  λ7L3
𝐼)𝑉[7] = 0, this gives: 

𝑣1
[7]

= �̃�1𝑣3
[7]

, 𝑣2
[7]

= �̃�2𝑣3
[7]

  , 𝑣3
[7]

 

and 𝑣3
[7]

 any non-zero real number  

Where: 

�̃�1= −
ℎ13

ℎ11
,  �̃�2 = −

ℎ23

ℎ22
. 

Let 𝑍[7] = (𝑧1
[7]

, 𝑧2
[7]

, 𝑧3
[7]

)𝑇 be the eigenvector connected with to the eigenvalue  λ7L3
= 0 of the 

matrix 𝐽7
 𝑇.Then, (𝐽7

 𝑇 −   λ7L3
𝐼) 𝑍[7] = 0. 

By solving this equation for , 𝑍[7] = (�̃�3𝑧3
[7]

, �̃�4𝑧3
[7]

, 𝑧3
[7]

)
𝑇

, 

and 𝑧3
[7]

 any non-zero real number. 

Where: 

�̃�3 = −
ℎ31̃

ℎ11̃
, �̃�4 = −

ℎ̃32

ℎ22̃
. 

Now, consider that: 

𝜕𝑓

𝜕𝑑
= 𝑓𝑑(𝑋, 𝑑) = (

𝜕𝑓1

𝜕𝑑
,
𝜕𝑓2

𝜕𝑑
,
𝜕𝑓3

𝜕𝑑
)

𝑇

= (0,0, −𝐿3)𝑇 . 

So, 𝑓𝑑(𝐻7, 𝑑) = (0,0, −�̃�3)
𝑇

and hence  

(𝑍[7])
𝑇

𝑓𝑑(𝐻7, 𝑑) = −�̃�3𝑧3
[7]

≠ 0. 

By substituting 𝑉[7] in (3) we get: 

𝐷2𝐹𝜇(𝐻7, 𝑑)(𝑉[7], 𝑉[7]) = (𝐴𝑖𝑗) 

𝐴11 = −
2𝑆1

𝐾1
+ (𝑣3

[7]
)

2

[
4𝑎1𝐶1𝐿1𝐿3

(𝑎1+𝐿1
2)

2 �̃�1 −
𝐶1

(𝑎1+𝐿1
2)

((𝑎1 + 𝐿1
2) − 2𝐿1

2 + (𝑎1 + 2𝐿1))�̃�1], 

𝐴21 = −2𝐶 + (𝑣3
[7]

)
2

[
𝐶2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)
�̃�1

2 − (
2𝑆2𝑓1

(1+𝑓1𝐿3)2 +
𝐶2𝑎2𝐿3

(1+𝑎2𝐿2)2(1+𝑎3𝐿3)2 +

𝐶2

(1+𝑎2𝐿2)(1+𝑎3𝐿3)2) �̃�1 +
2𝑆2𝑓1

2𝐿3

(1+𝑓1𝐿3)3], 



IHJPAS. 2025, 38 (1) 
 

440 
 

𝐴31 = (𝑣3
[7]

)
2

[
2𝑔1𝐿1𝐿3

𝑎1 + 𝐿1
2 �̃�1

2 + (
𝑔1(1 − 2𝐿1

2)

𝑎1 + 𝐿1
2 + 𝑔1 − 2𝜇) �̃�1 −

2𝑎2𝑔2𝐿3

(1 + 𝑎2𝐿2)3(1 + 𝑎3𝐿3)
�̃�2

2

−
𝑔2

(1 + 𝑎2𝐿2)2(1 + 𝑎3𝐿3)2
(𝑎2𝐿2 − 1)�̃�2 −

2𝑔2𝑎3𝐿3

(1 + 𝑎2𝐿2)(1 + 𝑎3𝐿3)3
] 

Hence, it was obtained by conditions 

(𝑍[7])
𝑇

[𝐷2𝐹𝜇(𝐻7, 𝑑)(𝑉[7], 𝑉[7])] = (𝑣3
[7]

)
2

𝑧3
[7]

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

This implies that system (1)  has no (PT) bifurcation at 𝐻7  where �̃� = 𝑑 ,and has a (SN) bifurcation 

at 𝐻7 with a parameter �̃� = 𝑑. 

 

4. Numerical simulation: 

       In the following section the system's (1) dynamic behaviour has been investigated . To validate 

the findings of the study and the influence of the parameters on the dynamic model, calculations can 

be performed for a number of sets of parameters with various initial points. Figure(1)(a-b) It appears 

that the system (1) has a positive global equilibrium point at the fictitious set of parameters (19). 

 𝑆1 = 1, 𝐾1 = 0 ∙ 5   , 𝐶1 = 0 ∙ 0001, 𝑎1 = 1.5, 𝑆2 = 0.5, 𝑓1 = 0 ∙ 001, 𝐶 = 0.5 

, 𝐶2 = 0 ∙ 0001, 𝑎2 = 1.5, 𝑎3 = 1.5, 𝑔1 = 0 ∙ 00001, 𝑔2 = 0 ∙ 00001, 

𝜇 = 0 ∙ 002, 𝑑 = 0 ∙ 001                                                                                                                 (19) 
Figure 1(a − b).Time series of the system's solution (1) start with initial point (1 , 2 , 0.9 ). (a) time series of the solution 

approaches to H7 =(0.499,0.999,0.737), (b) Solution of system (1). 

In order to investigate the effects of parameters, on the dynamical behavior of the system (1), the 

system (1) has been numerically solved using the data provided in (19) Results can be achieved by 

changing one parameter at a time. 

 summarizes how additional parameters affect dynamics.  

Table 2. the system's (1) point of bifurcation . 

Range of parameter The stable point The bifurcation point 

0.1 ≤ 𝑆1 < 1 𝐻7  

0.3 < 𝐾1 < 1.99 𝐻7  

0.0001 ≤ 𝐶1 < 1.79 

1.79 ≤ 𝐶1 ≤ 2 

𝐻7 

𝐻6 

𝐶1 = 1.79 

1.5 ≤ 𝑎𝑖 < 3 i=1,2,3 

 

 𝐻7  

0.4
0.5

0.6
0.7

0.8
0.9

1

0.8
1

1.2
1.4

1.6
1.8

2
0.7

0.75

0.8

0.85

0.9

0.95

First prey

(b)

Seconed prey

P
r
e
d
a
t
o
r

Stable point

( 0.499, 0.999,0.737

Initial point

(1,2,0.9)

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.5

1

1.5

2

2.5

Time

p
o
p
u
la

ti
o
n

(a)

 

 

L1

L2

L3



IHJPAS. 2025, 38 (1) 
 

441 
 

0.5 ≤ 𝑆2 < 2 

 

 𝐻7  

0.001 ≤ 𝑑 < 0.095 

0.095 ≤ 𝑑 ≤ 0.9 

 𝐻7 

 𝐻3 

d=0.095 

0.002 ≤ 𝜇 < 0.179 

0.0179 ≤ 𝜇 < 1 

 𝐻7 

 𝐻3 

𝜇 = 0.077 

0.0001 ≤ 𝐶2 < 0.34 

0.34 ≤ 𝐶2 < 1 

 

 𝐻7 

 𝐻4 

 

𝐶2 < 0.34 

 

0.001 ≤ 𝑓1 < 0.99 

 

 𝐻7  

0.5 ≤ 𝐶 < 1 𝐻7  

0.0001 ≤ 𝑔1 < 0.00001 𝐻7  

The effect of changing the parameter  𝐶1  in the vicinity 0.0001 ≤ 𝐶1 < 1.79 the solution gets closer 

to  H7, showing in Figure 2a , increasing the range further  1.79 ≤ 𝐶1 < 2  the solution gets closer 

to   H6, showing in Figure 2b . 

 Figure 2 (a-b). (a) Time series of the system's solution (1) with  C1 = 0.0001  ,which approaches to H7 =

(0.499,0.999,0.737), and (b) time series of the solution of system (1) with  C1 = 1.79  ,which approaches to H6 =

(0, 0.999,0.812). 

 

For the parameter 𝑑 in the vicinity 0.001 ≤ 𝑑 < 0.095 the solution gets closer to H7, showing in 

Figure3a, increasing the range further 0.095 ≤ 𝑑 < 0.9  the solution gets closer to H3, showing in 

Figure 3b.  

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time

P
o
p
u
la

ti
o
n

(b)

 

 

L1

L2

L3

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.5

1

1.5

2

2.5

Time

p
o
p
u
la

ti
o
n

(a)

 

 

L1

L2

L3

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time

p
o
p
u
la

ti
o
n

(b)

 

 

L1

L2

L3

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.5

1

1.5

2

2.5

Time

p
o
p
u
la

ti
o
n

(a)

 

 

L1

L2

L3



IHJPAS. 2025, 38 (1) 
 

442 
 

Figure3. (a-b). (a) Time series of the system's solution (1.1) with d = 0.001  ,which approaches to H7 =
(0.499,0.999,0.737), and (b) time series of the solution of system (1.1) with d = 0.095  ,which approaches to H3 =
(0.5,1,0). 

For the parameter 𝜇 in the vicinity 0.002 ≤ 𝜇 < 0.179 the solution gets closer to   H7, showing in 

Figure 4a, increasing the range further  0.179 ≤ 𝜇 < 1  the solution gets closer to H3, showing in 

Figure 4b. 

Figure 4. (a-b). (a) Time series of the system's solution (1.1) with μ = 0.002  ,which approaches to H7 =
(0.499,0.999,0.737), and (b) time series of the solution of system (1.1) with μ = 0.179  ,which approaches to H3 =
(0.5,1,0). 

For the parameter  𝐶2 in the vicinity 0.0001 ≤  𝐶2 < 0.34 the solution gets closer to   H7, showing in 

Figure 5a, increasing the range further   0.34 ≤  𝐶2 < 1  the solution gets closer to H4, showing in 

Figurer 5b. 

 

 

 

 

 

 

 

 

 

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time

P
o
p
u
la

ti
o
n

(b)

 

 

L1

L2

L3

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.5

1

1.5

2

2.5

Time

p
o
p
u
la

ti
o
n

(a)

 

 

L1

L2

L3

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time

P
o
p
u
la

ti
o
n

(b)

 

 

L1

L2

L3



IHJPAS. 2025, 38 (1) 
 

443 
 

  

 

 

 

 

 

 

 

 

Figure 5. (a-b). (a) Time series of the system's solution (1.1) with  C2 = 0.0001  ,which approaches to H7 =
(0.499,0.999,0.737), and (b) time series of the solution of system (1.1) with  C2 = 0.34  ,which approaches to H4 =
(0.499,0,0.736). 

5. Conclusion 

    In this article, we present the mathematical model of three differential equations for organisms that 

characterize the effect of anti-predation behavior in the mathematical model containing fear  have 

been analyzed, which are proposed .(LB) have been examined by changing a model's parameter in 

order to examine bifurcation curves' ability to predict dynamic behavior as well as its occurrence 

states of (SN) bifurcation occurring at point 𝐻7 and (TC) and (PF) bifurcation occurring at points 

𝐻3, 𝐻4, 𝐻6The bifurcation that occurs at point is determined  . Numerical simulation is used. With 

data given in equation (19). Which are summarized as follow: 

1) A periodic dynamics for the system (1) does not exist . 

2) The parameters  𝐶1, 𝑑, 𝐶2, 𝑎𝑛𝑑 𝜇 , have an important effect in the dynamics of the system (1). 

While the others parameters not effect on the bifurcation of system 1. 

 

Acknowledgments  

We would like to express our deep appreciation to the reviewers of this work and the publishers 

of the "Ibn Al-Haitham for pure and Science Journal" for this splendid opportunity. 

 

Conflict of Interest 

     The authors declare that there are no competing interests regarding the publication of this paper. 

 

Funding 

     This work is not supported by any the Foundation. 

 

Ethical Clearance 

     Ethics of scientific research were carried out in accordance with international conditions. 

 

References 

1. Majeed AA, Ali NH. The bifurcation of stage-structured prey-predator food chain model with refuge. Int J 

Mod Res Eng Technol. 2016. 

2. Madhi ZS, Hussain MAA. Bifurcation diagram of W(U_j; τ)-function with (p, q)-parameters. Iraqi J Sci. 

2022;667-674. https://doi.org/10.24996/ijs.2022.63.2.23 

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0

0.5

1

1.5

2

2.5

Time

p
o
p
u
la

ti
o
n

(a)

 

 

L1

L2

L3

https://doi.org/10.24996/ijs.2022.63.2.23


IHJPAS. 2025, 38 (1) 
 

444 
 

3. Huang M, Ji H, Sun J, Wei L, Zha X. Bifurcation-based stability analysis of photovoltaic-battery hybrid 

power system. IEEE J Emerg Sel Top Power Electron. 2017;5:1055-1067. 

https://doi.org/10.1109/JESTPE.2017.2681125 

4. 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. 

https://doi.org/10.1109/TSG.2019.2907600 

5. 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 

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

Baghdad Sci J. 2022;19:542. http://dx.doi.org/10.21123/bsj.2022.19.3.0542 

7. Xu C, Liao M. Bifurcation behaviours in a delayed three-species food-chain model with Holling type-II 

functional response. Appl Anal. 2013;92:2468-2486. https://doi.org/10.1080/00036811.2012.742187 

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

2017;6:1143-1146. https://doi.org/10.21275/ART2017872 

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

treatment. Ibn AL-Haitham J Pure Appl Sci. 2018;510-522. 

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

10. Xia Y, Huang X, Chen F, Chen L. Stability and bifurcation of a predator-prey system with multiple anti-

predator behaviors. J Biol Syst. 2024;32(2):889-919. https://doi.org/10.1142/S021833902450030X 

11. Prasad KD, Sasmal SK. Dynamics of anti-predator behavior and effect of fear on prey-predator model. J 

Biol Syst. 2022;30(4):887-912. https://doi.org/10.1142/S0218339022500322 

12. Harianto J, Suparwati T, Dewi ALP. Local stability dynamics of equilibrium points in predator-prey models 

with anti-predator behavior. J Ilmu Dasar. 2021;22(2):153-160. https://doi.org/10.19184/jid.v22i2.23991 

13. Naji RK, Majeed SJ. The dynamical analysis of a prey-predator model with a refuge-stage structure prey 

population. Int J Differ Equ. 2016;2010464. https://doi.org/10.1155/2016/2010464 

14. Naik PA, Eskandari Z, Shahkari HE, Owolabi KM. Bifurcation analysis of a discrete-time prey-predator 

model. Bull Biomathematics. 2023;1(2):111-123. https://doi.org/10.59292/bulletinbiomath.2023006 

15. 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. 2023;152:111345. 

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

16. Mezouaghi A, Djilali S, Bentout S, Biroud K. Bifurcation analysis of a diffusive predator-prey model with 

prey social behavior and predator harvesting. Math Methods Appl Sci. 2022;45(2):718-731. 

https://doi.org/10.1002/mma.7807 

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

2019;25(3):612-626. https://doi.org/10.1177/1077546318790871 

18. Tang B, Xiao Y. Bifurcation analysis of a predator-prey model with anti-predator behaviour. Chaos Solitons 

Fractals. 2015;70:58-68. https://doi.org/10.1016/j.chaos.2014.11.008 

19. Orrell D, Smith LA. Visualizing bifurcations in high dimensional systems: The spectral bifurcation 

diagram. Int J Bifurc Chaos. 2003;13:3015-3027. 

20. Zhao Y, Zhao L, Huang C. Stability and bifurcation analysis of a fractional predator-prey model involving 

two nonidentical delays. Math Comput Simul. 2021;181:562-580. 

https://doi.org/10.1016/j.matcom.2020.10.013 

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

structured with harvesting and toxicity. J Phys Conf Ser. 2019;1362:012155. https://doi.org/10.1088/1742-

6596/1362/1/012155 

22. Hadi HR, Majeed AA. The fear and anti-predator behavior effect on the dynamics of an ecological model 

with Holling-type IV and Crowley-Martin-type of functional responses. Am Inst Phys Conf Ser. 

2024;3097(1):080032. https://doi.org/10.1063/5.0209445 

https://doi.org/10.1109/JESTPE.2017.2681125
https://doi.org/10.1109/TSG.2019.2907600
https://doi.org/10.24996/ijs.2022.63.12.35
http://dx.doi.org/10.21123/bsj.2022.19.3.0542
https://doi.org/10.1080/00036811.2012.742187
https://doi.org/10.21275/ART2017872
https://doi.org/10.30526/2017.IHSCICONF.1810
https://doi.org/10.1142/S021833902450030X
https://doi.org/10.1142/S0218339022500322
https://doi.org/10.19184/jid.v22i2.23991
https://doi.org/10.1155/2016/2010464
https://doi.org/10.59292/bulletinbiomath.2023006
https://doi.org/10.1016/j.chaos.2021.111345
https://doi.org/10.1002/mma.7807
https://doi.org/10.1177/1077546318790871
https://doi.org/10.1016/j.chaos.2014.11.008
https://doi.org/10.1016/j.matcom.2020.10.013
https://doi.org/10.1088/1742-6596/1362/1/012155
https://doi.org/10.1088/1742-6596/1362/1/012155
https://doi.org/10.1063/5.0209445