EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5397 ISSN 1307-5543 – ejpam.com Published by New York Business Global Dynamical Analysis of a Three Species Food-Web Model With Extended Holling Type II Functional Response Bakhan Bahman Kamal1, Arkan Nawzad Mustafa1,∗ 1 Department of Mathematics, College of Education, University of Sulaimani, Sulaymaniyah 46001, Iraq Abstract. This paper deals with dynamical analysis of a food-web including three logistically growing interaction species, prey, intermediate predators and apex predators. The intermediate predator species predate the prey species according to the Holling type-II functional response, while the apex predator species predate both prey species and intermediate predator species ac- cording to extended Holling type II functional response for two prey species. Firstly, We conduct a thorough analytical examination of the system, demonstrating the positivity and boundedness of the solutions criteria for the persistence of the model are founded, four biologically possible steady states are determined and the local as well as the dynamics around the model’ steady states based on the parameters are also investigated. The occurrence of Hopf-bifurcation of thee model near all steady states, are discussed. Finally, with the help of MATLAB program, it is performed numerical simulations to support the evidence of our analytical results. 2020 Mathematics Subject Classifications: 37N25, 34D20, 37G99, 92D25, 92D40 Key Words and Phrases: Food web, functional response, permanence, stability analysis, Hopf-bifurcation 1. Introduction Mathematical models for the dynamics of interaction between two species of prey and predator, has been derived by many mathematician author [2, 3, 6, 9]. Food web models are important conceptual tool for illustrating the feeding relationships among species within a community, revealing species interactions and community structure, and understanding the dynamics of energy transfer in an ecosystem, Therefore two-species model has been extended to the three-species model by many authors [5, 8, 10, 18]. The important element to represent the dynamics relationship between predator pop- ulation and prey population is the functional response for the predator, which defined as ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5397 Email addresses: bakhan.kamal@univsul.edu.iq (B. B. Kamal), arkan.mustafa@univsul.edu.iq (A. N. Mustafa) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 2 of 21 the number of consumed prey per predator per unit time [9]. C.S. Holling [12] identi- fied three types of functional response Type I, Type II and Type- III. The most useful functional response is the Holling type II functional response, which is characterized by decelerating intake rate [12]. Many authors used this type of functional response for mod- eling the dynamics of interactions between predator and prey [4, 14]. Jha and Ghorai [14] proposed a prey-predator model with selective harvesting between the species using a Holling-type functional response. Hastings [11, 16] has studied the chaotic behavior of an ecology model including the three species. Khan et al. [15] investigated bifurcation analysis of a three-species discrete time. Sk et al. [22] explored the impact of the fear of predators in prey and shelter in a three-species food chain model with delays in hunt- ing cooperation. The authors [20, 21] have investigated the dynamic behavior of a three species system with a scavenger. Diana et al. [7] investigated the three species model’s dynamic behavior with logistic growth in which disease was included. The behavior of a three-species model with time delay and noise was analyzed stochastically by Danane and Torres [13]. Numerous writers have examined the three species model’s dynamics in the presence of both logistic and non-logistic growth as of now. It has been explored the fear effect and stability analysis of the food chain model [1, 17, 19]. In nature, there are many predator species that consume more than one species of prey as well as intermediate predator. For example, lions usually predate a number of large land-based animals, such as antelopes, buffaloes, crocodiles, giraffes, pigs, zebra, wild dogs and wildebeest. In 2022, the Holling type II functional response is extended to more than one prey species [9]. Therefore, this current work is considered and studying the following food web system with extended Holling type II functional response to two preys. dX dt = X [ r1 ( 1− X K ) − βY 1+βTX − α1Z 1+α1T1X+α2T2Y ] , dY dt = Y [ r2 ( 1− Y X ) − α2Z 1+α1T1X+α2T2Y ] , dZ dt = r3Z ( 1− Z X+Y ) , (1) where X(t), Y (t) and Z(t) represent the individual numbers of the prey, intermedi- ate predators and apex predators, respectively. The parameters are positive and their descriptions are given in Table. System (1) is derived on the following assumptions (i) All the species grow logistically. (ii) The apex predators predate both prey and intermediate predators according to ex- tended Holling type II functional response [8]. (iii) The prey population predated by intermediate predators according to Holling type II functional response. The paper organized as follows: in the next section some preliminaries on the model solution property, which are also needed in this work, are given. In the third section, the existence conditions of all feasible and possible steady stat points of system (1) are B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 3 of 21 found and their local as well as globally stability are investigated. In section five, the Hopfbifurcation near to each steady state point is studied. In section six, system (1) numerically solved to observe the impact of parameters and confirm the analytical results in this work. Finally in section eight. a brief conclusion on the total work is given. Where, all other parameter description in Table 1. Parameters Description r1, r2, r3 Intrinsic growth rate for prey, intermediate predators and apex predators, respectively K Carrying capacity for prey species β by intermediate predators’ predation rate T the predator’s average handling time of intermediate predators α1, α2 predator’s search efficiency of prey, intermediate predators, respectively T1, T2 predator’s average handling time of prey, intermediate predators, respectively Table 1: Parameter description of system (1). 2. Preliminaries and permanence In this section, it is proved some lemma on the model solutions, the lemmas are also needed to obtain the result in this paper. A permanent ecological model, imply that all the species in an ecosystem continues to exist, therefore, the criteria that make the model permanent is important, so in this section, the definition of permanence is reviewed and then it is proved that system (1) is permanent under certain conditions. Lemma 2.1. System (1), exhibits a unique solution that is non-negative and satisfy the following inequalities: lim t→∞ sup(X(t)) ≤ K (2) lim t→∞ sup(Y (t)) ≤ K (3) lim t→∞ sup(Z(t)) ≤ 2 K (4) Proof. Right side of system (1), are continuous and has partial derivatives on the space R3, and hence, system (1) satisfies the Lipschitzian condition. Therefore, by uniqueness Theorem, it has unique solution. Further, the time derivative of X,Y and Z are zero in Y Z - Plane, XZ - Plane and XY - plane, respectively. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 4 of 21 Therefore, if the solution initiates at a non-negative point, then the component X,Y and Z of the solution points of system (1), cannot cross any coordinates of the solution points. Hence components X,Y and Z of solution points is always non negative. The first equation in system (1) gives that dX dt ≤ r1X ( 1− X K ) . Solving above differential inequality, it gets lim t→∞ sup(X(t)) ≤ K. Apply above inequality at the second equation of system (1), it gets dY dt ≤ r2Y ( 1− X K ) . Again Solving above differential inequality, it gets lim t→∞ sup(Y (t)) ≤ K. So, inequalities (2) and (3) are guaranteed. Apply inequalities (2, 3) in the third equation of system (1) , it gets dZ dt ≤ r1Z ( 1− Z 2K ) So, lim t→∞ sup(Z(t)) ≤ 2K. This completes the proof. Lemma 2.2. (i) limt→∞ infX(t) ≥ K1, if r1 > (β + 2α1)K (5) (ii) limt→∞ inf Y (t) ≥ K2, if Inq.(5) and the following inequality are satisfied r2 > 2α2 (1 + α1T1K1) K (6) (iii) limt→∞ inf Z(t) ≥ K1 +K2, if both conditions (5) and (6) are provided. where, K1 and K2 are positive number to be determined in the proof. Proof. Suppose t approaches infinity, then B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 5 of 21 (i) Apply inequalities (2)–(4) in Lemma 2.1, at the first equation of system (1), it gets dX(t) dt ≥ X [ r1 ( 1− X K ) − (β + 2α1)K ] . Therefore under condition (5), it gets lim t→∞ infX(t) ≥ . [r1 − (β + 2α1)K]K r1 = K1 > 0 (ii) From condition(4) in Lemma 2.1 and inequality (7), it follows dY dt ≥ Y [ r2 ( 1− Y K1 ) − 2α2K 1 + α1T1K1 ] , so under condition (6), it is concluded that: lim t→∞ inf Y (t) ≥ [ 1− 2α2K r2 (1 + α1T1K1) ] K1 = K2 > 0 (iii) From (7) and (8), it follows dZ dt ≥ r3Z ( 1− Z K1 +K2 ) , and hence lim t→∞ inf Z(t) ≥ K1 +K2. Definition 2.1. [5] System (1) is said to be permanent if there exist positive constants a and b such that b ≥ max { lim t→∞ sup(X(t)), lim t→∞ sup(Y (t)), lim t→∞ sup(Z(t)) } ≥ min { lim t→∞ inf(X(t)), lim t→∞ inf(Y (t)), lim t→∞ inf(Z(t)) } ≥ a. Theorem 2.1. If the following conditions (5) and (6) in Lemma 2.2 are provided, then system (1), is permanent. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 6 of 21 Proof. From Lemma 2.1, it is obtained that min { lim t→∞ inf(X(t)), lim t→∞ inf(Y (t)), lim t→∞ inf(Z(t)) } ≥ min {K1,K2} . From Lemma 2.2, it follows 2K ≥ max { lim t→∞ sup(X(t)), lim t→∞ sup(Y (t)), lim t→∞ sup(Z(t)) } . The proof is completed. 3. Stability analysis This section including three subsections. In the first subsection, the existence condi- tions of all feasible and possible steady stat points of system (1) are determined and their local and globally stability are investigated in the second and third subsection, respectively. 3.1. Existence of steady states System (1) has at most four steady states, they are the only prey existence steady state S1 = (K, 0, 0), apex predatorfree steady state S2 = (S̄, S̄, 0), intermediate predators-free steady state S3 = (S̄, 0, S̄) and coexistence steady state is S4 = (X∗, Y ∗, Z∗) always exist, where S̄ = r1βTK − r1 − βK + √ (r1βTK − r1 − βK)2 + 4r21βTK 2r1βT , S̄ = r1α1T1K − r1 − α1K + √ (r1α1T1K − r1 − α1K)2 + 4r21α1T1K 2r1α1T1 , Z∗ = X∗ + Y ∗, Y ∗ = a1X ∗ + a2X ∗2 + a3X ∗3 and X∗ is root for the following function F (X) = r1 ( 1− X K ) − βY (X) 1 + βTX − α1(X + Y (X)) 1 + α1T1X + α2T2Y (X) , (10) with a1 = r1α2βT r2α1K , a2 = r1α2 r2α1K − r1α2βT r2α1 and a3 = r2β r2α1 + 1 − r1α2 r2α1 and Y (X) = a1X + a2X 2 + a3X 3. Note the steady states S1, S2 and S3 are always uniquely exists. But the existence criteria for the coexistence steady state S4 is determined in the following Theorem Theorem 3.1. The coexistence steady state S4 = (X∗, Y ∗, Z∗) exist uniquely, if ai ≥ 0, i = 2, 3 (11) B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 7 of 21 Proof. From (10), it gets F (0) = r1 > 0 and F (K) = − βY (K) 1 + βTK − α1(K + Y (K)) 1 + α1T1K + α2T2Y (K) < 0. Intermediate value Theorem, guarantees that F (X) has a root, namely X∗ in (0,K). Further dF (X) dX < 0, for all x ≥ 0, due to condition (11), that is F (X) is decreasing on the positive real line, and hence X∗ is unique root of F (X), and this completes the proof. 3.2. Local stability Here, local asymptotical stability (LAS) of all the steady states of system (1) is studied. Suppose λiX , λiY and λiZ represent the eigenvalues of the Variational matrix at Si, in the X-, Y -, and Z-directions, respectively; i = 1, 2, 3, 4. Then, (i) λ1X = −r1, λ0Y = r2 and λ1Z = r3 (ii) λ2Z = r3 and λ2X and λ2X are roots of the equation λ2 +A1λ+B1 = 0, where, A1 = r1 K S̄ + r2 − ( β2T S̄2 ) /(1 + βT S̄)2 and B1 = r2r1 K S̄ + ( r2βS̄ ) /(1 + βT S̄)2 (iii) λ3X , λ3Y and λ3Z are roots of the equation λ3 +A2λ 2 +B2λ+ r3C2 = 0, where, A2 = − (R1 +R5) , B2 = R1R5−R2R4−r3 (R1 +R3 +R5) and C2 = R1R5−R2R4−R3R4+R3R5 with R1 = −r1 K S̄ + ( α2 1T1S̄ 2 ) / ( 1 + α1T1S̄ )2 , R2 = − βS̄ (1 + βT s̄) − α1S̄/ ( 1 + α1T1S̄ ) , R3 = −α2S̄/ ( 1 + α1T1S̄ ) , R4 = −α2S̄/ ( 1 + α1T1S̄ )2 and R5 = r2 − ( α2S̄ ) / ( 1 + α1T1S̄ ) (iv) λ4X , λ4Y and λ4Z are roots of the equation λ3 +A3λ 2 +B3λ+ r3C3 = 0, where A3 = − (E1 + E5) , B2 = E1E5 − E2E4 − r3 (E1 + E3 + E5 + E6) and C3 = E1E5 + E1E6 − E2E4 − E2E6 − E3E4 + E3E5 B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 8 of 21 with E1 = r1 − 2 r1 K X∗ − βY ∗ (1 + βTX∗)2 − α1 (1 + α2T2Y ∗)Z∗ (1 + α1T1X∗ + α2T2Y ∗)2 , E2 = − βx∗ 1 + βTx∗ − α1 (1 + α1T1X ∗)Z∗ (1 + α1T1X∗ + α2T2Y ∗)2 , E3 = − α1x ∗ 1 + α1T1X∗ + α2T2Y ∗ E4 = r2 ( Y ∗ X∗ )2 − α2 (1 + α2T2Y ∗)Z∗ (1 + α1T1X∗ + α2T2Y ∗)2 , E5 = r2 − 2r2 Y ∗ X∗ − α2 (1 + a1T1X ∗)Z∗ (1 + α1T1X∗ + α2T2Y ∗)2 and E6 = − α2Y ∗ 1 + α1T1X∗ + α2T2Y ∗ Therefore, the following Theorem can be derived based on the above argument. Theorem 3.2. (i) The only prey existence steady state and the apex predator-free steady state are unstable. (ii) The intermediate predators-free steady state is LAS if and only if, all the following Routh-Hurwize criteria hold. A2 > 0, C2 > 0 and A2B2 > r3C2 (12) (iii) The coexistence steady state it is LAS if and only if, all the following Routh-Hurwize criteria hold. A3 > 0, C3 > 0 and A3B3 > r3C3 (13) 3.3. Global stability Global stability (or globally asymptotically stable GAS) means that any trajectories finally tend to the attractor of the system, regardless of initial conditions. Therefore, Most of biological systems, especially prey predator system, are needed to be globally stable. Since S1 and S2 are not LAS, so they cannot be GAS. However GAS for both S3 and S4 of the system (1) is established in Theorem 3.3 and Theorem 3.4, respectively. Theorem 3.3. Suppose that condition 5 holds, and then S3 = (S̄, 0, S̄) is GAS, if the following inequality holds: r2 < α2K1 1 + α1T1K + α2T2K (14) r1 + r3 K2 > α1 (1 + α1T1K1) 2 (15) B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 9 of 21 Proof. From Lemma 2.1, limt→∞ supX(t) ≤ K and limt→∞ supY (t) ≤ K. From Lemma 2.2 limt→∞ infX(t) ≥ K1, if condition 5 , holds. Therefore, the third equation f system (1) gives that lim t→∞ infX(t) ≥ K1 and lim t→∞ inf Y (t) ≥ K1. And hence, from the second equation of system (1), it is obtained that dY dt ≤ Y [ r2 − α2K1 1 + α1T1K + α2T2K − r2Y K ] . So, dY dt is negative due to condition 14 , consequently limt→∞ Y (t) = 0. Therefore as time approaches infinity, system (1) reduced to the following subsystem{ dX dt = X [ r1 ( 1− X K ) − α1Z 1+α1T1X ] = F (X,Z) dZ dt = r3Z ( 1− Z X ) = G(X,Z) (16) Consider now the function H(X,Z) = 1/XZ, clearly H is a continuously differentiable function. Further, ∂(HF ) ∂X + ∂(HG) ∂Z = − r1 KZ + α1 (1 + α1T1X)2 − r3 X2 < α1 (1 + α1T1K1) 2 − r1 + r3 K2 . Under the condition (15), it is clear that ∂(HF ) ∂X + ∂(HG) ∂Z does not change sign and is not identically zero. So, by BendixsonDulac criterion, there is no periodic curve in of the XZ-plane. Since (S̄, S̄) represent unique positive equilibrium point of the subsystem (16), so lim t→∞ X(t) = lim t→∞ Z(t) = S̄. This completes the proof. Theorem 3.4. Suppose that S4 = (X∗, Y ∗, Z∗) is exist, then it is GAS, if in addition of conditions (5), (6), the following inequalities hold: [ 2α1α2Z ∗ (T2 + T1) G2(X,Y ) + 2r2Y ∗ XX∗ − 2 ( β + β2TX∗) F2(X) ]2 < 4 [ r1 K − β2TY ∗ F2 (K1) − α2 1T1Z ∗ G2 (K1,K2) ] [ r2 K − α2 2T2Z ∗ G2 (K1,K2) ] (17) [ 2 ( α1 + α2 1T1X ∗ + α1α2T2Y ∗) G2(X,Y ) − 2r3 (X + Y ) ]2 < 2r3 K [ r1 K − β2TY ∗ F2 (K1) − α2 1T1Z ∗ G2 (K1,K2) ] B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 10 of 21 [ 2 ( α2 + α1α2T1X ∗ + α2 2T2Y ∗) G2(X,Y ) − 2r3 (X + Y ) ]2 < 2r3 K [ r2 K − α2 2T2Z ∗ G2 (K1,K2) ] , where, F2(X) = (1 + βTX) (1 + βTX∗) and G2(X,Y ) = (1 + α1T1X + α2T2Y ) (1 + α1T1X ∗ + α2T2Y ∗) (19) Proof. Consider the function L(X,Y, Z) = 2 [ X −X∗ −X∗ ln ( X X∗ )] + 2 [ Y − Y ∗ − Y ∗ ln ( Y Y ∗ )] + 2 [ Z − Z∗ − Z∗ ln ( Y Z∗ )] . Clearly, L(X,Y, Z) ∈ C1 ( R3 +, R ) With L (X∗, Y6, Z6) = 0 and L(X,Y, Z) > 0,for all (X,Y, Z) ∈ R3 + with (X,Y, Z) ̸= (X∗, Y6, Z6). Further, dL dt =− [ r1 K − β2TY ∗ F2(X) − α2 1T1Z ∗ G2(X,Y ) ] (X −X∗)2 − [ r2 X − α2 2T2Z ∗ G2(X,Y ) ] (Y − Y ∗)2 + [ 221α2Z ∗ (T2 + T1) G2(X,Y ) + 2r2Y ∗ XX∗ − 2 ( β + β2TX∗) F2(X) ] (Y − Y ∗) (X −X∗) − [ r1 K − β2TY ∗ F2(X) − α2 1T1Z ∗ G2(X,Y ) ] (X −X∗)2 − r3 (X + Y ) (Z − Z∗)2 − [ 2 ( α1 + α2 1T1X ∗ + α1α2T2Y ∗) G2(X,Y ) − 2r3 (X + Y ) ] (X −X∗) (Z − Z∗) − [ r2 X − α2 2T2Z ∗ G2(X,Y ) ] (Y − Y ∗)2 − r3 (X + Y ) (Z − Z∗)2 − [ 2 ( α2 + α1α2T1X ∗ + α2 2T2Y ∗) G2(X,Y ) − 2r3 (X + Y ) ] (Y − Y ∗) (Z − Z∗) . Condition (5, 6) guaranteeing that: limt→∞ infX(t) ≥ K1 and limt→∞ inf Y (t) ≥ K2. Therfore, B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 11 of 21 dL dt <− [ r1 K − β2TY ∗ F2 (K1) − α2 1T1Z ∗ G2 (K1,K2) ] (X −X∗)2 − [ r2 K − α2 2T2Z ∗ G2 (K1,K2) ] (Y − Y ∗)2 + [ 2α1α2Z ∗ (T2 + T1) G2(X,Y ) + 2r2Y ∗ XX∗ − 2 ( β + β2TX∗) F2(X) ] (Y − Y ∗) (X −X∗) − [ r1 K − β2TY ∗ F2 (K1) − α2 1T1Z ∗ G2 (K1,K2) ] (X −X∗)2 − r3 2K (Z − Z∗)2 − [ 2 ( α1 + α2 1T1X ∗ + α1α2T2Y ∗) G2(XY ) − 2r3 (X + Y ) ] (X −X∗) (Z − Z∗) − [ r2 K − α2 2T2Z ∗ G2 (K1,K2) ] (Y − Y ∗)2 − r3 2K (Z − Z∗)2 − [ 2 ( α2 + α1α2T1X ∗ + α2 2T2Y ∗) G2(X,Y ) − 2r3 (X + Y ) ] (Y − Y ∗) (Z − Z∗) . Provide both inequalities (17-19), it gets dL dt <− [√ r1 K − β2TY ∗ F2 (K1) − α2 1T1Z∗ G2 (K1,K2) + √ r2 K − α2 2T2Z∗ G2 (K1,K2) ]2 − [√ r3 2K − √ r2 K − α2 2T2Z∗ G2 (K1,K2) ]2 − [√ r1 K − β2TY ∗ F2 (K1) − α2 1T1Z∗ G2 (K1,K2) − √ r3 2K ]2 . So, dL dt is negative, and hence L2 is Lyapunov function with respect to S4 = (X∗, Y ∗, Z∗), this proves the Theorem. 4. Hopf-bifurcation Here, the occurrence of Hopf bifurcation in system1 near all steady states, are discussed as follows: From Theorem (3) it is observed that the only prey existence steady state S1 = (K, 0, 0), apex predator-free steady state S2 = (S̄, S̄, 0), are not LAS, Therefore there is no possibility to have a Hopf bifurcation near S1 and S2. The conditions that guaran- tee the occurring of Hopf bifurcation near, intermediate predators-free steady state S3 =( S̄, 0, S̄ ) and coexistence steady state S4 = (X∗, Y ∗, Z∗) are established in Theorem 4.1 and Theorem 4.2, respectively. Theorem 4.1. System (1) exhibits a Hopf bifurcation near intermediate predators-free steady state as the parameter r3 passes through the value. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 12 of 21 r3 = A2 (R1R5 −R2R4) C2 +A2 (R1 +R3 +R5) . (20) If one of the following conditions holds R1 +R5 < 0, R1R5 > R2R4 and R1 +R3 +R5 > −C2 A2 , R1 +R5 < 0, R1R5 < R2R4 and R1+R3+R5 < −C2 A2 , where A2, C2 and Ri; i = 1, 2, 3, 4, 5 are given in previous section. Proof. Recall that the eigenvalues of Variational matrix at S3 = (S̄, 0, S̄) satisfy λ3 +A2λ 2 +B2λ+ r3C2 = 0. (22) The Hopf bifurcation near intermediate predators-free steady state of system (1) occurs if and only if two roots of Eq. (22) have two complex conjugate eigenvalues with the third eigenvalue real and negative such that there exists a constant parameter value, say r3 satisfying: Ree (λ (r3))r3=r3 = 0 and [ dRe e (λ (r3)) dr3 ] r3=r3 ̸= 0 where λ is a complex root of Eq. (22) if r3 = r̄3, then A2B2 = r3C2 and hence equation () can be written ( λ2 +B2 ) (λ+A2) = 0 (23) Clearly, equation (23) has the following three roots: λ1 = i √ B2, λ2 = −i √ B2 and dλ3 = −A2. However, for all values of r3 in the neighborhood of r3, these roots can be written in general as λ1 = a (r3) + ib (r3) , λ2 = a (r3)− ib (r3) and λ3 = −A2 (r3) . Clearly, a (r3) = 0 which means that the first condition of Hopf bifurcation holds. Now, the proof will follows, if we can verify the above second condition (known as transversality condition) of Hopf bifurcation when Ree (λ (r3)) = a (r3). Thus, by substituting λ1 = a (r3) + ib (r3) in (22) and calculating the derivative with respect to the r3, it is obtained that{ D1 (r3) da(r3) dr3 −D2 (r3) db(r3) dr3 +D3 (r3) = 0 D1 (r3) db(r3) dr3 +D2 (r3) da(r3) dr3 +D4 (r3) = 0 (24) where D1 (r3) = 3 [a (r3)] 2 + 2A2 (r3) a (r3) +B2 (r3)− 3 [b (r3)] 2 B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 13 of 21 D2 (r3) = 6a (r3) b (r3) + 2A2 (r3) b (r3) D3 (r3) = a (r3) dB2 (r3) dr3 + C2 D4 (r3) = b (r3) dB2 (r3) dr3 Thus, by solving the linear system (24) for the unknown da(r3) dr3 ,it gets da (r3) dr3 = dRee (λ (r3)) dr3 = −D2 (r3)D4 (r3) +D1 (r3)D3 (r3) [D1 (r3)] 2 + [D2 (r3)] 2 So for r3 = r3, it is easy to verify that[ dRee (λ (r3)) dr3 ] r3=r3 ̸= 0 Also, condition (20) or (21) guarantee that A2 is positive for all values of r3, and hence λ3 = −A2 (r3) < 0. The proof is completed. In similar ways of proving above theorem we can get the following theorem. Theorem 4.2. System (1) exhibits a Hopf bifurcation coexistence steady state as the parameter r3 passes through the value r3 = A3 (E1E5 − E2E4) C3 +A3 (E1 + E3 + E5) (25) If one of the following condition holds E1 + E5 < 0, E1E5 > E2E4 and E1 + E3 + E5 + E6 > −C3 A3 E1 + E5 < 0, E1E5 < E2E4 and E1 + E3 + E5 + E6 < −C3 A3 where A3, C3 and Ei; i = 1, 2, 3, 4, 5, 6 are given in previous section 5. Numerical simulation In order to support the analytical finding in this paper via the phase portrait and time series[22] , some numerical simulations are performed; all the simulations are carried out through runga kutta method of order six method, using MATLAB. First, we assume two set of parameter values as given in (27) and (28). r1 = 0.5, r2 = 0.3, r3 = 0.15, K = 1000, β = 0.001 T = 1, α1 = 0.005, α2 = 0.004, T1 = 1 and T2 = 1 (27) B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 14 of 21 r1 = 1.2, r2 = 0.9, r3 = 0.15, K = 1000, β = 0.001 T = 1, α1 = 0.005, α2 = 0.004, T1 = 1 and T2 = 1 (28) The parameter values in (27) and (28), satisfy the condition for global stability of the intermediate predator-free steady state and coexistence steady state, respectively. Figure 1: Both the time series and the phase portrait shows that model 1 approaches intermediate predator free Steady state, when the parameter values are as given in (27). B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 15 of 21 Figure 2: Both the time series and the phase portrait shows that model 1 approaches coexistence Steady state, when the parameter values are as given in (28) Therefore, Fig. 1 and Fig. 2 confirm the analytical results regarding to global stability of intermediate predator-free steady state and coexistence steady state, respectively. The critical values r3 ≈ 0.11 and r3 ≈ 0.95 for the parameters given in (27) and (28), respectively. First let decrease the value of r3 to 0.1 and fixed another parameter as given in (27), then it is observed that, the dynamics of the system (1) induced a transition from the a stability situation to the state where the prey species and apex predators oscillate peri- odically oscillate periodically as illustrated in Fig.3. Further, for a range of values r3 ∈ [004,0.12], in Fig.4, bifurcation diagram with respect to r3 is drawn. From Fig. 4, the prey and apex predator populations show oscillatory for the further smaller values of r3; This indicates that decreasing r3 may induce a transition from the a stability situation to the state where the populations oscillate periodically. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 16 of 21 Figure 3: Both the time series and the phase portrait show periodic oscillations around intermediate predators- free steady state, when r3 = 0.1 and other parameter values are as given in (27). B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 17 of 21 Figure 4: The figure shows the bifurcation diagrams of system (1), which indicates the emergence of Hopf bifurcations as parameter r3 decreases, and other parameters are fixed as in(27). Again let further decrease the value of r3 to 0.07 and fixed other parameters given in (28), then it is observed that, the dynamics of the system (1) induced a transition from the a stability situation to the state where the population oscillate periodically oscillate periodically as illustrated in Fig.5. Further, for a range of values r3 ∈ [0, 0.2], in Fig.6, bifurcation diagram with respect to r3 is drawn. From Fig. 6, all the species show oscillatory for the further smaller values of r3; This indicates that decreasing r3 may induce a transition from the a stability situation to the state where the populations oscillate periodically. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 18 of 21 Figure 5: Both the time series and the phase portrait shows periodic oscillations around coexistence Steady state, when r3 = 0.07 and other parameter values are as given in (28). Note that Fig. 3 and Fig. 5 confirms the analytical finding in Theorem 6 and Theorem 7, respectively. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 19 of 21 Figure 6: The figure shows the bifurcation diagrams of system (1), which indicates the emergence of Hopf bifurcations as parameter r3 decreases, and other parameters are fixed as in(28). 6. Conclusion A three species food-web is considered in this paper, the species are, prey, interme- diate predators and apex predators, the intermediate predators predating the prey with the Holling type-II functional response, while the apex predators predating both prey and intermediate predators according to extended Holling type II functional response for two prey species. it is proved that under condition (5)and(6), the model is per- manent, it is explored that the model has biologically possible steady states, they are S1 = (K, 0, 0), S2 = (S̄, S̄, 0), S3 = (S̄, 0, S̄) and S4 = (X∗, Y ∗, Z∗), and also it is discov- ered that both S1 = (K, 0, 0), S2 = (S̄, S̄, 0) are unstable, but stability as well as global stability for both S3 = (S̄, 0, S̄) and S4 = (X∗, Y ∗, Z∗) based on the sample parameters are given. And also it is proved that critical values r3 and r3, make the occurrence of Hopf-bifurcation of thee model near S3 = (S̄, 0, S̄) and S4 = (X∗, Y ∗, Z∗), respectively. The model solved numerically by choosing suitable value of the model parameters as given in (27) and (28), the numerical solutions insulated in Fig. 1 and Fig.2, it is explained that both figures confirms analytical results regarding to the stability of both S3 = (S̄, 0, S̄) and S4 = (X∗, Y ∗, Z∗), respectively. Also we have shown that Fig 3, and Fig.4, confirm analytical results regard to Hopf-bifurcation near S3 = (S̄, 0, S̄) and S4 = (X∗, Y ∗, Z∗), respectively. Because they show periodic solution of the model. B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 20 of 21 References [1] Dakhil R. A. and Majeed S. J. Three-species lotka-volterra food chain model with fear effect and hunting cooperation. University of Thi-Qar Journal, 18(1), 2023. [2] I. Al-Darabsah, X. Tang, and Y. Yuan. A prey-predator model with migrations and delays. Discrete Cont. Dyn. Syst.-B, 21:737–761, 2016. [3] G. E. Arif, J. Alebraheem, and W. B. Yahia. Dynamics of predator-prey model under fluctuation rescue effect. Baghdad Sci. J., 20:1742–1750, 2023. [4] B. Divya and K. Kavitha. Dynamical behavior of three species model with holling type functional response in the presence of logistic growth. J. Appl. Sci. Eng., 28(1):101– 108, 2025. [5] S. Djilali and B. Ghanbari. Dynamical behavior of two predators-one prey model with generalized functional response and time-fractional derivative. Adv. Differ. Equations, 2021(1):1–19, 2021. [6] M. F. Elettreby. Two-prey one-predator model. Chaos, Solitons Fractals, 39:2018– 2027, 2009. [7] Diana A. F., Widowati W., and Tjahjana R. H. Stability analysis of lotka-volterra model for three species with disease. AIP Conference Proceedings, 2738(1), 2023. [8] A. G. Farhan, N. Khalaf, and T. Aldhlki. Dynamical behaviors of four-dimensional prey-predator model. Eur. J. Pure Appl. Math., 16(2):899–918, 2023. [9] S. Fatah, A. Mustafa, and A. Shilan. Predator and n-classes-of-prey model incorpo- rating extended holling type ii functional response for n different prey species. AIMS Mathematics, 8:5779–5788, 2022. [10] U. Ghosh, S. Sarkar, and B. Mondal. Study of stability and bifurcation of three species food chain model with nonmonotone functional response. Int. J. Appl. Comput. Math., 7:1–24, 2021. [11] A. Hastings and T. Powell. Chaos in a three-species food chain. Ecology, 72(3):896– 903, 1991. [12] C. S. Holling. Some characteristics of simple types of predation and parasitism. Can. Entomol., 91(7):385–398, 1959. [13] Danane J. and Torres D. F. Three-species predator-prey stochastic delayed model driven by lévy jumps and with cooperation among prey species. Mathematics, 11(7):1595, 2023. [14] P. Jha and S. Ghorai. Stability of prey-predator model with holling type response function and selective harvesting. J. Appl. Comput. Math., 6(3):1–7, 2017. [15] A. Khan, S. Qureshi, and A. Alotaibi. Bifurcation analysis of a three-species discrete- time predator-prey model. Alexandria Eng. J., 61(10):7853–7875, 2022. [16] A. Klebanoff and A. Hastings. Chaos in three-species food chains. J. Math. Biol., 32:427–451, 1994. [17] Faraj B. M., Sabir P. O., Salih D. T. M., and Hilmi H. Comment on the paper by jalal et al. [chaos, solitons and fractals 135 (2020) 109712]. Chaos, Solitons & Fractals, 187:115428, 2024. [18] S. Mishra and R. K. Upadhyay. Exploring the cascading effect of fear on the foraging B. B. Kamal, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (1) (2025), 5397 21 of 21 activities of prey in a three-species agroecosystem. Eur. Phys. J. Plus, 136:1–36, 2021. [19] Soumitra P., Tiwari P. K., Mishra A. K., and Wang H. Fear effect in a three- species food chain model with generalist predator. Mathematical Biosciences and Engineering, 21:1–33, 2023. [20] Khan A. Q. and Kazmi S. S. Bifurcations of a three-species prey-predator system with scavenger. Ain Shams Engineering Journal, 14(11):102514, 2022. [21] Khan A. Q. and Kazmi S. S. Dynamical analysis of a three-species discrete bio- logical system with scavenger. Journal of Computational and Applied Mathematics, 440:115644, 2024. [22] Sk N. Tiwari and P. K. Pal. A delay nonautonomous model for the impacts of fear and refuge in a three-species food chain model with hunting cooperation. Mathematics and Computers in Simulation, 192:136–166, 2022.