367 Β© 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 local Bifurcation of the Dynamic Behavior of Predator-Prey System with Refuge for both Species Intsar M. Kafi1* and Saad Naji2 1Department of Mathematics, College of Science, University of Baghdad, Iraq. 2Department of Mathematics, College of Science for Women, University of Baghdad, Iraq. *Corresponding Author. Received: 5 May 2023 Accepted: 4 September 2023 Published: 20 January 2025 doi.org/10.30526/38.1.3461 Abstract The main purpose of this paper is to study a predator – prey dynamical system consisting of three species prey, specialized predator and generalist predator namely H (t), I (t) and J (t) respectively, w food web and refuge for the prey and specialized predator population. The consider system has five equilibrium points 𝐴0 = (0, 0, 0), 𝐴1 = (1, 0, 0), 𝐴2 = (β„Ž βˆ’ , 𝑖 βˆ’ , 0), 𝐴3 = (β„ŽΜΏ, 0, 𝑗 ΜΏ), and the positive equilibrium point 𝐴4 = (β„ŽΜƒ, 𝑖̃, 𝑗̃ ).The stability and bifurcation of the equilibrium points was studied and the main influence was the qualitative behavior of the solution. It was found that 𝐴0 was unstable while the other equilibrium points are stable under condition so we study their bifurcation and show that 𝐴1, 𝐴2 and 𝐴3 are transcritical while 𝐴4 is saddle node bifurcation. Numerical simulations were used to illustrate the occurrence of local bifurcation of this model. Keywords: Local bifurcation, predator–prey, stability analysis, Lyapunov’s function, ecological, Refuge. 1. Introduction The mathematical study of changes in a dynamical system's qualitative asymptotic structure is known as bifurcation theory (1, 2). As well as its attempts to explain various phenomena that have been described in the natural sciences over the centuries. Where performing bifurcation analysis is often a powerful way to analyze the properties of such systems. The prey and predator model is an important topic at present as it is used to solve many problems in the ecological nature and other sciences. The prey system includes several interactions, such as competition co- existence and stage-structured (3). The system is also impacted by a number of other factors, such as shelter, sickness, and others. A bifurcation, which is the primary qualitative shift in the behavior of a dynamic system as a result of changing one of its coefficients, can occasionally https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ https://doi.org/10.30526/38.1.3462 https://orcid.org/orcid-search/search?searchQuery=0009-0002-4001-8873 mailto:intsarmkafi@gmail.com https://orcid.org/%200002-0009-6865-9798 mailto:saadnaji58@gmail.com IHJPAS. 2025, 38 (1) 368 emerge from variations in any parameter in the system, leading to complex behavior that leads to system instability. Local and global bifurcations were the two main classes that made up the bifurcation. Changes in the local stability parameters of equilibria or periodic orbits can be used to evaluate local bifurcation. Global bifurcation, on the other hand, happens when periodic orbits run into equilibrium. This leads to changes in the topology of the trajectories in phase space that cannot be contained within a limited region as is the case with local bifurcation. These bifurcations occur when a single parameter is changed (4-9). Perko (10), on the other hand, identified the prerequisites for local bifurcation, including saddle-node, transcritical, pitchfork, and period-doubling. Finally, Sotomayor's theorem(11) for local bifurcation was applied in this work to examine the occurrence of local bifurcation at equilibrium sites local bifurcation methods close to all equilibrium points, and a number of fundamental results (12-22). 2. Model formulation An ecological model was suggested for investigation in this section. A prey was included in the model, and its overall population density at time t is represented by the symbol H(t), engaging with a specialized predator I(t) is the population density at time t and generalist s assumed that the prey wa. It J(t)by 4is denoted tat time 4population density whosepredator with refuge and specialist predator is with the prey refuge. The following presumptions are now used to create the fundamental ecological model shown in Table 1: 𝑑𝐻 𝑑𝑇 = 𝜌H (1 - 𝐻 𝐾 ) – Ξ± (1-π‘š1) H I –β (1-π‘š2) H J 𝑑𝐼 𝑑𝑇 = 𝛼1 (1 - π‘š1) H I – 𝛾 (1-π‘š3) IJ –𝑑1 I (1) 𝑑𝐽 𝑑𝑇 = Ξ²1 (1 -π‘š2) H J –𝛾1 (1-π‘š3) IJ –𝑑2 J Table 1. The parameters of model (1) parameters Biological meaning 𝜌 > 0 intrinsic growth π‘˜ > 0 carrying capacity (in logistic growth) Ξ± > 0 maximum attack rate by specialist predator Ξ² > 0 maximum attack rate by generalist predator 𝑑1π‘Žπ‘›π‘‘ 𝑑2 natural death rate of specialist and generalist predator 𝛼1 > 0 maximum predation rate of the specialist predator over the prey Ξ²1 > 0 maximum predation rate of the generalist predator over the prey 𝛾 > 0 maximum attack rate by generalist predator on specialist predator 𝛾1 > 0 maximum predation rate of the generalist predator over the specialist predator 0< π‘š1 < 1, The refuge rates constants of the prey from the specialist predator 0< π‘š2 < 1 The refuge rates constants of the prey from the generalist predator 0< π‘š3 < 1 The refuge rates constants of the specialist predator from the generalist predator IHJPAS. 2025, 38 (1) 369 The above model has 13 parameters so, it is difficult to study all of them, there for reduce them a dimensionless variables and parameters are defined: π‘‘βˆ— = 𝜌 t , h = H k , i = I k , j = J k , r1 = Ξ±k(1βˆ’ π‘š1) 𝜌 , r2 = Ξ²k(1βˆ’ π‘š2) 𝜌 , r4 = Ξ³k(1βˆ’ π‘š3) 𝜌 , , π‘Ÿ5 = d1 𝜌 , r6 = Ξ²1k(1βˆ’ π‘š2) 𝜌 , π‘Ÿ7 = 𝛾1π‘˜(1βˆ’π‘š1) 𝜌 , r3 = Ξ±1k(1βˆ’ π‘š1) 𝜌 π‘Ÿ8 = d2 𝜌 . Now, for simplicity rename π‘‘βˆ— = 𝑑. So, the dimensional system (1) can be formulated as: π‘‘β„Ž 𝑑𝑑 = β„Ž[1 βˆ’ β„Ž βˆ’ r1𝑖 βˆ’ r2j] = 𝑓1( β„Ž, 𝑖 , 𝑗 ) 𝑑𝑖 𝑑𝑑 = 𝑖[ r3β„Ž βˆ’ r4𝑗 βˆ’ π‘Ÿ5 ] = 𝑓2( β„Ž, 𝑖 , 𝑗 ) 𝑑𝑗 𝑑𝑑 = 𝑗[ r6β„Ž + r7𝑖 βˆ’ r8] = 𝑓3( β„Ž, 𝑖 , 𝑗 ) } (2) With β„Ž( 0 ) β‰₯ 0 , 𝑖( 0 ) β‰₯ 0 π‘Žπ‘›π‘‘ 𝑗( 0 ) β‰₯ 0. In system (2) there are 8 parameters. All of the functions on system (2) right side are 𝐢2 (ℝ3, ℝ+). 𝑅+ 3 = {(β„Ž, 𝑖 , 𝑗 ) ∈ 𝑅3 ∢ β„Ž( 0 ) β‰₯ 0, 𝑖( 0 ) β‰₯ 0 , 𝑗( 0 ) β‰₯ 0 }. Therefore, these functions are Lipschitzian on 𝑅+ 3 , and as a result, system (2) has a unique and existing solution. Theorem 1 [Uniformly Boundedness]: 5 All the solutions of system ( 2 ) with nonnegative initial conditions are uniformly bounded. Proof: Let the solution of (2) be [β„Ž(𝑑) , 𝑖(𝑑) , 𝑗(𝑑)] the initial condition [β„Ž(0) , 𝑖(0) , 𝑗(0)] ∈ 𝑅+ 3 are nonnegative. Now, let HΜ‡(𝑑) = β„Ž(𝑑) + 𝑖(𝑑) + 𝑗(𝑑), 𝑑HΜ‡ 𝑑𝑑 < 2β„Ž βˆ’ (π‘Ÿ1 βˆ’ π‘Ÿ3)β„Žπ‘– βˆ’ (π‘Ÿ2 βˆ’ π‘Ÿ6) β„Žπ‘— βˆ’ (π‘Ÿ4 βˆ’ π‘Ÿ7) 𝑖𝑗 βˆ’ β„Ž βˆ’ π‘Ÿ5𝑖 βˆ’ π‘Ÿ8 𝑗. Now, ecologically r3 < r1, r6 < r2 and r7 < r4. 𝑑HΜ‡ 𝑑𝑑 < 2 βˆ’ Ξ΄ HΜ‡ , π‘€β„Žπ‘’π‘Ÿπ‘’ 𝛿 = π‘šπ‘–π‘› {1 , r5, r8}. Now, by solving this differential inequality for the initial value H (0) = 𝐻0 , we get that: HΜ‡(t) ≀ 2 Ξ΄ + (HΜ‡(0) βˆ’ 2 Ξ΄ ) eβˆ’Ξ΄ t Thus 0 ≀ HΜ‡(𝑑) ≀ 2 Ξ΄ as 𝑑 β†’ ∞. Hence system (2) has uniformly bounded. 3. Equilibrium points are existence and stable: There are maximum of five equilibrium points in System (2), which are listed below: β¦Ώ The point of equilibrium 𝐴0 = ( 0 ,0 ,0 ), it is always present it is referred to as the vanishing point, which is unstable and always existing. β¦Ώ The axial equilibrium5 point 𝐴1 = (1 ,0 ,0 ), existence without conditions additionally. IHJPAS. 2025, 38 (1) 370 Therefore, the characteristic equation of J (𝐴1) is as follows: 𝐽1 =𝐽(𝐴1) = [ βˆ’1 βˆ’ π‘Ÿ1 βˆ’ π‘Ÿ2 0 π‘Ÿ3 βˆ’ π‘Ÿ5 0 0 0 π‘Ÿ6 βˆ’ π‘Ÿ8 ]. (1. π‘Ž) (βˆ’1 βˆ’ πœ†) [(r3 βˆ’ π‘Ÿ5) βˆ’ πœ†] [(r6 βˆ’ π‘Ÿ8) βˆ’ πœ†] = 0. Which gives the eigenvalues of 𝐽1 𝑏𝑦: πœ†1β„Ž = βˆ’1 < 0, πœ†1𝑖 = r3 βˆ’ π‘Ÿ5 < 0 and πœ†1𝑗 = (r6 βˆ’ π‘Ÿ8) < 0 The equilibrium point 𝐴1 then becomes asymptotically stable under the following conditions: π‘Ÿ3 > π‘Ÿ5 , (3) π‘Ÿ8 > π‘Ÿ6 . (4) Otherwise, 𝐴1 is unstable. However, it is a saddle point. β¦Ώ The equilibrium7 point 𝐴2 (β„Ž βˆ’ , 𝑖 βˆ’ , 0) exists7 uniquely7 in 7𝐼𝑛𝑑. 𝑅+ 2 (Interior of 𝑅+ 2 ) of β„Ži βˆ’ plane provided that: π‘Ÿ3 > π‘Ÿ5. (5) Where: β„Ž βˆ’ = π‘Ÿ5 π‘Ÿ3 > 0. (6) 𝑖 βˆ’ = π‘Ÿ3 βˆ’ π‘Ÿ5 π‘Ÿ1 π‘Ÿ3 . (7) And the Jacobian matrix of system ( 2 ) at 𝐴2 can be written as 𝐽2 = 𝐽(𝐴2) = [ πœ‡π‘–π‘— ]3Γ—3 , ( 2. π‘Ž ) where: πœ‡11 = 1 βˆ’ 2β„Ž βˆ’ βˆ’ π‘Ÿ1 𝑖 βˆ’ , πœ‡12 = βˆ’π‘Ÿ1β„Ž βˆ’ < 0, πœ‡13 = βˆ’π‘Ÿ2β„Ž βˆ’ , πœ‡21 = π‘Ÿ3 𝑖 βˆ’ > 0, πœ‡22 = π‘Ÿ3β„Ž βˆ’ βˆ’π‘Ÿ5 , πœ‡23 = βˆ’π‘Ÿ4 𝑖 βˆ’ , , πœ‡31 = 0 , πœ‡32 = 0, πœ‡33 = π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ βˆ’ π‘Ÿ8. Consequently, the characteristic equation of J (𝐴2) is as follows: ( πœ‡33 -Ξ»)[ πœ†2 βˆ’ Θ‚)Ξ» + det (Θ‚)] = 0 where ∢ Θ‚ = [ 1 βˆ’ 2β„Ž βˆ’ βˆ’ π‘Ÿ1 𝑖 βˆ’ βˆ’π‘Ÿ1β„Ž βˆ’ π‘Ÿ3 𝑖 βˆ’ π‘Ÿ3β„Ž βˆ’ βˆ’π‘Ÿ5 ] , Then [(π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ βˆ’ π‘Ÿ8) βˆ’ Ξ»] [(1 βˆ’ 2β„Ž βˆ’ βˆ’ π‘Ÿ1 𝑖 βˆ’ βˆ’ Ξ»)(π‘Ÿ3β„Ž βˆ’ βˆ’π‘Ÿ5 βˆ’ Ξ») + π‘Ÿ1π‘Ÿ3 β„Ž βˆ’ 𝑖 βˆ’ ] = 0 Either, πœ†2β„Ž = π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ βˆ’ π‘Ÿ8 Which, as a result of the following condition, produces the first eigenvalues 𝐽2with negative real parts: π‘Ÿ8 > π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 (8) or, πœ†2 βˆ’ tr(Δ„)Ξ» + det (Δ„) = 0 Where: tr (Θ‚) = πœ†2𝑖 + πœ†2𝑗 = 1 βˆ’ 2β„Ž βˆ’ βˆ’ π‘Ÿ1 𝑖 βˆ’ + π‘Ÿ3β„Ž βˆ’ βˆ’π‘Ÿ5 = (1 + π‘Ÿ3β„Ž βˆ’ ) βˆ’ (2β„Ž βˆ’ + π‘Ÿ1 𝑖 βˆ’ +π‘Ÿ5) > 0, IHJPAS. 2025, 38 (1) 371 det (Θ‚) = πœ†2π‘–πœ†2𝑗 = (π‘Ÿ3β„Ž βˆ’ βˆ’π‘Ÿ5) (1 βˆ’ 2β„Ž βˆ’ βˆ’ π‘Ÿ1 𝑖 βˆ’ ) + π‘Ÿ1π‘Ÿ3 2π‘Ÿ5 𝑖 βˆ’ = (2π‘Ÿ5 + π‘Ÿ3)β„Ž βˆ’ + π‘Ÿ1π‘Ÿ5 𝑖 βˆ’ βˆ’ (2π‘Ÿ3 β„Ž βˆ’ 2 + π‘Ÿ5) > 0. Which, as a result of the following criteria, produces the second two eigenvalues of J2with negative real parts: 2β„Ž βˆ’ + π‘Ÿ1 𝑖 βˆ’ +π‘Ÿ5 < 1 + π‘Ÿ3β„Ž βˆ’ , (9) 2π‘Ÿ3 β„Ž βˆ’ 2 + π‘Ÿ5 < (2π‘Ÿ5 + π‘Ÿ3)β„Ž βˆ’ + π‘Ÿ1π‘Ÿ5 𝑖 βˆ’ . (10) Therefore, 𝐴2 is stable equilibrium point if conditions (8), (9) and (10) are satisfied. However otherwise, it is unstable. β¦Ώ The specialist predator free equilibrium point 𝐴3 =(β„ŽΜΏ, 0, 𝑗 ΜΏ) exists if the solutions to the following set of equations are positive: β„ŽΜΏ = π‘Ÿ8 π‘Ÿ6 > 0 (11) 𝑗̿ = π‘Ÿ6 βˆ’ π‘Ÿ8 π‘Ÿ2 π‘Ÿ6 (12) The equation (12) is positive, provided that: π‘Ÿ6 > π‘Ÿ8 . (13) The Jacobian matrix of system ( 2 ) at 𝐴2 can be written as: 𝐽3 = 𝐽(𝐴3) = [ πœ‚π‘–π‘— ]3Γ—3 , ( 3. π‘Ž ) where: πœ‚11 = 1 βˆ’ 2β„ŽΜΏ βˆ’ π‘Ÿ2𝑗,ΜΏ πœ‚12 = βˆ’π‘Ÿ1β„ŽΜΏ < 0 , πœ‚13 = βˆ’π‘Ÿ2β„ŽΜΏ < 0 πœ‚21 = 0 , πœ‚22 = π‘Ÿ3β„ŽΜΏβˆ’π‘Ÿ4 𝑗 ΜΏ βˆ’ π‘Ÿ5 , πœ‚23 = 0, πœ‚31 = π‘Ÿ6 𝑗,ΜΏ πœ‚32 = π‘Ÿ7 𝑗 ΜΏ πœ‚33 = π‘Ÿ6β„ŽΜΏ βˆ’ π‘Ÿ8. Characteristic equation for J (𝐴3) is then provided by: ( πœ‚22 - Ξ»)[ πœ†2 βˆ’ tr(Γ…)Ξ» + det (Γ…)] = 0 where ∢ Γ… = [ 1 βˆ’ 2β„ŽΜΏ βˆ’ π‘Ÿ2𝑗 ΜΏ βˆ’π‘Ÿ2β„ŽΜΏ π‘Ÿ6 𝑗 ΜΏ π‘Ÿ6β„ŽΜΏ βˆ’ π‘Ÿ8 ]. Then [(π‘Ÿ3β„ŽΜΏβˆ’π‘Ÿ4 𝑗̿ βˆ’ π‘Ÿ5) βˆ’ Ξ»] [(1 βˆ’ 2β„ŽΜΏ βˆ’ π‘Ÿ2𝑗̿ βˆ’ Ξ»)] ( π‘Ÿ6β„ŽΜΏ βˆ’ π‘Ÿ8 βˆ’ Ξ») + π‘Ÿ2π‘Ÿ6 β„ŽΜΏ 𝑗 ΜΏ = 0 Either, πœ†3β„Ž = π‘Ÿ3β„ŽΜΏβˆ’π‘Ÿ4 𝑗 ΜΏ βˆ’ π‘Ÿ5 Because of the following circumstance the first eigenvalues of J3 have negative real portions: π‘Ÿ3β„ŽΜΏ > π‘Ÿ4 𝑗 ΜΏ + π‘Ÿ5 (14) or, πœ†2 βˆ’ tr(Γ… )Ξ» + det ( Γ… ) = 0 Where, tr (Γ…) = πœ†3𝑖 + πœ†3𝑗 = 1 βˆ’ 2β„ŽΜΏ βˆ’ π‘Ÿ2𝑗̿ + π‘Ÿ6β„ŽΜΏ βˆ’ π‘Ÿ8 > 0 det (Γ…) = πœ†3π‘–πœ†3𝑗 = (1 βˆ’ 2β„ŽΜΏ βˆ’ π‘Ÿ2𝑗)ΜΏ( π‘Ÿ6β„ŽΜΏ βˆ’ π‘Ÿ8) + π‘Ÿ2π‘Ÿ6 β„ŽΜΏ 𝑗 ΜΏ IHJPAS. 2025, 38 (1) 372 = ( π‘Ÿ6 + 2 π‘Ÿ8)β„ŽΜΏ + π‘Ÿ2π‘Ÿ8 𝑗 ΜΏ βˆ’ (2π‘Ÿ6 β„ŽΜΏ 2 + π‘Ÿ8) > 0. Therefore, as a result of the following requirements, results in the second two eigenvalues of J3 having negative real portions: 2β„ŽΜΏ + π‘Ÿ2𝑗̿ βˆ’ π‘Ÿ6β„ŽΜΏ < 1, (15) 2π‘Ÿ3 β„Ž βˆ’ 2 + π‘Ÿ5 < (2π‘Ÿ5 + π‘Ÿ3)β„Ž βˆ’ + π‘Ÿ1π‘Ÿ5 𝑖 βˆ’ . (16) Therefore, 𝐴3 is stable equilibrium point if conditions (14), (15) and (16) are satisfied. On the other hand, it is unstable. β¦Ώ Finally, the positive (coexistence) equilibrium point 𝐴4 = (β„ŽΜƒ, 𝑖̃, 𝑗̃ )exists if the following system of equations has a positive solution: β„ŽΜƒ = ( π‘Ÿ4 + π‘Ÿ1 π‘Ÿ5) βˆ’ π‘Ÿ1 π‘Ÿ4 π‘Ÿ8 π‘Ÿ7⁄ ( π‘Ÿ4 + π‘Ÿ1 π‘Ÿ3) βˆ’ π‘Ÿ1 π‘Ÿ4 π‘Ÿ6 π‘Ÿ7⁄ 𝑖̃ = π‘Ÿ8 βˆ’ π‘Ÿ6β„ŽΜƒ π‘Ÿ7 , 𝑗̃ = π‘Ÿ3 β„ŽΜƒβˆ’ π‘Ÿ5 π‘Ÿ4 . Note that β„Ž = is positive, provided that: ( π‘Ÿ4 + π‘Ÿ1 π‘Ÿ5) < π‘Ÿ1 π‘Ÿ4 π‘Ÿ8 π‘Ÿ7⁄ and ( π‘Ÿ4 + π‘Ÿ1 π‘Ÿ3) < π‘Ÿ1 π‘Ÿ4 π‘Ÿ6 π‘Ÿ7⁄ . Or ( π‘Ÿ4 + π‘Ÿ1 π‘Ÿ5) > π‘Ÿ1 π‘Ÿ4 π‘Ÿ8 π‘Ÿ7⁄ π‘Žπ‘›π‘‘ ( π‘Ÿ4 + π‘Ÿ1 π‘Ÿ3) > π‘Ÿ1 π‘Ÿ4 π‘Ÿ6 π‘Ÿ7⁄ So, 𝑖̃ and 𝑗̃ are positive5, provided that: π‘Ÿ3 β„ŽΜƒ > π‘Ÿ5 and π‘Ÿ8 > π‘Ÿ6 β„ŽΜƒ respectively. For 𝐴4 = (β„ŽΜƒ, 𝑖̃, 𝑗̃), can be expressed as: 𝐽4 = 𝐽(𝐴4) = [ ñ𝑖𝑗 ]3Γ—3 , (4. π‘Ž) Where: Γ±11 = 1 βˆ’ 2β„ŽΜƒ βˆ’ π‘Ÿ1𝑖̃ βˆ’ π‘Ÿ2 𝑗̃, Γ±12 = βˆ’ π‘Ÿ1β„ŽΜƒ < 0 , Γ±13 = βˆ’ π‘Ÿ2β„ŽΜƒ < 0, Γ±21 = π‘Ÿ3𝑖̃ > 0 , 𝑛22 = π‘Ÿ3β„ŽΜƒ βˆ’ π‘Ÿ4 𝑗̃ βˆ’ π‘Ÿ5 , Γ±23 = βˆ’ π‘Ÿ4𝑖̃, Γ±31 = π‘Ÿ6 𝑗̃ > 0 , Γ±32 = π‘Ÿ7 𝑗 Μƒ> 0, Γ±33 = π‘Ÿ6β„ŽΜƒ + π‘Ÿ7 𝑖 Μƒ- π‘Ÿ8 A characteristic equation for J (𝐴4) is then provided by: πœ†3 + Ř1πœ† 2 + Ř2πœ† + Ř3 = 0, (4.b) where: Ř1 = βˆ’(Γ±11 + Γ±22 + Γ±33) , Ř2 = βˆ’[Γ±23Γ±32 βˆ’ Γ±22Γ±33 βˆ’ Γ±11(Γ±22 + Γ±33) + Γ±21Γ±12 + Γ±13Γ±31] Ř3 = βˆ’Γ±11(Γ±22Γ±33 βˆ’ Γ±23Γ±32) + Γ±12Γ±21Γ±33 βˆ’ Γ±12Γ±31Γ±23 βˆ’ Γ±13Γ±21Γ±32 βˆ’ Γ±13Γ±31Γ±22. Now, Ř1 > 0 and Ř2 > 0 provided that: 1 < 2β„ŽΜƒ + π‘Ÿ1𝑖̃ + π‘Ÿ2 𝑗,Μƒ (17) β„ŽΜƒ < π‘Ÿ4 𝑗̃ + π‘Ÿ5 , (18) π‘Ÿ6β„ŽΜƒ + π‘Ÿ7 𝑖̃ < π‘Ÿ8 (19) Also, βˆ† = Ř1 Ř2 - Ř3 > 0. By the following condition: β„ŽΜƒ < π‘Ÿ4 π‘Ÿ1π‘Ÿ3+2 π‘Ÿ4 (1 βˆ’ π‘Ÿ1 𝑖 Μƒ βˆ’ π‘Ÿ2 𝑗 Μƒ), (20) π‘Ÿ2π‘Ÿ6 β„ŽΜƒ 𝑗 Μƒ < 𝑀1 𝑀2 , (21) π‘Ÿ4π‘Ÿ6 𝑖 Μƒ 𝑗 Μƒ < π‘Ÿ1 β„ŽΜƒ(π‘Ÿ6 β„ŽΜƒ + π‘Ÿ7 𝑖 Μƒ - π‘Ÿ8) . (22) Where, 𝑀1 = 2β„ŽΜƒ + π‘Ÿ1𝑖̃ + π‘Ÿ2 𝑗̃ -1, 𝑀2 = π‘Ÿ6β„ŽΜƒ + π‘Ÿ7 οΏ½ΜƒοΏ½ - π‘Ÿ8. IHJPAS. 2025, 38 (1) 373 Using the Routh-Hurwitz criterion, however, allows each of the additional eigenvalues of eq. (4. 𝑏), have negative real parts if and only if 𝑅1 > 0, 𝑅3 > 0 and 𝑅1𝑅2 βˆ’ 𝑅3 > 0. Therefore, all of 𝐽(𝐴4) eigenvalues have a negative real portion if the additional criteria from (17) - (22) hence 𝐴4 is asymptotically stable locally. In contrast, it is unstable. 4. Local Bifurcation Analysis This section investigates the dynamical behavior of system (2) around each equilibrium point as a result of altering the parameter values. Remember that the existence of the system (2)'s non-hyperbolic equilibrium point is a required, but not sufficient, need for bifurcation. As a result, it is appropriate to apply the Sotomayor's Theorem for local bifurcation in the following theorems. Currently, in accordance with the Jacobian matrix9 of system (2). J = [ πœ•π‘“1 πœ•β„Ž πœ•π‘“1 πœ•π‘– πœ•π‘“1 πœ•π‘— πœ•π‘“2 πœ•β„Ž πœ•π‘“2 πœ•π‘– πœ•π‘“2 πœ•π‘— πœ•π‘“3 πœ•β„Ž πœ•π‘“3 πœ•π‘– πœ•π‘“3 πœ•π‘— ] (5. π‘Ž) where 𝑓𝑖 ; 𝑖 =1, 2, 3 are displayed on the system's right side (2) and πœ•π‘“1 πœ•β„Ž = 1 βˆ’ 2h βˆ’ π‘Ÿ1𝑖 βˆ’ π‘Ÿ2, πœ•π‘“1 πœ•π‘– = βˆ’ π‘Ÿ1β„Ž, πœ•π‘“1 πœ•π‘— = βˆ’ π‘Ÿ2β„Ž , πœ•π‘“2 πœ•β„Ž = π‘Ÿ3𝑖, πœ•π‘“2 πœ•π‘– = π‘Ÿ3β„Ž βˆ’ π‘Ÿ4 𝑗 - π‘Ÿ5 , πœ•π‘“2 πœ•π‘— = - π‘Ÿ4𝑖, πœ•π‘“3 πœ•β„Ž = π‘Ÿ6 𝑗, πœ•π‘“3 πœ•π‘– = π‘Ÿ7 𝑗, πœ•π‘“3 πœ•π‘— = π‘Ÿ6β„Ž + π‘Ÿ7 𝑗 - π‘Ÿ8. It is clear to 7verify that for any nonzero vector7 οΏ½Μ‡οΏ½ = (οΏ½Μ‡οΏ½1, οΏ½Μ‡οΏ½2, οΏ½Μ‡οΏ½3) 𝑇 we have: JοΏ½Μ‡οΏ½ =[ πœπ‘–π‘— ]3Γ—1 Where: 𝜏11 = (1 βˆ’ 2β„Ž βˆ’ π‘Ÿ1𝑖 βˆ’ π‘Ÿ2𝑗)οΏ½Μ‡οΏ½1 βˆ’ π‘Ÿ1β„ŽοΏ½Μ‡οΏ½2 βˆ’ π‘Ÿ2β„ŽοΏ½Μ‡οΏ½3, 𝜏21 = π‘Ÿ3𝑖�̇�1 + (π‘Ÿ3β„Ž βˆ’ π‘Ÿ4𝑗 βˆ’ π‘Ÿ5)οΏ½Μ‡οΏ½2 βˆ’ π‘Ÿ4𝑖�̇�3, 𝜏31 = π‘Ÿ6𝑗�̇�1 + π‘Ÿ7𝑗�̇�2 + (π‘Ÿ6β„Ž + π‘Ÿ7𝑖 βˆ’ π‘Ÿ8)οΏ½Μ‡οΏ½3. D2β„‰πœ‡(Ý, πœ‡)(οΏ½Μ‡οΏ½ , οΏ½Μ‡οΏ½) = [ οΏ½ΜˆοΏ½π‘–π‘— ]3Γ—1. (23) Where: �̈�11 = (βˆ’2οΏ½Μ‡οΏ½1 βˆ’ π‘Ÿ1οΏ½Μ‡οΏ½2 βˆ’ π‘Ÿ2οΏ½Μ‡οΏ½3)οΏ½Μ‡οΏ½1 βˆ’ π‘Ÿ1οΏ½Μ‡οΏ½1οΏ½Μ‡οΏ½2 βˆ’ π‘Ÿ2οΏ½Μ‡οΏ½1οΏ½Μ‡οΏ½3 , �̈�21 = π‘Ÿ3οΏ½Μ‡οΏ½1οΏ½Μ‡οΏ½2 + (π‘Ÿ3οΏ½Μ‡οΏ½1 βˆ’ π‘Ÿ4οΏ½Μ‡οΏ½3)οΏ½Μ‡οΏ½2 βˆ’ π‘Ÿ4οΏ½Μ‡οΏ½2οΏ½Μ‡οΏ½3, �̈�31 = π‘Ÿ6οΏ½Μ‡οΏ½1οΏ½Μ‡οΏ½3 + π‘Ÿ7οΏ½Μ‡οΏ½2οΏ½Μ‡οΏ½3 + (π‘Ÿ6οΏ½Μ‡οΏ½1 + π‘Ÿ7οΏ½Μ‡οΏ½2)οΏ½Μ‡οΏ½3. Where Ý = (β„Ž, 𝑖, 𝑗)𝑇 and πœ‡ is any bifurcation parameter. Theorems in the following the local bifurcation conditions near equilibrium points are established. 4.1 7Local 7bifurcation analysis7 near π‘¨πŸ: Theorem (2): If the value of the parameter π‘Ÿ3 passes through π‘Ÿ3̈ = π‘Ÿ3 then, system (2) at the axial equilibrium point 𝐴1 = (1, 0, 0) possesses: β€’ No saddle-node bifurcation9. β€’ 7Transcritical bifurcation7. Proof: According to the Jacobian matrix 𝐽(𝐴1) given by eq.(1. π‘Ž ): Zero eigenvalue exists for system (2) at equilibrium point 𝐴1= (1,0,0). (say πœ†1𝑖 = 0 ) at π‘Ÿ3 = π‘Ÿ3̈ , and the Jacobian matrix 𝐽1 with π‘Ÿ3 = π‘Ÿ3̈ becomes: IHJPAS. 2025, 38 (1) 374 𝐽1̈ = 𝐽1( π‘Ÿ3 = π‘Ÿ3̈ ) = [ βˆ’1 βˆ’ π‘Ÿ1 βˆ’ π‘Ÿ2 0 0 0 0 0 π‘Ÿ3 βˆ’ π‘Ÿ5 ], Now, let οΏ½Μ‡οΏ½[1] = (οΏ½Μ‡οΏ½1 [1], οΏ½Μ‡οΏ½2 [1] , οΏ½Μ‡οΏ½3 [1] ) 𝑇 be the eigenvector7 corresponding7 to the 7eigenvalue πœ†1𝑖 = 0. Thus ( 𝐽1̈ βˆ’ πœ†1𝑖𝐼) οΏ½Μ‡οΏ½ [1] = 0, which gives: 𝜌[1] = (βˆ’π‘Ÿ1οΏ½Μ‡οΏ½2 [1], οΏ½Μ‡οΏ½2 [1], 0) where οΏ½Μ‡οΏ½2 [1] and οΏ½Μ‡οΏ½2 [1] are any 7nonzero real number. Let Γ‡[1] = (Γ‡1 [1] , Γ‡2 [1] , Γ‡3 [1]) 𝑇 be the 7eigenvector 7associated with the 7eigenvalue πœ†1𝑖 = 0 of the 7matrix [ J̈1] 𝑇 . Then we have, ( 𝐽1̈ 𝑇 βˆ’ πœ†1𝑖𝐼) Γ‡ [1] = 0. By solving7 this 7equation Γ‡[1], We obtain, Γ‡[1] = (0, Γ‡2 [1], 0) 𝑇 where Γ‡2 [1] 7any 7nonzero real number. Now, βˆ‚π‘“ βˆ‚π‘Ÿ3 = π‘“π‘Ÿ3(𝐴1, π‘Ÿ3) = ( βˆ‚π‘“1 βˆ‚π‘Ÿ3 , βˆ‚π‘“2 βˆ‚π‘Ÿ3 , βˆ‚π‘“3 βˆ‚π‘Ÿ3 ) T = (0, h𝑖, 0)T. So, βˆ‚π‘“ βˆ‚π‘Ÿ3 (𝐴1, �̈�3) = (0, 0, 0) 𝑇 that's why (Γ‡)𝑇 βˆ‚π‘“ βˆ‚π‘Ÿ3 (𝐴1, �̈�3) = 0. Therefore, 7according to 7Sotomayor’s theorem7 the saddle7-7node 7bifurcation cannot occur. While the first 7condition of 7transcritical 7bifurcation is satisfied7. Now, since π·π‘“π‘Ÿ3( 𝐴1, �̈�3)οΏ½Μ‡οΏ½ [1] = ( 0 0 0 0 1 0 0 0 0 )( βˆ’π‘Ÿ1οΏ½Μ‡οΏ½2 [1] οΏ½Μ‡οΏ½2 [1] 0 ) = ( 0 βˆ’οΏ½Μ‡οΏ½2 [1] 0 ), (Γ‡[1]) 𝑇 [π·π‘“π‘Ÿ3( 𝐴1 , �̈�3)οΏ½Μ‡οΏ½ [1]] = (0, Γ‡2 [1], 0) 𝑇 ( 0 οΏ½Μ‡οΏ½2 [1] 0 ) = Γ‡2 [1]𝜌 Μ‡ 2 [1] β‰  0. Now, by 7substituting οΏ½Μ‡οΏ½[1] in (23), we get: 𝐷2𝑓( 𝐴1, �̈�3)(οΏ½Μ‡οΏ½ [1], οΏ½Μ‡οΏ½[1]) = ( 4π‘Ÿ1 2[οΏ½Μ‡οΏ½2 [1]]2 βˆ’ π‘Ÿ1π‘Ÿ3[οΏ½Μ‡οΏ½2 [1]]2 0 ). Hence, it is 7obtained that: (Γ‡[1]) 𝑇 𝐷2𝑓( 𝐴1 , �̈�3)(οΏ½Μ‡οΏ½ [1], οΏ½Μ‡οΏ½[1]) = (0 , Γ‡2 [1] , 0) 𝑇 ( 4π‘Ÿ1 2[οΏ½Μ‡οΏ½2 [1]]2 βˆ’ π‘Ÿ1π‘Ÿ3[οΏ½Μ‡οΏ½2 [1]]2 0 ) = βˆ’ π‘Ÿ1π‘Ÿ3[οΏ½Μ‡οΏ½2 [1]]2 Γ‡2 [1] β‰  0 Thus, 7according to 7Sotomayor’s 7theorem system ( 2 ) has 7transcrirtical 7bifurcation but not 7experience a 7saddl node 7bifurcation at 𝐴1 with the 7parameter π‘Ÿ3, where π‘Ÿ3 = π‘Ÿ3̈ . 4.2 Near by local bifurcation analysisπ‘¨πŸ ( 𝐑 βˆ’ , π’Š βˆ’ , 𝟎 ) : Theorem (3): Assume that the following conditions are met: €2 β‰  0, (24) €4 β‰  0, (25) πœ— = βˆ’2€1€3(€1 + π‘Ÿ1 + π‘Ÿ2€2) + π‘Ÿ3€1 βˆ’ π‘Ÿ4€2 + €2€4(π‘Ÿ6€1 + π‘Ÿ7) β‰  0 (26) IHJPAS. 2025, 38 (1) 375 €1 = h βˆ’ [π‘Ÿ1π‘Ÿ4 h βˆ’ 𝑖 βˆ’ βˆ’ π‘Ÿ2(π‘Ÿ3β„Ž βˆ’ βˆ’ π‘Ÿ3)] 𝑖 [π‘Ÿ4 (1 βˆ’ π‘Ÿ1 𝑖 βˆ’ ) βˆ’ (2π‘Ÿ4 + π‘Ÿ2π‘Ÿ3)h βˆ’ ] βˆ’ , €2 = ( 1 βˆ’ 2 h βˆ’ βˆ’ π‘Ÿ1 𝑖 βˆ’ )𝑖 βˆ’ βˆ’ π‘Ÿ1β„Ž βˆ’ π‘Ÿ2β„Ž βˆ’ €3 = π‘Ÿ3β„Ž βˆ’ βˆ’ π‘Ÿ5 π‘Ÿ1β„Ž βˆ’ βˆ’ , €4 = π‘Ÿ2β„Ž βˆ’ €3 + π‘Ÿ4 𝑖 βˆ’ π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ βˆ’ π‘Ÿ8 . Then system (2) at the equilibrium5 point 𝐴2 = ( h βˆ’ , 𝑖 βˆ’ ,0 ) with the parameter5 π‘Ÿ8 ̈ = π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ 5possesses: β€’ No 5saddle-5node bifurcation. β€’ Transcritical5 5bifurcation. Proof: According5 to the Jacobian5 matrix 𝐽(𝐴2) given by eq.(2. π‘Ž ) of system ( 2 ) at the 5equilibrium point 𝐴2 = ( h βˆ’ , 𝑖 βˆ’ ,0 ) has zero eigenvalue (say πœ†2𝑗 = 0 ) at π‘Ÿ8 = �̈�8 , and the 5Jacobian matrix 𝐽2 with �̈�8 = π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ becomes: 𝐽2̈ = 𝐽2(π‘Ÿ8 = �̈�8) = [πœ‡π‘–π‘—π‘–π‘—]3x3 , where, πœ‡π‘–π‘—Μ‡ ij = πœ‡π‘–π‘— for all i, j =1,2,3 except5 πœ‡π‘–π‘—Μ‡ 33 = π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ βˆ’ π‘Ÿ8=0. Let οΏ½Μ‡οΏ½[2] = (οΏ½Μ‡οΏ½1 [2] , οΏ½Μ‡οΏ½2 [2] , οΏ½Μ‡οΏ½3 [2] ) 𝑇 be the 5eigenvector corresponding5 to the eigenvalue πœ†2𝑗 = 0. Thus (𝐽2̈ βˆ’ πœ†2𝑗𝐼) οΏ½Μ‡οΏ½ [2] = 0, which gives5: οΏ½Μ‡οΏ½1 [2] = €1οΏ½Μ‡οΏ½2 [2] and οΏ½Μ‡οΏ½3 [2] = €2οΏ½Μ‡οΏ½2 [2] , where οΏ½Μ‡οΏ½2 [2] any nonzero, real number with €1, €2 and which and which are mentioned in the state of the theorem. Let Γ‡[2] = (Γ‡1 [2] , Γ‡2 [2] , Γ‡3 [2]) 𝑇 be the 4eigenvector 4associated with the 4eigenvalue πœ†2𝑀 = 0 of the 4matrix j ̈ 2 . Then we have ( j ̈ 2 βˆ’ Ξ»2wI) Γ‡ [2] = 0. By 4solving this 4equation for Γ‡[2], we 4obtain Γ‡[2] = (€3Γ‡2 [2] , Γ‡2 [2] , €4Γ‡2 [2]) 𝑇 , where Γ‡2 [2] any real numbers that are not zero, with €3, €4 𝑀hich are 54mentioned in the state of the theorem45. Now, consider: βˆ‚π‘“ βˆ‚π‘Ÿ8 (Ý, π‘Ÿ8) = π‘“π‘Ÿ8(Ý , π‘Ÿ8) = ( βˆ‚π‘“1 βˆ‚π‘Ÿ8 , βˆ‚π‘“2 βˆ‚π‘Ÿ8 , βˆ‚π‘“3 βˆ‚π‘Ÿ8 ) T = (0 , 0 , 𝑗)𝑇 So, βˆ‚π‘“ βˆ‚π‘Ÿ8 (A2 , π‘Ÿ8) = (0 , 0, 0 )𝑇 . That's why (Γ‡[2]) 𝑇 π‘“π‘Ÿ8(A2, �̈�8) = 0. IHJPAS. 2025, 38 (1) 376 Therefore9, according to Sotomayor9’s theorem there can be no saddle-node bifurcation. Although the first need for transcritical bifurcation has been satisfied. Now, since π·π‘“π‘Ÿ8( 𝐴2 , �̈�8)οΏ½Μ‡οΏ½ [2] = ( 0 0 0 0 0 0 0 0 βˆ’ 1 )( €1οΏ½Μ‡οΏ½2 [2] οΏ½Μ‡οΏ½2 [2] €2οΏ½Μ‡οΏ½2 [2] ) = ( 0 0 0 βˆ’β‚¬2οΏ½Μ‡οΏ½2 [2] ) (Γ‡[2]) 𝑇 [π·π‘“π‘Ÿ8( 𝐴2 , �̈�8)οΏ½Μ‡οΏ½ [2]] = (€3Γ‡2 [2] , Γ‡2 [2] , €4Γ‡2 [2]) ( 0 0 βˆ’β‚¬2οΏ½Μ‡οΏ½2 [2] ) = βˆ’β‚¬2€4Γ‡2 [2]οΏ½Μ‡οΏ½2 [2] β‰  0. Moreover, by 4substituting οΏ½Μ‡οΏ½[2] in (23), we get: 𝐷2𝑓(𝐴2 , �̈�8)(οΏ½Μ‡οΏ½ [2], οΏ½Μ‡οΏ½[2]) = [ βˆ’2€1(οΏ½Μ‡οΏ½2 [2])2 (€1 + π‘Ÿ1 + π‘Ÿ2€2) 2(οΏ½Μ‡οΏ½2 [2])2(π‘Ÿ3€1 βˆ’ π‘Ÿ4€2) 2€2(οΏ½Μ‡οΏ½2 [2])2(π‘Ÿ6€1 + π‘Ÿ7) ]. Hence, it is 4obtained that: (€3Γ‡2 [2] , Γ‡2 [2] , €4Γ‡2 [2]) (Γ‡[2]) 𝑇 𝐷2𝑓(𝐴2 , �̈�8)(οΏ½Μ‡οΏ½ [2] , οΏ½Μ‡οΏ½[2]) = 2πœ—(οΏ½Μ‡οΏ½2 [2])2Γ‡2 [2] Where, €1, €3 π‘Žπ‘›π‘‘ €4 are mentioned in the theorem's state. As a result of condition (26), we get that: (Γ‡[2]) 𝑇 𝐷2𝑓(𝐴2 , �̈�8)(οΏ½Μ‡οΏ½ [2] , οΏ½Μ‡οΏ½[2]) β‰  0 . Thus, 4according to 4Sotomayor’s 4theorem system (2) has a 4transcritical bifurcation4 at the 4equilibrium point 𝐴2 = (h βˆ’ , 𝑖 βˆ’ , 0 ) with the 4parameter �̈�8 = π‘Ÿ6β„Ž βˆ’ + π‘Ÿ7 𝑖 βˆ’ . 4.3 Local4 bifurcation4 analysis4 near π‘¨πŸ‘ = (οΏ½ΜΏοΏ½, 𝟎, 𝒋 ΜΏ): Theorem (4): Assume that the following criteria are fulfilled: π‘Ÿ4 > π‘Ÿ3 β„Ž ΜΏ βˆ’ π‘Ÿ5 𝑗 ΜΏ , (27) νœ€2Μ€ β‰  0, (28) ℝ β‰  0 (29) where: νœ€1Μ€ = π‘Ÿ1π‘Ÿ6β„ŽΜΏ + π‘Ÿ7𝑗(ΜΏ1 βˆ’ 2β„ŽΜΏ βˆ’ π‘Ÿ2) β„ŽΜΏ[π‘Ÿ2 + π‘Ÿ1(π‘Ÿ6β„ŽΜΏ βˆ’ π‘Ÿ8)] , νœ€2Μ€ = 1 βˆ’ π‘Ÿ2 𝑗 ΜΏ βˆ’ (2 βˆ’ π‘Ÿ2νœ€1Μ€)β„ŽΜΏ π‘Ÿ1 β„ŽΜΏ , νœ€3Μ€ = π‘Ÿ7𝑗̿ π‘Ÿ1β„ŽΜΏ , ℝ = Γ‡3 [3][€2(π‘Ÿ6€1 + π‘Ÿ7) βˆ’ €1€3(€1 + π‘Ÿ1 + π‘Ÿ2€2)] + Γ‡2 [3][π‘Ÿ3€1 βˆ’ π‘Ÿ4€2]. Then system ( 2 ) at the 4equilibrium point 𝐴3 = (β„ŽΜΏ, 0, 𝑗 ΜΏ) with the 4parameter value: �̈�4 = π‘Ÿ3β„ŽΜΏ βˆ’ π‘Ÿ4𝑗 ΜΏ has a transcritical bifurcation, but a saddle4 βˆ’ node cannot 4occur at 𝐴3. Proof: The characteristic equation represented by eq. ( 3. π‘Ž ) of 𝑠ystem (2) at the equilibrium point 𝐴3 has zeroeigenvalue (say πœ†3𝑖 = 0 ) at π‘Ÿ4 = �̈�4 and the Jacobian matrix 𝐽3 with 4parameter π‘Ÿ4 = �̈�4 becomes: IHJPAS. 2025, 38 (1) 377 𝐽3̈ = 𝐽3(π‘Ÿ4 = �̈�4) = [�̃�𝑖𝑗 ]3Γ—3 , Where: �̃�𝑖𝑗 = πœ‚π‘–π‘— for all i , j = 1,2,3 except �̃�𝑖𝑗 = π‘Ÿ3β„ŽΜΏ βˆ’ �̈�4𝑗̿ βˆ’ π‘Ÿ4 = 0. Let οΏ½Μ‡οΏ½[3] = ((οΏ½Μ‡οΏ½2 [3] , οΏ½Μ‡οΏ½2 [3] , οΏ½Μ‡οΏ½3 [3] ) 𝑇 be the eigenvector which follows the eigenvalue be the eigenvector which follows the eigenvalue. πœ†3𝑠 = 0. Thus (𝐽3̈ βˆ’ πœ†3𝑠𝐼) οΏ½Μ‡οΏ½ [3] = 0 , which gives: οΏ½Μ‡οΏ½[3] = (οΏ½Μ‡οΏ½1 [3], νœ€2Μ€ οΏ½Μ‡οΏ½1 [3] , νœ€1Μ€ οΏ½Μ‡οΏ½1 [3]) 𝑇 , where οΏ½Μ‡οΏ½1 [3] any nonzero4 number with, νœ€1Μ€ which are mentioned4 in the state4 of the theorem. Let Γ‡[3] = (Γ‡1 [3] , Γ‡2 [3] , Γ‡3 [3] ) 𝑇 become the eigenvector linked to the eigenvalue πœ†3𝑠 = 0 of the matrix 𝐽3̈. Then we have (𝐽3̈ 𝑇 βˆ’ πœ†3𝑖𝐼) Γ‡ [3] = 0. By solving this equation for Γ‡[3], we obtain: Γ‡[3] = (νœ€3Μ€ Γ‡3 [3], Γ‡2 [3], Γ‡3 [3] ) 𝑇 where Γ‡3 [3] any nonzero number9 with νœ€3Μ€ those are referred to in the theorem's state. Now, consider: βˆ‚π‘“ βˆ‚π‘Ÿ4 (Ý, π‘Ÿ4) = π‘“π‘Ÿ4(Ý , π‘Ÿ4) = ( βˆ‚π‘“1 βˆ‚π‘Ÿ4 , βˆ‚π‘“2 βˆ‚π‘Ÿ4 , βˆ‚π‘“3 βˆ‚π‘Ÿ4 ) T = (0,βˆ’π‘–π‘—, 0)𝑇 So, βˆ‚π‘“ βˆ‚π‘Ÿ4 (A3 , π‘Ÿ4) = (0 , 0, 0 ) 𝑇 . And hence (Γ‡[3]) 𝑇 π‘“π‘Ÿ4(A3, �̈�4) = 0. Therefore, the saddle-node bifurcation is ruled out by Sotomayor's theorem. While the first transcritical bifurcation condition has been satisfied. Now, sinceοΏ½Μ‡οΏ½[3] = [οΏ½Μ‡οΏ½1 [3], νœ€2Μ€ οΏ½Μ‡οΏ½1 [3] , νœ€1Μ€ οΏ½Μ‡οΏ½1 [3]] π·π‘“π‘Ÿ4( 𝐴3 , �̈�4)οΏ½Μ‡οΏ½ [3] = ( 0 0 0 0 βˆ’ 𝑗 0 0 0 0 )( οΏ½Μ‡οΏ½1 [3] νœ€2Μ€ οΏ½Μ‡οΏ½1 [3] νœ€1Μ€ οΏ½Μ‡οΏ½1 [3] ) = ( 0 βˆ’π‘— νœ€2Μ€ οΏ½Μ‡οΏ½1 [3] 0 ) (Γ‡[2]) 𝑇 [π·π‘“π‘Ÿ4( 𝐴2 , �̈�4)οΏ½Μ‡οΏ½ [2]] = (€3Γ‡3 [3] , Γ‡2 [3], Γ‡3 [3]) ( 0 βˆ’π‘— νœ€2Μ€ οΏ½Μ‡οΏ½1 [3] 0 ) = βˆ’νœ€2Μ€ Γ‡2 [3]𝑗 ΜΏ οΏ½Μ‡οΏ½1 [3]. So, by condition (28), we obtain that: (Γ‡[2]) 𝑇 [π·π‘“π‘Ÿ4( 𝐴2 , �̈�4)οΏ½Μ‡οΏ½ [2]] β‰  0 Moreover, by 4substituting οΏ½Μ‡οΏ½[3] in (23), we get: 𝐷2𝑓(𝐴3 , �̈�4)(οΏ½Μ‡οΏ½ [3], οΏ½Μ‡οΏ½[3]) = [ βˆ’2€1(οΏ½Μ‡οΏ½2 [2])2 (€1 + π‘Ÿ1 + π‘Ÿ2€2) 2(οΏ½Μ‡οΏ½2 [2])2(π‘Ÿ3€1 βˆ’ π‘Ÿ4€2) 2€2(οΏ½Μ‡οΏ½2 [2])2(π‘Ÿ6€1 + π‘Ÿ7) ]. Hence, it is 4obtained that: (€3Γ‡3 [3] , Γ‡2 [3] , Γ‡3 [3]) IHJPAS. 2025, 38 (1) 378 (Γ‡[2]) 𝑇 𝐷2𝑓(𝐴3 , �̈�4)(οΏ½Μ‡οΏ½ [3] , οΏ½Μ‡οΏ½[3])=2(οΏ½Μ‡οΏ½2 [2])2 ℝ Where, €1, €2, €3 and ℝ are mentioned in the state of the theorem. As a result of condition (29), we get that: (Γ‡[2]) 𝑇 𝐷2𝑓(𝐴2 , �̈�8)(οΏ½Μ‡οΏ½ [2] , οΏ½Μ‡οΏ½[2]) β‰  0 . Thus, 4according to 4Sotomayor’s 4theorem system (2) has a 4transcritical bifurcation4 at the 4equilibrium point 𝐴3 = (β„ŽΜΏ, 0, 𝑗 ΜΏ) with the 4parameter �̈�4 = π‘Ÿ3β„ŽΜΏ βˆ’ π‘Ÿ4𝑗.ΜΏ 4.4 Local bifurcation4analysis near π‘¨πŸ’ = (οΏ½ΜƒοΏ½, οΏ½ΜƒοΏ½, 𝒋̃ ): Theorem (5): Assume that the following criteria are fulfilled: €1= βˆ’π‘Ÿ3 π‘Ÿ6 , €2= π‘Ÿ8βˆ’π‘Ÿ6hΜƒ βˆ’(π‘Ÿ7οΏ½ΜƒοΏ½+𝒋̃€1) π‘Ÿ6 οΏ½ΜƒοΏ½ , €3 = βˆ’π‘Ÿ3 οΏ½ΜƒοΏ½ π‘Ÿ6 jΜƒ , €4= (π‘Ÿ6hΜƒ+π‘Ÿ7οΏ½ΜƒοΏ½ βˆ’π‘Ÿ8) €3βˆ’ π‘Ÿ4 𝒋̃ π‘Ÿ2 hΜƒ β‰  0, (30) 1 = 2 hΜƒ + π‘Ÿ1𝑖̃ + π‘Ÿ2 jΜƒ, (31) €2€4 [€2 + π‘Ÿ1€1 + π‘Ÿ2] β‰  €1[π‘Ÿ3€2 + π‘Ÿ4] + €4[π‘Ÿ6€2 + π‘Ÿ7€1]. (32) Then system ( 2 ) at the 4equilibrium point 𝐴4 = (β„ŽΜƒ, 𝑖̃, 𝑗̃ ) with the 4parameter value: π‘Ÿ1̈ = 1 βˆ’ 2 hΜƒ βˆ’ π‘Ÿ2 jΜƒ 𝑖̃ , has a saddle βˆ’ node bifurcation, but neither a transcritical4nor a pitchfork bifurcation4 at 𝐴4. Proof: The characteristic equation given by eq. (4.a) if 𝐻4 = 0 and 𝐴4 becomes a non- hyperbolic equilibrium point, of system (2) having zero eigenvalue (say πœ†4β„Ž = 0). The Jacobian matrix for system (2) at equilibrium point 𝐸4 with parameter( π‘Ÿ1 = π‘Ÿ1̈) clearlybecomes: 𝐽4̈ = 𝐽4( π‘Ÿ1 = π‘Ÿ1̈) = [�̿�𝑖𝑗 ]3Γ—3 , where, �̿�𝑖𝑗 = Γ±11 for all i, j =1, 2, 3 except Γ±11 which is given by: Γ±11 = 1 βˆ’ 2β„ŽΜƒ βˆ’ π‘Ÿ1̈ 𝑖̃ βˆ’ π‘Ÿ2 jΜƒ Let οΏ½Μ‡οΏ½[4] = (οΏ½Μ‡οΏ½1 [4] , οΏ½Μ‡οΏ½2 [4] , οΏ½Μ‡οΏ½3 [4] ) 𝑇 be the eigenvector9 that follows the 9eigenvalue πœ†4β„Ž = 0 . Thus (𝐽4̈ βˆ’ πœ†4β„ŽπΌ) οΏ½Μ‡οΏ½ [4] = 0, which gives: οΏ½Μ‡οΏ½[4] = (€2οΏ½Μ‡οΏ½3 [4] , €1οΏ½Μ‡οΏ½3 [4], οΏ½Μ‡οΏ½3 [4]) 𝑇 , where οΏ½Μ‡οΏ½3 [4] any non-zero value that has the €1 and €2 conditions given in the theorem. Let Γ‡[4] = (Γ‡[4] , Γ‡[4] , Γ‡[4]) 𝑇 be the eigenvector4 associated4 with the 4eigenvalue πœ†4𝑖 = 0 of the matrix 𝐽4̈. Next, we have (𝐽4̈ 𝑇 βˆ’ πœ†4𝑖𝐼) Γ‡ [4] = 0. By solving this 4equation for βˆ…[4], we obtain: Γ‡[4] = (€4Γ‡2 [4] , Γ‡2 [4], €3Γ‡2 [4] ) 𝑇 where Γ‡2 [4] any 4nonzero 4number with €3 and €4 which are 4mentioned in the state4 of the theorem. Now, IHJPAS. 2025, 38 (1) 379 πœ•π‘“ πœ•π‘Ÿ1 = π‘“π‘Ÿ1(Ý , π‘Ÿ1) = ( πœ•π‘“1 πœ•π‘Ÿ1 , πœ•π‘“2 πœ•π‘Ÿ1 , πœ•π‘“3 πœ•π‘Ÿ1 , ) 𝑇 = (βˆ’β„Žπ‘–, 0, 0)𝑇 , So, π‘“π‘Ž1(𝐴4 , π‘Ÿ1̈) = (βˆ’β„ŽΜƒ 𝑖̃ , 0, 0) 𝑇 , and hence, it is 4obtained that: (Γ‡[4]) 𝑇 π‘“π‘Ÿ1(𝐴4 , π‘Ÿ1̈) = βˆ’β„ŽΜƒ 𝑖̃ €4Γ‡2 [4] . In the context of condition (30), we thus have that: (Γ‡[4]) 𝑇 π‘“π‘Ÿ1(𝐴4 , π‘Ÿ1̈) β‰  0. 𝐷2𝑓(𝐴4 , π‘Ÿ1̈)(οΏ½Μ‡οΏ½ [4] , οΏ½Μ‡οΏ½[4]) = [ βˆ’2 €2 (οΏ½Μ‡οΏ½3 [4])2[€2 + π‘Ÿ1€1 + π‘Ÿ2] 2 €1 (οΏ½Μ‡οΏ½3 [4])2[π‘Ÿ3€2 + π‘Ÿ4] 2 (οΏ½Μ‡οΏ½3 [4])2[π‘Ÿ6€2 + π‘Ÿ7€1] ] . Hence, it is 4obtained that: (Γ‡[4]) 𝑇 𝐷2𝑓(𝐴4 , π‘Ÿ1̈)(οΏ½Μ‡οΏ½ [4] , οΏ½Μ‡οΏ½[4]) = βˆ’2 (οΏ½Μ‡οΏ½3 [4])2Γ‡2 [4] (€2€4 [€2 + π‘Ÿ1€1 + π‘Ÿ2] βˆ’ €1[π‘Ÿ3€2 + π‘Ÿ4] βˆ’ €4[π‘Ÿ6€2 + π‘Ÿ7€1]). Therefore, in accordance with condition (32) we get that: (Γ‡[4]) 𝑇 𝐷2𝑓(𝐴4 , π‘Ÿ1̈)(οΏ½Μ‡οΏ½ [4] , οΏ½Μ‡οΏ½[4]) β‰  0 Sotomayor's theorem is used to show that system (2) has a saddle-node bifurcation 𝐴4 = (β„ŽΜƒ, 𝑖̃, j Μƒ ) at π‘Ÿ1̈. 5. Numerical Simulation This section numerically analyzes System (2)'s dynamical behavior for a given set of parameters and various initial point sets. These are the study's objectives, including: 1. Analyze the impact of changing the value of each parameter on the system's dynamic behavior (2). 2. Verify the analytical results that were found. The following hypothetical collection of parameters is found to satisfy the stability requirements of the positive equilibrium point, system ( 2 ) has a globally asymptotically9 stable positive equilibrium point as shown in Figure πŸ“. 𝟏 (𝒂, 𝒃, 𝒄), (5.1) π’“πŸ = 𝟎. πŸŽπŸ“ , π’“πŸ = 𝟎. 𝟏 , π’“πŸ‘ = 𝟎. πŸŽπŸ‘, π’“πŸ’= 𝟎. 𝟏 , π’“πŸ“= 𝟎. 𝟎𝟏 , π’“πŸ” = 𝟎. πŸŽπŸ” , π’“πŸ• = 𝟎. πŸŽπŸ• , π’“πŸ– = 𝟎. 𝟏 IHJPAS. 2025, 38 (1) 380 Figure1: The time series of the solution of system (2) started from the three different initial point (0.8, 0.7, 1.64), (1.2, 0.7, 1.64) and (0.2, 0.4, 0.64) for the data given by (5.1), (a) the trajectories of h as a function of time, (b) the trajectories of i as a function of time, (c) the trajectories of j as a function of time. As the solution of system (2) approaches asymptotically to the positive equilibrium point A4= (0.4544, 0.6567, 2. 1722) beginning with three distinct starting points, Figure1. clearly demonstrates that system (2) has a globally asymptotically stable, and this is supporting our obtained analytical results. We will now talk about how system's (2) parameter settings affect the system's dynamical behavior. The system is numerically solved for the data in (5.1) by changing one parameter at a time, sometimes even two, and the results are shown below. The system will get closer to the point of positive equilibrium A4 in the interior of the positive quadrant of the hij-plane as shown in Figure2. a1 when the predation rate on a prey is varied from the specialist predator in the range 0.0001 < r1 ≀ 0.595 while maintaining other parameters as data given in (5.1), r1 has a usual value of 0.15. As shown in Figure2.a2 for average value r1=0.9, it is seen that the solution of system (2) approaches asymptotically to the equilibrium point A2 in the range 0.595 < r1 < 2 Figure2.a1: The time series of the solution of system (2) approaches asymptotically to the positive equilibrium point A4= (0.8182, 0.7275, 0.7271) in the interior of 𝑅+ 3 . For the data in (5.1) with r1 = 0.15. Figure2.a2: The time series of the solution of system (2) approaches asymptotically to the positive equilibrium point A2= (0.3333, 0.7407, 0) in the interior of 𝑅+ 3 . For the data in (5.1) with r1 = 0.9. Additionally, changing the specialist predator's mortality death rate between 0.0001 and 0.03 while maintaining the other parameters according to the data in (5.1) results in the specialized predator going extinct, and Figure3.e1 illustrates how system (2) approaches asymptotically to the positive equilibrium pointA4 with a typical value of r5=0.01, however, for a typical value of 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 Time Po pu la tio n (1.1) h i j 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 Time Po pu la tio n (1.2) h i j 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Time Po pu lat ion (1.3) h i j 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Time Po pu lat ion (1.3) h i j 0 1 2 3 4 5 6 7 8 9 10 x 10 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time Po pu la tio n (8.1) h i j IHJPAS. 2025, 38 (1) 381 r5=0.1, as shown in Figure3.e2, at 0.03< r5 <0.15 approaches asymptotically to the axial equilibrium point A1, additionally, for a typical value ofr5=0.5, Figure3.e3 shows that for 0.15≀ r5 <0.6 approaches asymptotically to the equilibrium point A3. Figure3.e1: The time series of the solution of system (2) approaches asymptotically to the positive equilibrium point A4= (0.99043, 0.6535, 0.6304) in the interior of 𝑅+ 3 . For the data in (5.1) with r5 = 0.01. Figure3.e2: The time series of the solution of system (2) approaches asymptotically to the positive equilibrium point A1= (1, 0, 0) in the interior of 𝑅+ 3 . For the data in (5.1) with r5 = 0.1. Figure3.e3 time series of the solution of system (2) approaches asymptotically to point A3= (0.6667, 0, 3.3333). The system (2) still approaches asymptotically to the equilibrium point A2 despite changing the parameter r6 which represents the conversion rate from the prey to the generalist predator in the range 0.0001 ≀ r6 < 0.015 this causes extinction in the prey, however in additional for 0.015≀ r6<0.1 approaches asymptotically to the positive equilibrium point A4, as shown in Figure4. g2, for typical value r6=0.08. Figure4.g1: The time series of the solution of system (2) approaches asymptotically to the positive equilibrium point A2= (0.3392, 13.2187, 0) in the interior of 𝑅+ 3 . For the data in (5.1) with r6 = 0.013. Figure4.g2: The time series of the solution of system (2) approaches asymptotically to the positive equilibrium point A4= (0.7311, 4.1682, 0.6044) in the interior of 𝑅+ 3 . For the data in (5.1) with r6 = 0.08 6. Discussion Starting with the hypothetical set of data provided by eq. (5.1), system (2) has been numerically solved for several sets of initial points and various sets of parameters, and the following observations are obtained: 1. When approaching globally stable locations via Int. 𝑅+ 3 techniques system (2) only has two types of attractors. The system (2) approaches asymptotically to the globally stable positive point𝐴4 = ( 0.4544, 0.6567 , 2.1722 )for the set parameter value specified in (5.1). 2. The positive equilibrium point 𝐴4 being approached by the solution of system (2) as the assault rate on a victim from the specialized predator r1 increases to 0.15 while maintaining the 0 1 2 3 4 5 6 7 8 9 10 x 10 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Time Po pu lat ion (7.1) h i j 0 1 2 3 4 5 6 x 10 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time Po pu lat ion (14) h i j 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 5 0 0.5 1 1.5 2 2.5 3 Time Po pu lat ion (10.1) h i j 0 0.5 1 1.5 2 2.5 3 x 10 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time Po pu lat ion (8.2) h i j 0 1 2 3 4 5 6 7 8 9 10 x 10 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time Po pu lat ion (8.1) h i j IHJPAS. 2025, 38 (1) 382 other parameters as in eq. (5.1), but if r1= 0.59 it can be seen that the solution to system (2) approaches the equilibrium point A2 asymptotically; r1= 0.595 is a bifurcation point. 3. The trajectory changed from the axial point A1 to the equilibrium point A2 in the range 0.0001 ≀ r3 < 0.006, and from the equilibrium point A2 to the equilibrium point A4 in the range 0.006 ≀ r3 < 0.0125 consequently the bifurcation points for the parameter r3 are at π‘Ÿ3= 0.006 and π‘Ÿ3 = 0.01 respectively. 4. As the attack rate of generalist predator on specialist predator π‘Ÿ4 keeping the rest of parameters as in eq.(5.1), the solution of system (2) approaches to the positive point 𝐴4, if π‘Ÿ4 =1.5, it is shown that the solution of system (2 ) approaches asymptotically to the equilibrium point 𝐴3, indicating that π‘Ÿ4=1.5 is a bifurcation point. 5. The natural death rate of specialist predator r5 the solution of system (2) advances from the positive equilibrium point A4 to the axial equilibrium point A1with 0.03 < r5 < 0.15 keeping the other parameters as in eq.(5.1), and from the axial equilibrim point A1 to the equilibrium point A3 in the ring 0.15≀ r5 <0.6, the parameter π‘Ÿ5 = 0.045 is bifurcation point. 6. As the conversion9 rate of prey to the generalist predator π‘Ÿ6 decreasing, and keeping the rest parameters values as in eq. (5.1 ) the solution of system (2) approaches the equilibrium point A4, while for 0.0001 ≀ r6 < 0.015 , the generalist predator population revives and then the trajectory changed from the point A2 to the positive equilibrium point A4, while for 0.015 ≀ r6 < 0.1 thus, the parameter π‘Ÿ7 = 0.015 is a bifurcation point. 7. In light of the conversion rate of specialist predator to the generalist predator π‘Ÿ7 decreasing to 0.0002 and keeping the rest parameters values as in eq. (5.1) the solution of system (2) approaches the equilibrium point A2, while for 0.00001 ≀ π‘Ÿ7 < 0.005 the generalist predator population revives and when 0.005≀ π‘Ÿ7 <0.2, the trajectory changed from the point A2 to the positive equilibrium point A4, the parameter π‘Ÿ7 =0.005 is a bifurcation point. 7. Conclusion This study used an ecological mathematical model that includes a predator-prey model and a food web, as well as a population of prey and a population of specialized predators as refuges. Additionally, this model includes linear types of functional reactions for the predation of creatures that were not protected. 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 Sciences. Conflict of Interest There are no conflicts of interest. Funding There is no funding for the article. References 1. Mondal N, Barman D, Alam S.A.M. Impact of adult predator incited fear in a stage structured prey– predator model. Environ Dev Sustain. 2021;23(6):9280-9307. IHJPAS. 2025, 38 (1) 383 2. Xie Y, Lu J, Li Y. Stability and bifurcation of a delayed generalized fractional-order prey–predator model with interspecific competition. Appl Math Comput. 2019. 3. Wang Z, et al. Stability and bifurcation of a delayed generalized fractional-order prey–predator model with interspecific competition. Appl Math Comput. 2019;347. 4. Naji RK, Majeed SJ. A prey – predator model with a refuge –stage structure prey population. Int J Differ Equ. 2016; 2016:2010464:10-25. 5. Guckenheimer J. Dynamical systems and bifurcations of vector fields. Appl Math Sci. 1986. 6. Kadhim ZJ, Majeed AA, Naji RK. The bifurcation analysis of a stage-structured prey food web model with refuge. Iraq J Sci. 2016;Special Issue, Part A:139-155. 7. Alabacy ZKH, Majeed AA. The local bifurcation analysis of two preys stage structured predator model with anti-predator behavior. J Phys Conf Ser. 2022;012061. 8. Kafi EM, Majeed AA. The local bifurcation of an eco-epidemiological model in the presence of stage-structured with refuge. Iraqi J Sci. 2020:2087-2105. 9. Wiggins S. Introduction to applied nonlinear dynamical systems and chaos. Springer-Verlag, New York. 1990. 10. Mortoja SG, Panja P, Mondal SK. Dynamics of a predator-prey model with stage-structure on both species and anti-predator behavior. Inform Med Unlocked. 2018;10:50-57. 11. Perko L. Differential equations and dynamical systems. Springer Science & Business Media. 2013. 12. Carr JCW, Hale J. Abelian integrals and bifurcation theory. J Differ Equations. 1985;59:413-436. 13. Hallam TG, De Luna JT. Effects of toxicants on populations: a qualitative approach III. Environmental and food chain pathways. Academic Press Inc (London) Ltd. 1984. 14. Hastings A, Powell T. Chaos in a three-species food chain. Ecology. 1991;72(3):896-903. 15. Molla H, Sarwardi S, Haque M. Dynamics of adding variable prey refuge and an Allee effect to a predator–prey model. Alexandria Eng J. 2021. 16. Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Anim Ecol. 1975;44:331-340. 17. Chakraborty K, Jana S, Kar TK. Global dynamics and bifurcation in a stage structured prey-predator fishery model with harvesting. Appl Math Comput. 2012;218(18):9271-9290. 18. Latha HR, Rama Prasath A. Chaos based dimensional logistic map for image security. J Crit Rev. 2020;7(15). 19. Chen LJ, et al. Qualitative analysis of predator prey models with Holling type II functional response incorporating a constant prey refuge. Nonlinear Anal Real World Appl. 2010;11:246-252. 20. Abdulkadhim MM, Mohsen AA, Al-Husseiny HF. Stability analysis and bifurcation for a bacterial meningitis spreading with stage structure: mathematical modeling. Iraqi J Sci. 2023;21/5. 21. Mohsen AA, Aaid IA. Stability of a prey-predator model with SIS epidemic disease in predator involving Holling type II functional response. IOSR J Math. 2015;11(2):38-53.