EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6026 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Comparative Study on the Impact of Infection on the Stability of Eco-epidemiological Models Kawa Ahmad Hassan 1 Department of Mathematics, College of Science, University of Sulaimani, Sulaymaniyah 46001, Iraq Abstract. This study investigates the stability of eco-epidemiological systems by analyzing a de- terministic ordinary differential equation (ODE) model of predator-prey interactions with disease transmission. The model incorporates both Beddington-DeAngelis and Holling Type IV functional responses. Addressing a critical gap in understanding coexistence dynamics under disease pres- sure, we employ Lyapunov stability theory and Routh-Hurwitz criteria to demonstrate that fatal prey infection universally precludes three-species coexistence, regardless of functional response type. Our results reveal that system dynamics exhibit sensitive dependence on infection and pre- dation rates, showing heightened extinction risks due to prey defense mechanisms. These findings advance theoretical ecology by quantifying disease-predation tradeoffs and provide predictive sta- bility thresholds with direct implications for ecosystem management and conservation strategies in disease-affected environments. 2020 Mathematics Subject Classifications: 92D30, 92D25, 34D20, 37N25, 92D40 Key Words and Phrases: Eco-epidemiological Model, Prey-predator, Local Stability, Beddington- DeAngelis, Holling Type IV 1. Introduction Ecology and epidemiology, though distinct, share several common features and are increasingly integrated into the field of eco-epidemiology, which addresses both ecological and epidemiological dimensions, highlighting the profound influence of disease on ecologi- cal systems from both mathematical and ecological standpoints [1], [2]. The introduction of a parasite into a population of hosts and its subsequent behavior are intricately shaped by interactions with other community members, particularly predators. The intensity of predation plays a critical role in determining community structure and ecosystem charac- teristics. Within host-parasite systems, predation can have a substantial impact, modify- ing the population dynamics of both hosts and parasites, potentially acting as a barrier to parasite establishment [3], [4]. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6026 Email address: kawa.hassan@univsul.edu.iq (K. A. Hassan) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 2 of 20 In theoretical investigations of host-parasite-predator interactions, predator behavior is frequently modeled in a simplified manner. A key component in these models is the functional response, defined as the rate at which a predator consumes prey per unit of time [4]. Recent advances in discrete-time modeling [5] and conformable derivative ap- proaches [6] have expanded our understanding of these dynamics. This response captures characteristics of both predator and prey behavior, shaped by factors such as prey es- cape mechanisms, habitat complexity, and handling time. The present study introduces a host-parasite-predator model to explore the influence of infection rate, attack rate, and the nature of the predator functional response on the system’s behavior [7]. This research contributes original insights through an eco-epidemiological model that analyzes dual predation (susceptible versus infected prey) and disease transmission to predators, while comparing Beddington-DeAngelis and Holling Type IV functional re- sponses. Our analysis identifies critical stability thresholds that reveal how predator be- havior mediates disease-driven ecosystem collapse, advancing both theoretical ecology and conservation applications. Venturino (1995) examined SIS models in which disease transmission occurs within the prey population, incorporating logistic growth for both prey and predator populations, with predators exclusively feeding on infected prey. Hethcote et al. (2004) explores a predator-prey model modification to include logistic growth in the prey population and an SIS disease framework, under the assumption that infected prey is more vulnerable to predation and serves as a resource for predator population growth [2], [8]. However, these models possess certain limitations, as predators generally consume both infected and susceptible prey. Feeding on infected prey can adversely affect predator populations, often resulting in negative growth. Whenever a parasite causes infection in the prey population, the infection can transmit to the predator population via their interactions, potentially causing the collapse of both populations when the infection is lethal. Our work introduces and analyzes a mathematical model incorporating these eco-epidemiological features. Investigating and comparing the model’s dynamics to identify critical system parameters and their ranges is our objective, enabling the prediction of various theoretical outcomes arising from the interplay among susceptible prey, infected prey, and predators. In this study, two functional responses are employed. The first is the Beddington- DeAngelis functional response, an extension of the Holling Type II response, which ac- counts for mutual interference among predators. Our model describes how the rate of prey consumption by predators is influenced not only by prey density but also by the density of the predators themselves. The key feature is the inclusion of an additional term that accounts for the reduction in predation efficiency due to interactions among predators, such as competition or interference. This functional response is particularly useful in eco- logical modeling as it provides a more realistic representation of predator-prey dynamics, especially in environments where predator density significantly influences predation rates [9], [10], [11]. The following form expresses the Beddington-DeAngelis functional response: φ(S, P ) = αS S + bP + c (1) K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 3 of 20 The second functional response utilized is the Holling Type IV, a model employed in ecology to characterize how a predator’s rate of prey consumption varies in response to changes in prey density. Unlike the simpler Holling Types I, II, and III, which describe linear, hyperbolic, and sigmoidal responses respectively, Holling Type IV incorporates a more complex scenario where the predator’s consumption rate initially increases with prey density, then decreases after reaching a peak. This decrease can be due to factors like predator satiation or increased handling time as prey becomes more abundant [12], [13]. This type of functional response is particularly useful in modeling situations where predators experience diminishing returns in prey capture efficiency at high prey densities, reflecting more realistic ecological interactions. The generalized Holling Type IV functional response [14], [15], [16] is given as: φ(S) = αS 1 + a1S + a2S2 , a2 > 0, a1 > −2 √ a2, α > 0. (2) The Holling Type IV response function φ(S), also referred to as the Monod-Haldane response function [14], is positive for S > 0 and increases within the interval S ∈ [0, 1/ √ a2]. The function reaches its maximum value at S = 1/ √ a2, decreases for S > 1/ √ a2, and approaches zero as S tends toward infinity. The phenomenon, when prey populations grow large enough, they become more capable of avoiding, protecting, or defending themselves against predators, is captured in this response function [17], [18]. To formulate our model, which extends classical predator-prey frameworks by integrat- ing disease dynamics in prey, we consider a system where predators preferentially consume infected individuals due to their reduced escape ability. The Beddington-DeAngelis func- tional response accounts for predator interference, while the Holling Type IV functional response captures prey defense mechanisms at high densities. This formulation aligns with empirical observations where disease alters predation efficiency and justifies the inclusion of biologically meaningful terms such as: λIS (disease transmission), γIP (predation on infected prey), we depended on the model used in [19], based on the assumptions in [19] and the functional responses we use in this work as discussed in Section 2 and 3. 2. First model analysis In this section, the boundedness, finding equilibria, and their stability are described with respect to the Beddington-DeAngelis functional response. The eco-epidemiological model incorporating the Beddington-DeAngelis functional response:  dS dT = rS ( 1− S+I K ) − λIS − αSP S+bP+c , dI dT = λIS − γIP − d1I, dP dT = αmSP S+bP+c −mγIP − d2P, (3) K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 4 of 20 where system (3) is examined under the initial conditions S(0) > 0, I(0) > 0, and P (0) > 0. Here, S represents the susceptible population, I corresponds to the infected popula- tion, and P indicates the predator population at time T . To maintain biological validity, all parameters in the model are positive, and their specific biological interpretations are detailed in Table 1. Table 1: Biological interpretations of all parameters. Parameter Biological Meaning r Intrinsic growth rate K Environmental carrying capacity λ Infectious contact rate α Maximal relative increase of predation b Magnitude of interference among predators c Half-saturation constant γ Attack rate on infected prey d1 Natural death rate d2 Natural mortality rate of predator m Conversion efficiency The following theorem demonstrates that the linear combination of susceptible prey, infected prey, and predator populations remains below a finite value, indicating the bound- edness of solutions of system (3). Theorem 1. The solution (S(t), I(t), P (t)) is uniformly bounded for initial condition (S0, I0, P0) ∈ R3 +. Proof. Consider the function W (t) = S(t) + I(t) + 1 mP (t), which is defined from R+ 0 → R+ 0 and is differentiable on the interval (0, t). The derivative dW (t) dT along the trajectory of the system (2.1) can be written as: dW (t) dT = rS ( 1− S + I K ) −λIS− αSP S + bP + c +λIS−γIP−d1I+ αSP S + bP + c −γIP−d2 m P. Now, for any a > 0, we have dW (t) dT + aW (t) = rS − rS2 K − rSI K − 2γIP − d1I − d2 m P + aS + aI + a m P ≤ S ( r + a− r K S ) − (d1 − a)I − ( d2 − a m ) P ≤ S ( r + a− rS K ) , where 0 < a < min{d1, d2}. The maximum value of the expression S ( r + a− rS K ) is K(r+a)2 4r . K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 5 of 20 Hence, dW (t) dT + aW (t) ≤ K(r + a)2 4r , ∀t ∈ (0, T ). Let H(t, y) = K(r + a)2 4r − ay, which satisfies the Lipschitz condition. Clearly, dW (t) dT = K(r + a)2 4r − aW (t) = H(t,W (t)), ∀t ∈ (0, T ). Let dX dT = H(t,X) = K(r + a)2 4r − aX and X0 = W (S0, I0, P0). The linear ordinary differential equation has a solution X(t) = M a ( 1− e−at ) +X0e −at. Since X(t) is bounded on (0, T ), by the comparison theorem [20], W (t) ≤ X(t) = M a ( 1− e−at ) +W (S0, I0, P0)e −at, ∀t ∈ (0, TM ). Thus, for t → ∞, we have 0 < W (t) < M a . Hence, all solutions are bounded uniformly on R+ 0 . For non-dimensionless system (3), we apply the transformation: s = S K , i = I K , p = P K , T = t λK , to obtain the following dimensionless system: ds dt = es (1− (s+ i))− si− n1sp s+bp+d , di dt = si− wip− d3i, dp dt = mn1sp s+bp+d −mwpi− d4p, (4) where e = r λK , n1 = α λK , d = c K , w = γ λ , d3 = d1 λK , and d4 = d2 λK . 2.1. Equilibria The system (4) exhibits the following equilibrium points: • E0 = (0, 0, 0) and E1 = (1, 0, 0) always exist. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 6 of 20 • E2 = ( d3, e(1−d3) e+1 , 0 ) exists where d3 < 1. • E3 = (s̄, 0, p̄) exists in the interior of R3 + if there is a positive solution to the following set of equations: e(1− s̄) = n1p̄ s̄+ bp̄+ d , (5) mn1s̄ s̄+ bp̄+ d = d4. (6) From equation (6), we get: p̄ = (mn1 − d4)s̄− dd4 bd4 . (7) By substituting the value of p̄ into equation (5), a set of polynomial equations is obtained: bms̄2 + (mn1 − bm− d4)s̄− dd4 = 0. (8) Equation (8) has a positive root: s̄ = (d4 + bm−mn1) + √ (d4 + bm−mn1)2 + 4bmdd4 2bm . So, E3 = (s̄, 0, p̄) exists in the interior of R3 + if and only if p̄ is given by equation (7) and s̄ satisfies: 0 < dd4 mn1 − d4 < s̄ < 1. (9) • The interior equilibrium point E∗ = (s∗, i∗, p∗) exists if and only if there is a positive solution to the following nonlinear equations: e(1− (s+ i))− i− n1p F = 0, (10) s− wp− d3 = 0, (11) mn1s F −mwi− d4 = 0, (12) where F = s+ bp+ d. From equation (11), we get: p∗ = s∗ − d3 w . Substituting the value of p∗ into equation (12) yields: i∗ = (wmn1 − d4(w + b))s∗ − d4(dw − bd3) mw((w + b)s∗ + (dw − bd3)) . K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 7 of 20 By substituting p∗ and i∗ into equation (10), we get a polynomial equation: P (s∗) = A1s ∗2 +A2s ∗ +A3, where A1 = mwe(w + b), A2 = mwe(dw − bd3) + (e+ 1)(mwn1 − d4(w + b) +mwn1 −mwe(w + b)), A3 = (e+ 1)(−d4(dw − bd3))−mwn1d3 −mwe(dw − bd3). Clearly, the equation P (s∗) = 0 has a positive root s∗ if dw − bd3 > 0. Hence, E∗ = (s∗, i∗, p∗) exists if: d3 < 1, dw − bd3 > 0, and s∗ > max { d3, d4(dw − bd3) mwn1 − wd4 − bd4 } . (13) 2.2. Stability analysis The general variational matrix of the system at (s, i, p) is J = [uij ]3×3, where: u11 = s ( −e+ n1p F 2 ) + e(1− s− i)− i− n1p F , u12 = −s(e+ 1), u13 = −n1s(s+ d) F 2 , u21 = i, u22 = s− wp− d3, u23 = −wi, u31 = mn1p(bp+ d) F 2 , u32 = −mwp, u33 = −mn1bsp F 2 + mn1s F −mwi− d4, and F = s+ bp+ d. Theorem 2. The system (4) is unstable around E0 for all parameter values, as evidenced by the eigenvalues: λ1 = e > 0, λ2 = −d3 < 0, and λ3 = −d4 < 0. Theorem 3. The system (4) is locally asymptotically stable around E1 = (1, 0, 0), if d3 > 1 and mn1 < (1 + d)d4. Proof. The eigenvalues of J(E1) are: λ1 = −e, λ2 = 1− d3, and λ3 = mn1 1 + d − d4. Then E1 is stable where d3 > 1 and mn1 < (1 + d)d4. Additionally, E1 is globally asymptotically stable if condition (9) does not hold. Theorem 4. The system (4) is locally asymptotically stable around E2 = ( d3, e(1−d3) e+1 , 0 ) , if: n1 < d3 + d md3 ( mwe(1− d3) e+ 1 + d4 ) . (14) K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 8 of 20 Proof. The variational matrix about the equilibrium point E2 = ( d3, e(1−d3) e+1 , 0 ) yields: λ23 = mn1d3 d3 + d − mwe(1− d3) e+ 1 − d4, where the eigenvalues λ21 and λ22 are roots of the polynomial equation: λ2 + ed3λ+ ed3(1− d3) = 0, and satisfy the following relations: λ21 + λ22 = −ed3 < 0 and λ21 · λ22 = ed3(1− d3) > 0 (since d3 < 1). Thus, all eigenvalues possess negative real parts provided that condition (14) is satisfied. Theorem 5. The system (4) is locally asymptotically stable around E3 = (s̄, 0, p̄) if one of the following conditions holds: d3 ≥ 1 and d4 m < n1 < e p̄ (s̄+ bp̄+ d)2, (15) d3 < 1 and d4 m ( d wp̄+ d3 + 1 ) < n1 < e p̄ (s̄+ bp̄+ d)2. (16) Proof. By substituting E3 in the variational matrix, we obtain J(E3) = [aij ]3×3, where: a11 = s̄ ( n1p̄ (s̄+ bp̄+ d)2 − e ) , a12 = −(e+ 1)s̄ (< 0), a13 = −n1s̄(s̄+ d) (s̄+ bp̄+ d)2 (< 0), a22 = s̄− wp̄− d3, a31 = mn1p̄(bp̄+ d) (s̄+ bp̄+ d)2 (> 0), a32 = −mwp̄ (< 0), a33 = −bmn1s̄p̄ (s̄+ bp̄+ d)2 (< 0), a21 = a23 = 0. The characteristic equation associated with this variational matrix can be expressed as: (λ− a22)[(λ− a11)(λ− a33)− a13a31] = 0, or equivalently: (λ− a22)[λ 2 + (a11 + a33)λ+ a11a33 − a13a31] = 0. The eigenvalues of J(E3) satisfy the following relations: λ32 = a22, λ31 + λ33 = a11 + a33, λ31 · λ33 = a11a33 − a13a31. We analyze the following cases: K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 9 of 20 • If d3 ≥ 1, then λ32 = a22 = s̄ − wp̄ − d3 < 0 (since s̄ < 1 from (2.7)). If n1 < e p̄(s̄+ bp̄+ d)2, then a11 < 0, and we achieve λ31 + λ33 < 0 and λ31 · λ33 > 0. Hence, for d3 ≥ 1 with d4 m < n1 < e p̄(s̄ + bp̄ + d)2, E3 exists and is locally asymptotically stable. • If d3 < 1, then λ32 = a22 = s̄ − wp̄ − d3 < 0 where dd4 mn1−d4 < s̄ < wp̄ + d3. We obtain d4 m ( d wp̄+d3 + 1 ) < n1. As in the previous case, if n1 < e p̄(s̄ + bp̄ + d)2, then a11 < 0, and we achieve λ31 + λ33 < 0 and λ31 · λ33 > 0. Hence, for d3 < 1 with d4 m ( d wp̄+d3 + 1 ) < n1 < e p̄(s̄+ bp̄+d)2, E3 exists and is locally asymptotically stable. Theorem 6. For all parameter values, the equilibrium point E∗ is unstable. Proof. By substituting E∗ in the variational matrix, we obtain J(E∗) = [bij ]3×3, where: b11 = s∗ [ n1p ∗ F ∗2 − e ] , b12 = −(e+ 1)s∗ < 0, b13 = −n1s ∗(s∗ + d) F ∗2 < 0, b21 = i∗ > 0, b22 = 0, b23 = −wi∗ < 0, b31 = mn1p ∗(bp∗ + d) F ∗2 > 0, b32 = −mwp∗ < 0, b33 = −bmn1s ∗p∗ F ∗2 < 0, and F ∗ = s∗ + bp∗ + d. The characteristic equation of J(E∗) can be expressed as: λ3 +D1λ 2 +D2λ+D3 = 0, where: D1 = −(b11 + b33), D2 = b11b33 − b13b31 − b23b32 − b12b21, D3 = b11b23b32 − b23b31b12 − b13b21b32 + b12b21b33. By substituting the values of bij into D3, we obtain: D3 = ( n1s ∗p∗ F ∗2 − s∗e ) mw2F ∗p∗−mwn1s ∗(s∗ + d)i∗p∗ F ∗2 + mn1(e+ 1)s∗p∗i∗ F ∗2 (bs∗ − w(bp∗ + d)) . (17) Since s∗ = wp∗ + d3 and simplifying equation (17), we obtain: D3 = −mw2es∗i∗p∗ + mwn1s ∗i∗p∗ (F ∗)2 [wp∗ − (s∗ + d)] + mn1(e+ 1)s∗i∗p∗ (F ∗)2 [b(wp∗ + d3)− wbp∗ − dw] = −mw2es∗i∗p∗ − mwn1es ∗i∗p∗ (F ∗)2 [d3 + d] + mn1(e+ 1) (F ∗)2 s∗i∗p∗ [bd3 − dw] From the existence condition (13) of E∗, bd3−dw < 0. Thus, D3 becomes negative, which violates the Routh-Hurwitz criterion for stability. Therefore, the equilibrium point E∗ is unstable for all parameter values. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 10 of 20 3. Second model analysis In this section, the boundedness, finding equilibria, and their stability are described with respect to the Holling type IV functional response. Eco-Epidemiological Model with Holling Type IV Functional Response The eco-epidemiological model with Holling Type IV functional response is given by: dS dT = rS ( 1− S+I K ) − λIS − αSP 1+a1S+a2S2 , dI dT = λIS − γIP − d1I, dP dT = αmSP 1+a1S+a2S2 −mγIP − d2P, (18) where the system (18) is analyzed under the initial conditions S(0) > 0, I(0) > 0, and P (0) > 0. For this functional response, we use parameters a1 and a2 instead of parameters b1 and b2 used in system (3). The rest of the parameters remain the same. The parameters a1 and a2 represent handling times and inhibitory effects, respectively. Theorem 7. The solution (S(t), I(t), P (t)) is uniformly bounded for initial conditions (S0, I0, P0) ∈ R3 +. Proof. As in Theorem 1. By applying the same transformation as in the previous model, the system (18) can be expressed in the following dimensionless form: ds dt = es(1− (s+ i))− si− n2sp c1+c2s+s2 , di dt = si− wip− d3i, dp dt = mn2sp c1+c2s+s2 −mwpi− d4p, (19) where: n2 = α λa2k2 , c1 = 1 a2k2 , c2 = a1 a2k , e = r λk , w = γ k , d3 = d1 λk , and d4 = d2 λk . 3.1. Equilibria The system (19) exhibits the following equilibrium points: • E′ 0 = (0, 0, 0), E′ 1 = (1, 0, 0) always exist. • E′ 2 = ( d3, e(1−d3) e+1 , 0 ) exists where d3 < 1. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 11 of 20 • E± 3 = (s±3 , 0, p ± 3 ) exist where s±3 = mn2 − c2d4 ± √ B 2d4 , p±3 = em d24 [(1 + c2)d4 −mn2] s ± 3 + c1d4, where B = (mn2 − c2d4) 2 − 4d24c1, under the following condition: d4 m (c2 + 2 √ c1) < n2 < d4 m (1 + c2). (20) • The interior equilibrium point E′ ∗ = (s̄, ī, p̄) exists if and only if there exists a positive solution to the system: e(1− (s+ i))− i− n2p G = 0, (21) s− wp− d3 = 0, (22) mn2s G −mwi− d4 = 0, (23) where G = c1 + c2s+ s2. From equation (21), we have s + i ≤ 1, and from equation (22), we have p̄ = s̄−d3 w . Then, p < 1 if: s̄ > d3. (24) From equations (21) and (23), we obtain: emwG− emwGs̄− (e+ 1)(mn2s̄− d4G)−mwn2p̄ = 0, (25) mn2s̄− d4G = mwGī. (26) By substituting equation (26) into (25), we obtain the cubic polynomial: P (s̄) = A3s̄ 3 −A2s̄ 2 −A1s̄−A0, where: A0 = ( emw2 + (e+ 1)wd4 ) c1 +mwn2d3 > 0, A1 = w(e+ 1)(d4c2 −mn2) +mw(ewc2 − ewc1 − n2), A2 = emw2(1− c2) + (e+ 1)wd4, A3 = emw2. Clearly, P (0) = −A0 < 0 and: P (c1) = c1w(e+1) [mn2 − d4c2 − d4c1 − d4]+mw [ ewc1(c 2 1 + c1c2 − c2 − 1)− n2(c1 − d3) ] . Hence, P (c1) > 0 if: n2 > max { d4(c2 + c1 + 1) m , ewc1(1 + c1 − c1c2 − c21) c1 − d3 } . (27) Since P (s̄) is a continuous function on [0, c1], P (0) < 0, and P (c1) > 0, by the Intermediate Value Theorem, there exists s̄, 0 < s̄ < c1, such that P (s̄) = 0. Thus, E∗ = (s̄, ī, p̄) exists if conditions (24) and (27) hold. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 12 of 20 3.2. Stability Analysis The variational matrix of the system (19) is J ′ = [u′ij ]3×3, where: u′11 = s [ n2p(c2 + 2s) G2 − e ] + e(1− (s+ i))− i− n2p G , u′12 = −s(e+ 1) < 0, u′13 = −n2s G < 0, u′21 = i > 0, u′22 = s− wp− d3, u′23 = −wi < 0, u′31 = mn2p(c1 − s2) G2 , u′32 = −mwp < 0, u′33 = mn2s G −mwi− d4. Theorem 8. The system (19) is unstable around E′ 0 = (0, 0, 0) for all parameter values. Proof. The eigenvalues of E′ 0 are λ1 = e > 0, λ2 = −d3 < 0, and λ3 = −d4 < 0. Thus, E′ 0 is unstable. Theorem 9. The system (19) is locally asymptotically stable around E′ 1 = (1, 0, 0), if d3 > 1 and n2 < d4 m (1 + c1 + c2). Proof. The eigenvalues of E′ 1 are λ1 = −e < 0, λ2 = 1 − d3, and λ3 = mn2 1+c1+c2 − d4. Thus, E′ 1 is stable if d3 > 1 and n2 < d4 m (1 + c1 + c2). Theorem 10. The system (19) is locally asymptotically stable around E′ 2 = ( d3, e(1−d3) e+1 , 0 ) , if and only if: n2 < d4 md3 (c1 + c2d3 + d23). Proof. From the variational matrix at E′ 2, we have the eigenvalue: λ3 = mn2d3 c1 + c2d3 + d23 − emw(1− d3) e+ 1 − d4, which is negative if: n2 < d4 md3 (c1 + c2d3 + d23). The other eigenvalues λ1 and λ2 are the roots of: λ2 + ed3λ+ ed3(1− d3) = 0, with λ1 + λ2 = −ed3 < 0 and λ1λ2 = ed3(1 − d3) > 0. Hence, λ1 and λ2 are either negative real numbers or complex conjugates with negative real parts. By the Routh- Hurwitz criterion, E′ 2 is locally asymptotically stable. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 13 of 20 Theorem 11. The system (19) is locally asymptotically stable around E− 3 = (s−3 , 0, p − 3 ), if one of the following holds: d3 ≥ 1, s−3 < √ c1, and n2 < ep−3 (c2 + 2s−3 ) G2 3 , d3 < 1, s−3 < min{ √ c1, wp − 3 + d3}, and n2 < ep−3 (c2 + 2s−3 ) G2 3 . Proof. The variational matrix around E− 3 = (s−3 , 0, p − 3 ) is: s−3 ( n2p − 3 (c2+2s−3 ) G2 3 − e ) −(e+ 1)s−3 −n2s − 3 G3 0 s−3 − wp−3 − d3 0 mn2p − 3 (c1−(s−3 )2) G2 3 −mwp−3 0  , where G3 = c1 + c2s − 3 + (s−3 ) 2. The eigenvalue in the i-direction is λ2 = s−3 − wp−3 − d3, and for: J ′ 13 = s−3 (n2p − 3 (c2+2s−3 ) G2 3 − e ) −n2s − 3 G3 mn2p − 3 (c1−(s−3 )2) G2 3 0  , we have: tr(J ′ 13) = s−3 ( n2p − 3 (c2 + 2s−3 ) G2 3 − e ) , det(J ′ 13) = mn2 2p − 3 s − 3 (c1 − (s−3 ) 2) G3 3 . For d3 ≥ 1: We have s−3 < 1, so s−3 < d3. Hence, λ2 = s−3 − wp−3 − d3 < 0. Also, det(J ′ 13) > 0 and tr(J ′ 13) < 0 where s−3 < √ c1 and n2 < ep−3 (c2+2s−3 ) G2 3 , respectively. For d3 < 1: When s−3 < min{√c1, wp − 3 + d3}, then λ2 = s−3 − wp−3 − d3 < 0 and det(J ′ 13) > 0. Moreover, if n2 < ep−3 (c2+2s−3 ) G2 3 , then tr(J ′ 13) < 0. Theorem 12. The system (19) is unstable around E+ 3 = (s+3 , 0, p + 3 ). Proof. From the existence condition (20) of E+ 3 = (s+3 , 0, p + 3 ), mn2 > d4(c2 + 2 √ c1), and: s+3 = mn2 − c2d4 + √ (mn2 − c2d4)2 − 4c1d24 2d4 > mn2 − c2d4 2d4 . Thus: s+3 > c2d4 + 2d4 √ c1 − c2d4 2d4 = √ c1. Therefore, s+3 > √ c1, so det(J ′ 13) is negative, which means that λ1 or λ3 is positive. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 14 of 20 Theorem 13. The system (19) is unstable around E′ ∗ = (s̄, ī, p̄) for all parameter values. Proof. From the variational matrix J ′(E′ ∗) = [b′ij ]3×3, the characteristic polynomial corresponding to J ′(E′ ∗) is: λ3 +H1λ 2 +H2λ+H3 = 0, where: H1 = −b′11, H2 = b′12b ′ 21 + b′23b ′ 32 + b′13b ′ 31, H3 = b′11b ′ 23b ′ 32 − b′31b ′ 12b ′ 23 − b′13b ′ 21b ′ 32. And: b′11 = s̄ [ n2p̄(c2 + 2s̄) Ḡ2 − e ] , b′12 = −s̄(e+ 1) < 0, b′13 = −n2s̄ Ḡ < 0, b′21 = ī > 0, b′22 = 0, b′23 = −wī < 0, b′31 = mn2p̄(c1 − s̄2) Ḡ2 , b′32 = −mwp̄ < 0, b′33 = 0, where Ḡ = c1 + c2s̄+ s̄2. Now if b′11 > 0, then H1 = −b′11 is negative, which means one of the eigenvalues is positive. This implies that the equilibrium E′ ∗ = (s̄, ī, p̄) is unstable. If b′11 > 0, then the stability of E′ ∗ = (s̄, ī, p̄) depends on the sign of b′31. For b ′ 31, from the existence condition of E′ ∗. Since 0 < s̄ < c1 and 0 < s + i < 1, it follows that s̄ < 1. Hence, s̄2 < s̄ and s̄2 < c1, which implies: b′31 = (c1 − s̄)2 > 0. Then, H3 is negative. Thus, by the Routh-Hurwitz criterion, one of the eigenvalues of the Jacobian matrix around E′ ∗ is positive, and therefore E′ ∗ = (s̄, ī, p̄) is unstable. 4. Numerical simulation Based on the conditions for existence and the analysis of stability for the equilibrium points, the parameters d3, n1, and n2 are recognized to be important in models (4) and (19). However, a direct comparison of the dynamics between models (3) and (18) in the d3-n1 and d3-n2 parameter planes is not feasible. To address this, we first reformulate all the conditions outlined in the respective theorems using the original system parameters, as outlined in Table 2, and then compare the results in the α-λ parameter plane. For both the Beddington-DeAngelis and Holling type IV functional responses, it is noted that the trivial equilibrium point E0 invariably remains an unstable saddle for all parameter values. If λ < d1 K , E1 is asymptotically stable if α < d2 m ( 1 + c K ) for the Beddington-DeAngelis functional response and α < d2 m ( a2K + 1 K + a1 ) for the Holling type IV functional response. It is important to observe that the net reproductive ratio, R0, is defined as R0 = λS d1 for each case. The number of new infections generated directly from a single infected prey is represented by R0. Whenever R0 < 1, the disease tends to die out, whereas if R0 > 1, K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 15 of 20 Table 2: Stability Conditions for Equilibria Equilibria Beddington-DeAngelis Holling Type IV E0 Unstable Unstable E1 λ < d1 K and α < d2 m ( 1 + c K ) λ < d1 K and α < d2 m ( a2K + 1 K + a1 ) E2 λ > d1 K and α < d1 + cλ md1 ( mrγ(1− d1 λK ) λ(c+K) + d2 ) λ > d1 K and α < d2 m ( λ+ a1d1 + a2d 2 1 λ ) E3 • λ ≤ d1 K and d2 m < α < r p̄(s̄+ bp̄+ c K )2 • λ > d1 K and d2 m ( cλ γKp̄+d1 + 1 ) < α < r p̄(s̄ + bp̄ + c K )2 With respect to E− 3 : • λ ≤ d1 K , s−3 < 1 K √ a2 , and α < ra22K 2p−3 (a1 + 2a2Ks−3 ) (1/K + a1s − 3 + a2K(s−3 ) 2)2 • λ > d1 K , s−3 < min { 1 K √ a2 , 1 λK (λγp−3 + d1), and α < ra22K 2p−3 (a1 + 2a2Ks−3 ) (1/K + a1s − 3 + a2K(s−3 ) 2)2 With respect to E+ 3 : Unstable E∗ Unstable Unstable the disease persists and becomes endemic within the host population. Additionally, it is evident that the net reproductive ratio increases proportionally with the susceptible population, S. Consequently, if the basic reproductive ratio remains below 1 even at the maximum host population level K (i.e., λK d1 < 1 or λ < d1 K ), the infection is unable to spread within the host population. From a biological perspective, this means that when both the infection rate and the rate at which susceptible prey are encountered are sufficiently low, the infected and predator populations cannot persist. As a result, the system reaches an equilibrium where only healthy prey remains. To validate and illustrate the observed results, we provide an example using the fol- lowing fixed parameter values: r = 4, K = 100, b = 0.8, c = 20, a1 = −0.1, a2 = 0.01, γ = 0.06, d1 = 0.2, m = 0.6, d2 = 0.1, (28) while varying only one ecological parameter, α, and one epidemiological parameter λ. Based on the specified parameter values, it is evident that the equilibrium E1 remains K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 16 of 20 stable when α is below 0.2 for the Beddington-DeAngelis functional response and below 0.151667 for the Holling type IV functional response, provided that λ < 0.002. By setting λ = 0.001 and selecting α = 0.15 for Beddington-DeAngelis and α = 0.1 for Holling type IV, we observe that all trajectories originating from the initial condition (50, 50, 50) converge to the equilibrium point where the susceptible prey population S persists as a stable equilibrium. This convergence demonstrates that E1 is globally asymptotically stable for both functional responses. Naturally, in the absence of both predation and infection, the equilibrium density of the susceptible population asymptotically approaches its carrying capacity K. The stability of E1 is further illustrated in Figures 1 and 2 for the Beddington-DeAngelis and Holling type IV functional responses, respectively. 0 25 50 75 100 125 150 175 200 Time 0 20 40 60 80 100 S, I, P S I P 50 60 70 80 90 100S 0 10 20 30 40 50 I 0 10 20 30 40 50 P Figure 1: Stability of E1 with Beddington-DeAngelis functional response. For the specified set of parameter values (28), it is observed that the stability of the equilibrium E2 requires the value of λ to exceed 0.002. By selecting λ = 0.008, it is fur- ther determined that α must remain below 0.75 for the Beddington-DeAngelis functional response and below 0.00633 for the Holling type IV functional response. Consequently, with λ = 0.008 and α set to 0.7 and 0.006 for the Beddington-DeAngelis and Holling type IV responses, respectively, all trajectories originating from the initial condition (50, 50, 50) converge to the predator-free equilibrium E2. At this equilibrium, the susceptible and in- fected prey populations coexist in a stable state. This convergence demonstrates that E2 is globally asymptotically stable for both functional responses. The stability of E2 is visually confirmed in Figures 3 and 4, which depict the dynamics for the Beddington-DeAngelis and Holling type IV functional responses, respectively. In the case of a lower infection rate, it is observed that the stability of the equilibrium E3 requires the value of λ to be less than 0.002. By selecting λ = 0.001, it is further determined that α should lie within the interval (0.1666, 283) for the Beddington-DeAngelis functional response and be less than 866 for the Holling type IV functional response. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 17 of 20 0 25 50 75 100 125 150 175 200 Time 0 20 40 60 80 100 S, I, P S I P 50 60 70 80 90 100S 0 10 20 30 40 50 I 0 10 20 30 40 50 P Figure 2: Stability of E1 with Holling-type IV functional response. 0 25 50 75 100 125 150 175 200 Time 0 10 20 30 40 50 60 70 S, I, P S I P 30 40 50 60 70S 20 30 40 50 60 I 0 5 10 15 20 25 P Figure 3: Stability of E2 with Beddington-DeAngelis functional response. Consequently, with λ = 0.001 and α = 0.3, all trajectories converge to the disease-free equilibrium E3, where the susceptible prey and predator populations coexist in a stable state. This demonstrates that E3 is globally asymptotically stable for both functional responses under the same infection and attack rates. The stability of E3 is illustrated in Figures 5 and 6, which depict the dynamics for the Beddington-DeAngelis and Holling type IV functional responses, respectively. For both functional responses, it is evident that the interior equilibrium, representing the coexistence of all three species, remains unstable across all parameter values. This K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 18 of 20 0 25 50 75 100 125 150 175 200 Time 0 20 40 60 80 S, I, P S I P 20 30 40 50 60 70 80 90S 10 20 30 40 50 60 I 0 10 20 30 40 50 P Figure 4: Stability of E2 with Holling-type IV functional response. 0 50 100 150 200 250 300 350 400 Time 0 25 50 75 100 125 150 175 S, I, P S I P 0 10 20 30 40 50S 0 10 20 30 40 50 I 40 60 80 100 120 140 160 P Figure 5: Stability of E3 with Beddington-DeAngelis functional response. equilibrium acts as a hyperbolic saddle point. 5. Conclusion The instability of the interior equilibrium in an ecological system arises from the detri- mental effects of infected prey on the predator population, which hinders their coexis- tence. Our model demonstrates that infection can disrupt previously stable predator- prey dynamics, while predators have the potential to disrupt otherwise stable interactions between hosts and parasites. These dynamics are critically influenced by the infection K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 19 of 20 0 50 100 150 200 250 300 350 400 Time 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 S, I, P S I P 0 1 2 3 4 5 6 7 8S 0 1 2 3 4 5 I 6 8 10 12 14 16 18 20 P Figure 6: Stability of E− 3 with Holling-type IV functional response. rate and the predator’s attack rate on susceptible prey. The analysis demonstrates that achieving coexistence in host-parasite-predator systems, where prey is affected by a lethal disease, is not feasible. Instead, only specific configurations such as healthy, disease-free, predator-free, or oscillating disease-free systems can be attained by precisely controlling these two critical parameters. Additionally, the biological control paradox is not inher- ently present in eco-epidemiological models, and our framework can be extended to various host-parasite-predator systems by incorporating different functional responses. References [1] K. P. Hadeler and H. I. Freedman. Predator-prey populations with parasitic infection. Journal of Mathematical Biology, 27(6):609–631, 1989. [2] H. W. Hethcote, W. Wang, L. Han, and Z. Ma. A predator–prey model with infected prey. Theoretical Population Biology, 66(3):259–268, 2004. [3] Julian Heidecke, Andrea Lavarello Schettini, and Joacim Rocklöv. West nile virus eco-epidemiology and climate change. PLOS Climate, 2(5):1–26, 05 2023. [4] O. D. Salomón and M. G. Quintana. Eco-epidemiological studies to develop integrated vector surveillance of leishmaniasis vectors in the americas. One Health Implement. Res., 2(2):45–55, 2022. [5] Juan F. Navarro M. Berkal. Qualitative behavior of a two-dimensional discrete-time prey-predator model. Computational and Mathematical Methods, 3(6):e1193, 2021. [6] M. Berkal M. B. Almatrafi. Stability and bifurcation analysis of predator-prey model with Allee effect using conformable derivatives. Journal of Mathematics and Com- puter Science, 36(3):299–316, 2025. [7] D. March, E. Susser, A. Cheskis Gelman, and M. Charles Gelman Professor. The eco- in eco-epidemiology. International Journal of Epidemiology, 35(6):1379–1383, 2006. K. A. Hassan / Eur. J. Pure Appl. Math, 18 (2) (2025), 6026 20 of 20 [8] C. Tannoia, E. Torre, and E. Venturino. An incubating diseased-predator ecoepidemic model. Journal of Biological Physics, 38(4):705–720, 2012. [9] S. Zhang and L. Chen. A study of predator–prey models with the beddington– deangelis functional response and impulsive effect. Chaos, Solitons & Fractals, 27(1):237–248, 2006. [10] Y. Shao and W. Kong. A predator–prey model with beddington–deangelis functional response and multiple delays in deterministic and stochastic environments. Mathe- matics, 10(18):3378, 2022. [11] H. H. Rahman and K. A. Hassan. Fractional order system dynamical behaviors with beddington-deangelis functional response. Passer Journal of Basic and Applied Sciences, 4(2):170–179, 2022. [12] B. Wang and X. Li. Modeling and dynamical analysis of a fractional-order predator– prey system with anti-predator behavior and a holling type iv functional response. Fractal and Fractional, 7(10):722, 2023. [13] D. B. Prakash and D. K. K. Vamsi. Stochastic time-optimal control and sensitivity studies for additional food provided prey-predator systems involving holling type- iv functional response. Frontiers in Applied Mathematics and Statistics, 9:1122107, 2023. [14] J. F. Andrews. A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates. Biotechnology and Bioengineering, X:707–723, 1968. [15] R. M. Etoua and C. Rousseau. Bifurcation analysis of a generalized gause model with prey harvesting and a generalized holling response function of type iii. Journal of Differential Equations, 249(9):2316–2356, 2010. [16] A. Arsie, C. Kottegoda, and C. Shan. A predator-prey system with generalized holling type iv functional response and allee effects in prey. Journal of Differential Equations, 309:704–740, 2022. [17] H. Freedman and G. Wolkowicz. Predator-prey systems with group defence: The paradox of enrichment revisited. Bulletin of Mathematical Biology, 48(5–6):493–508, 1986. [18] W. A. Foster and J. E. Treherne. Evidence for the dilution effect in the selfish herd from fish predation on a marine insect. Nature, 293(5832):466–467, 1981. [19] P.K. Roy N.Bairagi and J. Chattopadhyay. Role of infection on the stability of a predator-prey system, 2007. [20] G. Birkhoff and G.-C. Rota. Ordinary Differential Equations. Wiley, 4th edition, 1989.