East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

ARTICLE

ARTICLE INFO

Volume 7(1), 2026

https:

Bifurcation Analysis of Eco-Epidemiological

Mathematical Model with Saturated

Incidence Rate and General Holling Type

Response Function

//dx.doi.org/10.4314/eajbcs.v7i1.1S

ARTICLE HISTORY

Received: 16 November, 2025

Accepted: 30 April, 2026

Published Online: 10 June, 2026

Solomon Molla Alemu¹, Tesfaye Tefera Mamo², Mohammed Yiha

Dawed^3,∗�

CITATION

Alemu et.al (2026) Bifurcation Analysis

of Eco-Epidemiological Mathematical

Model with Saturated Incidence Rate

and General Holling Type Response

Function. East African Journal of

¹Addis Ababa Science and Technology University, Department of Mathematics, Addis Ababa, Ethiopia,

²Debre Berhan University, Department of Mathematics, Debre Berhan, Ethiopia,

³Hawassa University, Department of Mathematics, Hawassa, Ethiopia

Biophysical and Computational

Sciences Volume 7(1), 2026. .https://dx.

doi.org/10.4314/eajbcs.v7i1.1S.1-17

^∗Corresponding author: mohammedyiha@hu.edu.et

Abstract

OPEN ACCESS

This paper presents a bifurcation analysis of an Eco-epidemiological model with saturated incidence

rate and general Holling-Type functional responses. The model describes a predator–prey system in

which the prey population is infected by a communicable disease, and the predator feeds on both

susceptible and infected individuals. Fundamental properties of the system, including existence and

uniqueness, positivity, and boundedness of solutions, are established to ensure biological feasibility.

Equilibrium points are identified and their stability is examined. The basic reproduction number '

0

This work is licensed under the Creative

Commons open access license (CC

BY-NC 4.0).

is derived to determine threshold conditions for disease persistence. Using Sotomayor’s theorem,

transcritical and Hopf bifurcations are rigorously verified. The results indicate that increasing the

inhibition rate stabilizes the system and promotes coexistence, whereas higher transmission rates

destabilize equilibria and generate sustained oscillations. Numerical simulations and bifurcation

diagrams validate the analytical findings, demonstrating transitions between stable steady states and

periodic dynamics.

East African Journal of Biophysical and

Computational Sciences (EAJBCS) is

already indexed on known databases

like AJOL, DOAJ, CABI ABSTRACTS and

FAO AGRIS.

Keywords: Eco-epidemiology, Saturated incidence rate, Bifurcation, General Holling Type, Emergent

carrying capacity

essential approach for investigating and understanding the transmission

1 Introduction

and control of infectious diseases.

Numerous researchers (e.g.,

Hugo and Simanjilo (2019) and Sieber et al. (2014)) have explored

predator–prey models incorporating disease dynamics, highlighting how

infections within the prey and/or predator populations can signiﬁcantly

inﬂuence the ecological interactions and system stability. The primary

focus of eco-epidemiological models revolves around how infections

impact species mortality, decrease reproduction rates, the nature of

contamination, changes in population size, the eradication or control

of epidemic outbreaks, the persistence and the overarching dynamics

of the diseased species (Sieber et al., 2014). Saifuddin et al. (2016)

demonstrated that, under an explicit carrying capacity, susceptible and

infected prey exhibit identical competitive abilities, whereas under an

emergent carrying capacity, infected prey compete less eﬀectively than

susceptible ones in the presence of disease. Biswas et al. (2015) examined

a modiﬁed Lotka-Volterra system that incorporates the prey infection

propagation term based on the mass action law, while Haldar et al.

(2021) focused on standard incidence within predator-prey interactions.

Liu et al. (1987) proposed an epidemiological model characterized by

a nonlinear incidence rate. Gumel and Moghadas (2003) formulated

In applied mathematics, mathematical modeling serves as an essential

tool for investigating real-world problems across diverse disciplines,

including biology, epidemiology, and ecology (Bezabih et al., 2021).

Numerous researchers have demonstrated that the dynamic interactions

between predator and prey populations can be eﬀectively analyzed

using the tools of mathematical ecology (Das, 2016; Demir, 2019).

Building upon the foundational works of Lotka (1925) and Volterra

(1927), various sophisticated predator–prey models have been developed

to describe complex ecological interactions under diﬀerent realistic

scenarios (Ghanbari, 2021; Sieber et al., 2014). Furthermore, Anderson

and May (1986) established a pioneering framework that integrates the

epidemiological models of Kermack and McKendrick (Brauer, 2005)

with classical Lotka–Volterra predator–prey dynamics. As a result,

recent decades have been marked by a growing body of research

devoted to analyzing the dynamical behavior of eco-epidemiological

models (Biswas et al., 2015).

Since conducting experiments is

often impractical or unethical, mathematical modeling has become an

Alemu et.al (2026)

1

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

a tritrophic dynamics that incorporating a distinct saturating incidence

term to more accurately capture complex transmission dynamics, while

Ruan and Wang (2003) extended this line of research by examining

an epidemic model that integrates essential system with a saturating

To address these gaps, the study proposes a novel eco-epidemiological

model integrating saturated incidence, generalized Holling-type

predation on susceptible and infected prey, and an emergent carrying

capacity framework with distinct competition eﬀects. This integrated

approach strengthens theoretical understanding of system stability,

persistence, and complex population oscillations.

incidence term to investigate the overall system behavior.

Their

approach is deemed more justiﬁable because it considers behavioral

changes and the crowding eﬀect among infected individuals, thus

preventing the contact rate from becoming unbounded by selecting

appropriate parameters (Maiti et al., 2019). Hu et al. (2017) analyzed

a discrete-time eco-epidemiological framework, focusing on the system

dynamic behavior under a Holling type-II incidence function in place

of the bilinear incidence rate. Following these inﬂuential studies have

incorporated disease transmission into prey and/or predator populations

under various incidence mechanisms, including mass action, standard

incidence, and nonlinear forms. Among these, saturated incidence

rates have attracted considerable attention because they incorporate

behavioral changes and crowding eﬀects, thereby preventing unrealistic

unbounded transmission when the infected population becomes large.

Such formulations provide a more biologically realistic representation of

disease spread.

The remainder of this paper is organized as follows: Section 2 presents

the mathematical formulation of the model; Section 3 establishes

fundamental dynamical properties; Section 4 is devoted to stability

and bifurcation analysis; Section 5 provides numerical simulations

that support the analytical ﬁndings; Finally, the concluding section

summarizes the main results and discusses their ecological implications.

2 Mathematical Model

In this section, we investigate the eco-epidemiological dynamics to

explore the inﬂuence of a saturated incidence function on the sustainable

coexistence of two interacting species within the same ecosystem. Let ꢀ(C)

and ꢁ(C) denote the prey and predator densities at time C, respectively.

The model is formulated based on the following biological assumptions:

From an ecological perspective, predator–prey dynamics are strongly

inﬂuenced by the prey’s response to predation, while the predator

population, in turn, directly or indirectly regulates the prey population

(Panja, 2020). In order to accurately characterize the responsiveness

of predation rates to variations in prey biomass across diﬀerent

population densities, ecologically realistic functional responses have been

formulated that explicitly incorporate prey behavioral patterns. The

following functional responses are developed: Beddington-DeAngelis (Li

& Takeuchi, 2011), Crowley-Martin (Maiti et al., 2019), General Holling

type (Dawed et al., 2020), Michaelis-Menten type (HT-II), Holling type

III, Holling type IV (which came later) (Holling, 1959). Holling responses

are commonly categorized into speciﬁc forms (Type I–IV), each with

distinct ecological characteristics. However, in this study, the use of

the term “General Holling-Type functional responses” is intentionally

and methodologically justiﬁed. We mean either of these forms or

combinations of them:

The total prey population is divided into two compartments

ꢀ(C) = ꢀ (C) + ꢀ (C)

1

2

where ꢀ (C) and ꢀ (C) represent the susceptible and infected prey

1

2

populations, respectively.

The researchers assume that the lifespan of infected prey is shorter than

that of susceptible prey (Haldar et al., 2021). The susceptible prey

population ꢀ (C) follows logistic growth in the absence of predation

and disease. ¹Furthermore, both susceptible and infected prey share

limited environmental resources. However, they do not possess identical

competitive abilities. To capture this ecological feature, we incorporate

distinct competition coeﬃcients representing emergent carrying capacity:

0G

1 + G

0G²

1 + G²

0G

5(G) = 0G, ,(G) =

,

ℎ(G) =

,

A(G) =

,

1

denotes intra-speciﬁc competition among susceptible prey, while 1

2

1

1 + 1G + 2G²

represents inter-speciﬁc competition between susceptible and infected

prey (Ghanbari, 2021; Sieber et al., 2014). Thus, the logistic growth of

susceptible prey is given by

where, 0 is attack rate, 1 is a half saturation constant and 2 is the

measure of the predator tolerance to the prey to attack. Haque and

Venturino (2007) studied an eco-epidemic model in which the predator

population is infected and predation follows a ratio-dependent functional

response. Moreover, Kooi et al. (2011) also have discussed on tritrophics

food web eco-epidemiological system with predator infection, where the

infection transmitted among predators follow a hybrid response function

as Holling type-IV functional response and Beddington–DeAngelis type

functional response (Li & Takeuchi, 2011). Capasso and Serio (1978)

introduced an interaction term to account for the saturation eﬀect in

large infectious populations. Consequently, incorporating saturation

in disease transmission (Cai & Li, 2010) becomes particularly relevant

in eco-epidemiological models when the number of infectives is high.

Real-world predation involves complex mechanisms (Wayesa et al., 2024,

2025) such as prey refuge, handling time, predator interference, and

adaptive feeding, which can be captured using general Holling-type

functional responses. However, most eco-epidemiological models rely on

simpliﬁed predation terms and standard disease transmission functions,

with limited attention given to combining generalized predation

dynamics and saturated incidence. Key research gaps include:

3ꢀ

3C

1

= A_ꢀꢀ₁(1 − 1 ꢀ₁− 1 ꢀ ) .

1

2

1

The disease spreads among prey solely through direct contact. Infected

prey do not recover or acquire immunity; instead, they are removed from

the system through predation, disease-induced mortality at rate ꢀ, and

natural death at rate ꢁ .

1

We assume that susceptible prey become infected according to a

nonlinear saturated incidence function

ꢂꢀ ꢀ

2

1 +¹Bꢀ

2

as proposed in Maiti et al. (2019). Here, ꢂꢀ represents the force

2

1

of infection rate, while

accounts for behavioral changes and

1 + Bꢀ

2

crowding eﬀects among infected individuals. This formulation prevents

the transmission rate from becoming unbounded for large infected

populations (Ruan & Wang, 2003).

• Lack of systematic analysis of the combined eﬀects of saturated

disease transmission and general Holling-type responses on system

stability.

Ecologically, infected prey is generally more vulnerable to predation due

to its weakened physiological condition. To capture this phenomenon, we

incorporate distinct general Holling-Type functional responses, Φ_ꢀ(ꢀ )

and Φ_ꢀ(ꢀ ), which represent the predator’s consumption of susce¹ptib¹le

• Limited exploration of how these nonlinear mechanisms drive

qualitative changes such as transcritical and Hopf bifurcations.

2

• Few models incorporating emergent carrying capacity with

unequal competition between susceptible and infected prey.

and infected prey biomass, respectively.

Accordingly, the susceptible prey dynamics are given by

• The absence of a rigorous analytical framework linking the

basic reproduction number, stability switching, and bifurcation

dynamics Wang et al. (2016) under such generalized conditions.

ꢂꢀ ꢀ

3ꢀ

3C

1

2

1

= A_ꢀꢀ₁(1 − 1 ꢀ₁− 1 ꢀ ) −

− Φ_ꢀ(ꢀ )ꢁ,

1

2

1

1 + Bꢀ

2

Alemu et.al (2026)

2

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

while the infected prey dynamics are described by

3.1 Positivity of the solution

ꢂꢀ ꢀ

3ꢀ

3C

2

1 +¹Bꢀ

2

Let us denote R³= {((, ꢃ, %) ∈ R³: ( > 0, ꢃ ≥ 0, % > 0}, the positive

=

− (ꢁ₁+ ꢀ)ꢀ₂− Φ (ꢀ )ꢁ.

2

ꢀ

2

octants of the so⁺lution of our model system (3).

2

The model assumes a specialist predator population ꢁ(C) that feeds

on both susceptible and infected prey, with predation governed by

general Holling-type functional responses. Accordingly, the predator’s

population dynamics are formulated based on these generalized

predation interactions.

T heorem 1. The non-negative octant in R³is remain positive under the

dynamics for the model (3).

Proof. We want to verify

3ꢁ

3C

= ꢂ Φ (ꢀ )ꢁ + ꢂ Φ (ꢀ )ꢁ − ꢁ ꢁ,

1

2

ꢀ

2

1

(()) > 0,

ꢃ()) > 0,

%()) > 0,

for all , ) ≥ 0,

where ꢂ and ꢂ denote the conversion eﬃciencies of susceptible and

1

2

Rewrite the system (3) in the form

infected prey into predator biomass, respectively, and ꢁ represents the

2

natural mortality rate of the predator.

ꢀ

ꢁ

ꢄꢃ

1 + ꢅꢃ

2

ꢆ # (()%

3(

3)

1

(

= ( 1 − ( − ꢃꢃ −

−

= (& ((, ꢃ, %),

1

The descriptions of state variables and parameters are provided in Table 1.

All parameters are assumed to be positive. Hence, based on the above

assumptions, the governing eco-epidemiological model takes the form

(

ꢂ

ꢃ

ꢄ(

ꢆ # (ꢃ)%

3ꢃ

3)

ꢃ

ꢅ

= ꢃ

− ꢇ −

= ꢃ& ((, ꢃ, %),

2

1 + ꢅꢃ

ꢀ

ꢁ

ꢄ



ꢂꢀ ꢀ

3%

3)

3ꢀ

3C

1

2

1



= % #₍(() + #_ꢃ(ꢃ) − ꢈ = %& ((, ꢃ, %).

= A_ꢀ

ꢀ

1 − 1 ꢀ₁− 1 ꢀ

−

− Φ_ꢀ(ꢀ )ꢁ,

3

1

2

1

1 + Bꢀ

2





From the above expression and the initial conditions (4), we have:

ꢂꢀ ꢀ

3ꢀ

3C

1 +¹Bꢀ

2

(1)

(2)

=

− (ꢁ₁+ ꢀ)ꢀ₂− Φ (ꢀ )ꢁ,

2

ꢀ

2



2



3ꢁ

3C



= ꢂ Φ (ꢀ )ꢁ + ꢂ Φ (ꢀ )ꢁ − ꢁ ꢁ,

1

2

ꢀ

2

1



ꢂ

ꢃ

¹

)

(()) = (₀exp

& ((, ꢃ, %)3D

,

1

with initial conditions

0

ꢂ

ꢃ

¹

ꢀ (0) = ꢀ⁰> 0, ꢀ (0) = ꢀ⁰≥ 0, ꢁ(0) = ꢁ⁰> 0.

)

1

2

1

2

ꢃ()) = ꢃ₀exp

& ((, ꢃ, %)3D

2

,

0

ꢂ

¹

)

%()) = % exp

& ((, ꢃ, %)3D

3

.

2.1 Non-Dimensionalization

0

Non-Dimensionalization simplify and make the equations easier to

interpret. The transformation equations could be:

As, the initial conditions (4) and the exponential form are positive, thus,

1

B

1

all the state variables (()), ꢃ()) and %()) are positive ∀) ≥ 0. Therefore,

ꢀ

=

(, ꢀ

A_ꢀ

=

₁ꢃ, ꢁ = ₁%, C =

), Φ_ꢀ

1

2

1

A_ꢀ

1

every solutions of the mathematical model 3 are positive.

1

A_ꢀ

(ꢀ ) = ¹#₍((), and Φ_ꢀ(ꢀ ) = ¹#_ꢃ(ꢃ).

1

2

ꢂ

2

1

3.2 Bounded behavior of trajectories

Thus, the scaled form of the dynamical system is

ꢀ

ꢁ

ꢄ(ꢃ

1 + ꢅꢃ

 3(



T heorem 2. All possible solution of the dynamical system (3) are consistently

= ( 1 − ( − ꢃꢃ

−

− ꢆ # (()%,

1

(



3)

bounded in R³and enter in the invariant zone

+





ꢄ(ꢃ

3ꢃ

(3)

=

− ꢇꢃ − ꢆ # (ꢃ)%,

2

ꢃ



3)

1 + ꢅꢃ

= #₍(()% + #_ꢃ(ꢃ)% − ꢈ%,



nꢀ

ꢁ



Σ = (()), ꢃ()), %()) ∈ R³: 0 < ( ≤ max{( , 1},

0

3%

3)

+



ꢆ

ꢇꢇ



2

(1 + ꢈ − <)

0 < ꢉ ≤ max

, ꢉ

(5)

0

ꢂ

1

B

1

ꢁ₁+ ꢀ

4(ꢈ − <)

where ꢃ = ², ꢆ

=

, ꢆ

2

=

, ꢄ =

, ꢅ =

, ꢇ =

,

1

ꢂ

1 A

1

A_ꢀ

1

2

ꢀ

1

ꢁ

2

where ꢉ()) = (()) + ꢃ()) + ꢆ %()), 0 < # (ꢃ) ≤ <.

1

ꢃ

ꢈ =

, and

A_ꢀ

1

((0) = ( > 0, ꢃ(0) = ꢃ ≥ 0, %(0) = % > 0.

(4)

0

Proof. As established in Proposition (1), the solutions (()), ꢃ()), and %())

of system (3) remain positive for all ) ≥ 0. Considering the ﬁrst equation

of the model (3), it follows that

3 Mathematical Model Analysis

ꢀ

ꢁ

ꢄ(ꢃ

1 + ꢅꢃ

3(

3)

= ( 1 − ( − ꢃꢃ

−

− ꢆ # (()% ≤ ((1 − ()

1

(

The analysis of mathematical models in eco-epidemiology provides

valuable insights into disease transmission dynamics, host–pathogen

interactions, and the ecological feedback mechanisms within the

system. Such analysis helps to explore system behavior, identify critical

parameters, and examine aspects like stability, bifurcation, and possible

outcomes, including disease outbreaks. Furthermore, properties such as

the existence and uniqueness of solutions, positivity, boundedness, as

well as permanence, persistence, and numerical simulations of the model

system 3, will be studied.

This directly leads to

ꢈ

ꢉ₋₁

ꢀ

ꢁ

1

(

0

(()) ≤ 1 +

− 1 4⁻⁾

=

.

(

₀+ (1 − ( )4⁻⁾

(

0

Therefore,

lim sup (()) ≤ max{( , 1}.

0

)→∞

Alemu et.al (2026)

3

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Table 1: The state variables and parameters description

Ecological Meaning

Variables/Parameters

Dimension

ꢀ

ꢁ

Susceptible prey density

Per Area

Per time

Area

1

2

Infected prey density

Predator density

Response functions

Φ_ꢀ(resp.Φ_ꢀ

A_ꢀ

)

1

2

Natural propagation rate of susceptible prey

Intra-speciﬁc competition coeﬃcient among susceptible prey

Inter-speciﬁc competition coeﬃcient between susceptible and infected prey

Transmission rate

1

2

ꢂ

B

Area

Per time

Area coverage

No unit

Per time

Inhibition rate

ꢂ

Proportion of susceptible prey into predator

Proportion of infected prey into predator

Natural death rates of infected prey/predator

Disease induced mortality rate

1

2

ꢁ /ꢁ

1

2

ꢀ

Thus, ((C) is bounded. To show other state variables ꢃ()) and %()) are

Proof. Let the right parts of the dynamical system (3) be denoted by ꢄ =

bounded we consider

( 5 , 5 , 5 ). Since 5 , 5 and 5 are continuous function, ꢄ = ( 5 , 5 , 5 )

2

3

1

is¹continuous function²in seve³ral variables, that is, ꢄ ∈ ꢂ¹(R³).¹Th²us,³ꢄ

satisfy the Lipschitz condition with respect to G in ꢅ. Hence, t⁺he solution

of system (3) exists. The locally Lipschitz condition of ꢄ is veriﬁed using

ꢉ = ( + ꢃ + ꢆ %

1

By diﬀerentiating ꢉ with respect to time ), we obtain

% 5₈

, 8, 9 = 1, 2, 3 to be continuously bounded within the domain ꢅ

3ꢉ

3)

3(

3)

3ꢃ

3)

3%

%G₉

=

+

+ ꢆ

¹3)

3(

3)

(Bezabih et al., 2021). We note that ꢄ ∈ ꢂ¹(R³, !8?) in ꢅ and 5

=

,

ꢀ

ꢁ

ꢄ(ꢃ

1 + ꢅꢃ

1

+

= ( 1 − ( − ꢃꢃ

−

− ꢆ # (()%

1

(

% 5₈

3ꢃ

3%

3)

5

=

and 5

=

. To show

, 8, 9 = 1, 2, 3 to be continuously

ꢂ

ꢃ

2

3

3)

%G₉

ꢄ (ꢃ

1 + ꢅꢃ

2

+

− ꢇꢃ − ꢆ # (ꢃ)%

2

ꢃ

bounded. Now we get

ꢄ

ꢅ

+ ꢆ #₍(()% + #_ꢃ(ꢃ)% − ꢈ%

1

ꢀ

ꢁ

= ( 1 − ( − ꢃ(ꢃ − ꢇꢃ − ꢆ # (ꢃ)% + ꢆ # (ꢃ)% − ꢆ ꢈ%

2

1

ꢃ

ꢄ

ꢅ

% 5

ꢄꢃ

1 + ꢅꢃ

1

= 1 − 2( − ꢃꢃ −

− ꢆ #⁰(()% ≤ 1,

ꢀ

ꢁ

1

(

%(

% 5

≤ ( 1 − ( − ꢃ(ꢃ − ꢇꢃ + ꢆ # (ꢃ)% − ꢆ ꢈ%

1

ꢃ

ꢂ

ꢃ

ꢄ(

1

= (1 + ꢈ)( − (²− ꢃ(ꢃ − (ꢇ − ꢈ)ꢃ + ꢆ # (ꢃ)% − ꢈ(( + ꢃ + ꢆ %)

= −ꢃ( −

= −( ꢃ +

,

2

1

ꢃ

%ꢃ

(1 + ꢅꢃ)

ꢊ

ꢂ

ꢃ

≤ (1 + ꢈ)( − (²+ ꢆ # (ꢃ)% − ꢈꢉ

1

3 5

1

3ꢃ

ꢄ(

ꢃ

This implies

= ( ꢃ +

< ∞, as ( and ꢃ are bounded,

2

(1 + ꢅꢃ)

ꢊ

3 5

1

3%

3 5

1

ꢊ

The general Holling Type response function #_ꢃ(ꢃ) is bounded, say, #_ꢃ(ꢃ) ≤

= −ꢆ # (() implies

= −ꢆ # (() = ꢆ # (() < ꢆ # ,

1

(

3%

<, 0 < < < ꢈ ≤ ꢇ. Then after simpliﬁcation we arrive

as # (() ≤ # ∈ R,

1

(

2

3ꢉ

3)

(1 + ꢈ − <)

% 5

2

ꢄꢃ

ꢄ

ꢅ

≤ (1 + ꢈ − <)( − (²− (ꢈ − <)ꢉ ≤

− (ꢈ − <)ꢉ.

=

≤

,

4

%(

% 5

1 + ꢅꢃ

ꢄ(

2

We can thus conclude that

− ꢇ − ꢆ #⁰(ꢃ)% ≤

< ∞,

2

ꢃ

2

%ꢃ

(1 + ꢅꢃ)

ꢂ

ꢃ

2

(1 + ꢈ − <)

ꢉ()) ≤

(1 + ꢈ − <)

− ꢉ 4^{−(ꢈ−<))}.

as ( and ꢃ are bounded,

−

0

4(ꢈ − <)

ꢊ

3 5

2

3%

3 5

2

3%

ꢊ

= −ꢆ # (ꢃ) =⇒

= −ꢆ # (ꢃ) = ꢆ # (ꢃ) < ꢆ # ,

2

ꢃ

As a result, we ﬁnd that

as # (ꢃ) ≤ # ∈ R,

= #₍⁰(()% < ∞,

2

ꢃ

ꢆ

ꢇ

2

(1 + ꢈ − <)

% 5

3

lim sup ꢉ()) ≤ max

, ꢉ

.

0

4(ꢈ − <)

)→∞

%(

% 5

3

= #⁰_ꢃ(ꢃ)% < ∞,

Hence, ꢉ()) remains bounded for all ) ≥ 0, which implies that the other

%ꢃ

% 5

state variables are also bounded. Consequently, all solutions of system (3)

3

are uniformly bounded on [0, ∞).

= # (() + # (ꢃ) − ꢈ < # (() + # (ꢃ) ≤ #₁+ # < ∞.

2

(

ꢃ

(

ꢃ

%%

As these all are continuous and bounded, ꢄ satisfy the locally Lipschitz

condition. Therefore, the unique solution of the system (3) is veriﬁed as

3.3 Existence & Uniqueness

it is explained in Allen et al. (2007) and Hale (2009).

T heorem 3. Let ꢄ = ( 5 , 5 , 5 ). If ꢄ satisﬁes the Lipschitz condition and has

continuous ﬁrst partial deriv²atives with respect to G in a domain ꢅ, then ꢄ(), G)

1

3

3.4 Equilibrium points

is locally Lipschitz in G. Consequently, for any initial point () , G ) ∈ ꢅ, there

0

exists a unique solution G(), ) , G ) of the system

0

The ﬁxed points of the dynamical system (3) are the roots of a nonlinear

system of equations.

3G

3)

= ꢄ(), G), G() ) = G ,

0

ꢀ

ꢁ

ꢄ(ꢃ

1 + ꢅꢃ

( 1 − ( − ꢃꢃ

−

− ꢆ # (()% = 0,

(6)

1

(

which passes through () , G ).

0

Alemu et.al (2026)

4

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

ꢄ(ꢃ

1 + ꢅꢃ

2ꢅꢇ

ꢄ

− ꢇꢃ − ꢆ # (ꢃ)% = 0,

(7)

(8)

• If ꢇ > ꢄ (i.e., Δ > 0) and

+ ꢄ + ꢃ > ꢅ (i.e Δ > 0), then (11) has

2

3

2

ꢃ

no positive root meaning that there is no feasible equilibrium point

#₍(()% + #_ꢃ(ꢃ)% − ꢈ% = 0.

ꢆ².

Hence, the extinction ﬁxed point is ꢆ⁰(0, 0, 0), the axial ﬁxed point is

ꢆ¹(1, 0, 0).

• If ꢇ < ꢄ (i.e., Δ < 0), then there exists a unique equilibrium point

3

The predator free equilibrium point ꢆ²is obtained by the intersection

ꢆ².

ꢂ

ꢃ

3(

3)

point of the zero growth isocline of susceptible

= 0 and the zero

2ꢅꢇ

ꢄ

ꢂ

ꢃ

• If ꢇ > ꢄ (i.e., Δ > 0) and

+ ꢄ + ꢃ < ꢅ (i.e., Δ < 0), then

3ꢃ

3)

3

2

growth isocline of infected species

= 0 where % = 0. That is,

equation (11) has two positive roots, consequently two equilibrium

points ꢆ₁², and ꢆ₂².

ꢀ

ꢁ

ꢄ(^∗ꢃ^∗

1 + ꢅꢃ^∗

∗

(

1 − (^∗− ꢃꢃ^∗

−

ꢄ(^∗ꢃ^∗

1 + ꢅꢃ^∗

= 0,

(9)

− ꢇꢃ^∗= 0,

(10)

Disease free equilibrium point

The infection free ﬁxed point of the form ꢆ³((^∗, 0, %^∗) is solution of

From equation (10), we get

non-linear system

ꢇꢅ

ꢃ^∗

ꢇ

(^∗=

+

ꢀ

ꢁ

ꢄ

∗

(

1 − (^∗− ꢆ # ((^∗)%^∗= 0 and # ((^∗)%^∗− ꢈ%^∗= 0.

1

(

Substitute this equation in (9), after simpliﬁcation we arrived

ꢀ

ꢁ

Δ ꢃ^∗2+ Δ ꢃ^∗+ Δ = 0,

(11)

∗

(

1 − (^∗

1

2

3

This gives #₍((^∗) = ꢈ and %^∗=

.

ꢂ

ꢃ

ꢇꢅ

ꢄ

2ꢅꢇ

ꢄ

ꢇ

ꢈꢆ

where Δ = ꢅ ꢃ +

> 0, Δ

=

+ ꢄ + ꢃ − ꢅ and Δ

=

− 1.

1

2

3

ꢄ

The positive roots ꢃ^∗in the quadratic equation above is possible provided

Table 2 provides an explanation of the illness free equilibrium point’s

existence criteria. where ℎ is half saturation constant and F denote

predator attack rate.

that the discriminant of an equation is positive, that is, Δ²₂− 4Δ Δ > 0

and follow from Descartes’ rule of sign. We have the following r¹esults:

2

Table 2: Existence conditions of the disease free ﬁxed point.

HT

HT-I

HT-II

HT-III

ℎ²ꢈ

HT-IV

p

q

2

4($−ꢈ)+ 4²(ꢈ−$) −4ꢈℎ4

ꢈ

ℎꢈ

∗

(

$−ꢈ

2ꢈ

√

$

$ − ꢈ

Conditions

ꢈ < $

ℎꢈ < $ − ꢈ and $ > ꢈ

ℎ²ꢈ < $ − ꢈ and $ > ꢈ

(ꢈ − $)4 > 2 ꢈℎ

The basic reproduction number, '

Thus,

0

V (ꢃ) ≈ ꢇꢃ + ꢆ #⁰(0)%^∗ꢃ.

2

ꢃ

According to Layek (2015), the basic reproduction number is the average

number of new infections from a single sick individual in a community

that is completely susceptible over the course of the infectious period. It is

used to predict whether the epidemic will spread or die out (Omar et al.,

2024). To compute the basic reproduction number, we consider only the

infected compartment of system (3)

The linearized equation becomes

ꢋ

ꢌ

ꢄ

ꢅ

3ꢃ

3)

=

ꢄ(^∗− ꢇ + ꢆ #⁰(0)%^∗ꢃ.

2

ꢃ

ꢄ(ꢃ

3ꢃ

=

− ꢇꢃ − ꢆ # (ꢃ)%.

(12)

2

ꢃ

3)

1 + ꢅꢃ

Hence, the new infection rate is

Following the next-generation matrix approach, we write

ꢄ = ꢄ(^∗,

3ꢃ

3)

= ℱ (ꢃ) − V (ꢃ),

and the total removal rate is

where the new infection and the transition (removal) terms represent

+ = ꢇ + ꢆ #⁰(0)%^∗.

2

ꢃ

ꢄ(ꢃ

1 + ꢅꢃ

ℱ (ꢃ) =

, and V (ꢃ) = ꢇꢃ + ꢆ # (ꢃ)%.

2

ꢃ

By the next-generation method,

Since ' measures the invasion of infection when ꢃ is small, we linearize

the syst⁰em around ꢃ = 0. Using Taylor expansion,

'

= ꢄ+⁻¹

.

0

ꢄ(ꢃ

1 + ꢅꢃ

= ꢄ(ꢃ(1 − ꢅꢃ + $(ꢃ²)).

Therefore,

ꢄ(^∗₀

keeping only ﬁrst-order terms gives ℱ (ꢃ) ≈ ꢄ(ꢃ. Similarly, expanding

'

=

.

(13)

#_ꢃ(ꢃ) near ꢃ = 0, #_ꢃ(ꢃ) ≈ #⁰(0)ꢃ.

0

ꢇ + ꢆ # (0)%^∗

2

ꢃ

Alemu et.al (2026)

5

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Coexistence Equilibrium Point

4.1 Local stability analysis

T heorem 4. The system admits a coexistence equilibrium ꢆ^∗((^∗, ꢃ^∗, %^∗) if the

4.1.1 Stability nature near ꢆ⁰(0, 0, 0)

following conditions hold:

ꢄ(^∗

1 + ꢅꢃ^∗

(^∗+ ꢃꢃ^∗< 1,

> ꢇ, #₍((^∗) + #_ꢃ(ꢃ^∗) = ꢈ,

'

> 1

0

.

1

0

©

-

ª

®

ꢇ(ꢆ⁰) =

^-0 −ꢇ

-

Proof. To determine the coexistence equilibrium point ꢆ^∗((^∗, ꢃ^∗, %^∗), we

0

−ꢈ

«

¬

set

3(

3)

3ꢃ

3)

3%

3)

= 0,

= 0.

Thus, ꢊ = 1 > 0, ꢊ = −ꢇ < 0, and ꢊ = −ꢈ < 0. Hence, the trivial ﬁxed

point ꢆ⁰is unstable². Biologically, this³indicates that total extinction of the

1

Thus the equilibrium point satisﬁes the algebraic equations

populations is impossible.

ꢄ(^∗ꢃ^∗

(^∗(1 − (^∗− ꢃꢃ^∗) −

− ꢆ # ((^∗)%^∗= 0,

(3.9)

1

(

1 + ꢅꢃ^∗

4.1.2 System behavior near ꢆ¹(1, 0, 0)

ꢄ(^∗ꢃ^∗

− ꢇꢃ^∗− ꢆ # (ꢃ^∗)%^∗= 0,

(3.10)

(3.11)

2

ꢃ

1 + ꢅꢃ^∗

#₍((^∗)%^∗+ #_ꢃ(ꢃ^∗)%^∗− ꢈ%^∗= 0.

ꢄ(^∗

1 + ꢅꢃ^∗

−1 −ꢃ − ꢄ −ꢆ #₍(1)

ꢄ − ꢇ

¹0

©

-

ª

®

From the ﬁrst, second and third equations we have (^∗+ꢃꢃ^∗< 1,

>

ꢇ(ꢆ¹) =

0

ꢇ, and #₍((^∗) + #_ꢃ(ꢃ^∗) = ꢈ, respectively.

0

#₍(1) − ꢈ

«

¬

Note that

ꢈꢆ ꢄ #⁻¹(ꢈ)

1

(

'

=

ꢀ

ꢁ .

0

ꢇꢈꢆ₁+ ꢆ #⁰(0)#⁻¹(ꢈ) 1 − #₍⁻¹(ꢈ)

The eigenvalues are ꢊ = −1, ꢊ = ꢄ − ꢇ, and ꢊ = #₍(1) − ꢈ. Thus,

1

2

the axial ﬁxed point is locally asymptotically sta³ble whenever ꢄ < ꢇ

and #₍(1) < ꢈ. This has a biological implication that susceptible prey

population survive alone whenever no disease in the environment and

without predator whenever the conditions holds.

2

ꢃ

(

ꢄ(^∗

1 + ꢅꢃ^∗

ꢄ(^∗

ꢇ(1 + ꢅꢃ^∗)

> ꢇ ⇐⇒

> 1.

At equilibrium (using (^∗= #₍⁻¹(ꢈ) we obtain

4.1.3 System behavior near consumer free fixed point ꢆ²((^∗, ꢃ^∗, 0)

ꢄ(^∗

ꢇ(1 + ꢅꢃ^∗)

= ' .

0

T heorem 5. The consumer-free ﬁxed point ꢆ²((^∗, ꢃ^∗, 0) is locally

asymptotically stable if the following conditions hold

ꢄ(^∗

> ꢇ ⇐⇒

'

> 1.

0

1 + ꢅꢃ^∗

(i) #₍((^∗) + #_ꢃ(ꢃ^∗) < ꢈ,

(ii) ꢇ + 2(^∗+ ꢃꢃ^∗+

ꢄꢃ^∗

ꢄ(^∗

< 1 +

,

1 + ꢅꢃ^∗

(1 + ꢅꢃ^∗)²

ꢂ

ꢃ

ꢄꢃ^∗

4 Stability and Bifurcation

Analysis

(iii)

1 − 2(^∗− ꢃꢃ^∗−

1 + ꢅꢃ^∗

ꢂ

ꢃ

ꢄ(^∗

×

+

− ꢇ

(1 + ꢅꢃ^∗)²

ꢂ

ꢃ

ꢄ(^∗

ꢄꢃ^∗

1 + ꢅꢃ^∗

By examining sign of the derivative matrix’s eigenvalues, we can

determine the stability of a ﬁxed points as in Dawed et al. (2020). The

system (3) has a stable ﬁxed point ꢆ^∗((^∗, ꢃ^∗, %^∗) if all characteristic roots of

the Jacobian matrix, ꢇ(ꢆ^∗),

ꢃ(^∗+

> 0.

(1 + ꢅꢃ^∗)²

ꢇ₁₁ꢇ₁₂ꢇ₁₃

^-ꢇ₂₁ꢇ₂₂ꢇ₂₃

ꢇ₃₁ꢇ₃₂ꢇ₃₃

©

ª

®

-

ꢇ(ꢆ^∗) =

(14)

Proof. The community matrix of the model (3) at ꢆ²is given by

-

«

¬

have negative real part where

ꢄꢃ^∗

ꢄ(^∗

(1 + ꢅꢃ^∗)²

ꢄꢃ^∗

1 + ꢅꢃ^∗

ꢄ(^∗

ꢇ

= 1 − 2(^∗− ꢃꢃ^∗−

− ꢆ #⁰((^∗)%^∗, ꢇ = −ꢃ(^∗−

,

1 − 2(^∗− ꢃꢃ^∗−

−ꢃ(^∗−

−ꢆ # ((^∗)

11

1

12

1 + ꢅꢃ^∗

©

-

ª

®

(

1

(

(1 + ꢅꢃ^∗)

2

ꢄꢃ^∗

ꢄ(^∗

2

ꢄꢃ^∗

ꢄ(^∗

− ꢇ − ꢆ #⁰(ꢃ^∗)%^∗,

ꢇ(ꢆ ) =

ꢇ

= −ꢆ # ((^∗), ꢇ

=

, ꢇ

22

=

13

1

21

2

− ꢇ

−ꢆ # (ꢃ^∗)

(

2

ꢃ

1 + ꢅꢃ^∗

ꢃ

(1 + ꢅꢃ^∗)²

1 + ꢅꢃ^∗

2

(1 + ꢅꢃ^∗)

= −ꢆ # (ꢃ^∗), ꢇ = #⁰((^∗)%^∗, ꢇ = #⁰(ꢃ^∗)%^∗, and ꢇ = #₍((^∗) +

(

ꢃ

0

# ((^∗) + # (ꢃ^∗) − ꢈ

«

¬

23

2

31

32

33

ꢃ

(

ꢃ

#_ꢃ(ꢃ^∗) − ꢈ.

The associated auxiliary equation is 34C(ꢇ(ꢆ²) − ꢊꢃ ) = 0, that is

3

The corresponding characteristic equation is det(ꢇ(ꢆ^∗) − ꢊꢃ ) = 0, that is

3

ꢊ

∗

ꢋ − ꢊ

ꢌ

−ꢆ # (( )

1

(

ꢊ

^ꢊꢇ − ꢊ

ꢇ₁₂

ꢇ₂₂− ꢊ

ꢇ₃₂

ꢇ₁₃

ꢇ₂₃

∗

ꢊ

11

ꢇ₂₁

ꢇ₃₁

ꢌ

ꢋ − ꢊ −ꢆ # (ꢃ ) = 0,

2

ꢃ

ꢊ

= 0

(15)

ꢊ

0

ꢋ − ꢊ

3

ꢇ₃₃− ꢊ^ꢊ

ꢊ

Alemu et.al (2026)

6

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

ꢄꢃ^∗

ꢄ(^∗

(1 + ꢅꢃ^∗)²

where, ꢋ = 1 − 2(^∗− ꢃꢃ^∗−

, ꢋ

=

− ꢇ, ꢋ = #₍((^∗) +

we write

1

2

3

1 + ꢅꢃ^∗

ꢄ(^∗

ꢄꢃ^∗

ꢊ

= (ꢇ + ꢆ #⁰(0)%^∗)('₀− 1).

1

2

# (ꢃ^∗) − ꢈ, ꢌ = −ꢃ(^∗−

< 0 and ꢌ

=

> 0. Thus, one

ꢃ

1

2

ꢃ

1 + ꢅꢃ^∗

(1 + ꢅꢃ^∗)²

eigenvalue at ꢆ²((^∗, ꢃ^∗, 0) is ꢊ = ꢋ = # ((^∗) + # (ꢃ^∗) − ꢈ. Now, ꢊ is

1

3

1

(

ꢃ

Thus, if ' < 1, then ꢊ < 0 and the infection cannot invade. The others

two eigenvalues are computed from the matrix

0

1

negative if #₍((^∗) + #_ꢃ(ꢃ^∗) < ꢈ. The rest two eigenvalues are found from

the matrix

ꢄꢃ^∗

1 + ꢅꢃ^∗

ꢄ(^∗

(1 + ꢅꢃ^∗)²

1 − 2(^∗− ꢆ #⁰((^∗)%^∗−ꢆ # ((^∗)

1 − 2(^∗− ꢃꢃ^∗−

−ꢃ(^∗−

1

(

©

-

ª

®

©

-

ª

®

ꢇ =

2

#₍⁰((^∗)%^∗

0

¯

ꢇ(ꢆ ) =

.

«

¬

ꢄꢃ^∗

1 + ꢅꢃ^∗

ꢄ(^∗

− ꢇ

(1 + ꢅꢃ^∗)²

«

¬

By the Routh-Hurwitz stability rule, th₀e two eigenvalues o₀f ꢇ possess

negative real parts provided that ꢆ ꢈ# ((^∗)%^∗> 0 (i.e., #₍((^∗) > 0)

¯

Using the Routh–Hurwitz criterion, the two eigenvalues of ꢇ are negative

1

(

and #₍⁰((^∗)(^∗2− (2ꢈ + #₍⁰((^∗))(^∗+ ꢈ < 0. Thus, if conditions (1)–(3) are

satisﬁed, we conclude that the model system (3) is locally asymptotically

stable at the disease-free ﬁxed point ꢆ³. The predator eating eﬃciency is

so high whenever conditions (1)–(3) are satisﬁed. The predator will only

eat healthy prey because there is no infected prey present.

in their real parts provided that

ꢄꢃ^∗

1 + ꢅꢃ^∗

ꢄ(^∗

ꢇ + 2(^∗+ ꢃꢃ^∗+

< 1 +

and

(1 + ꢅꢃ^∗)²

ꢂ

ꢃ ꢂ

ꢃ

∗

ꢄꢃ

ꢄ(

∗

1 − 2( − ꢃꢃ −

− ꢇ +

∗

∗ 2

1 + ꢅꢃ

(1 + ꢅꢃ )

ꢂ

ꢃ

∗

ꢄ(

ꢄꢃ

∗

ꢃ( +

> 0.

∗

∗ 2

1 + ꢅꢃ

(1 + ꢅꢃ )

Therefore, we infer that the model system (3) is locally asymptotically

stable at the predator free equilibrium point ꢆ²as long as the conditions

4.1.5 Global stability analysis using the Bendixson-Dulac

theorem

(i), (ii), and (iii) hold.

T heorem 7. If the Infection Free Equilibrium point ((^∗, 0, %^∗) is locally

4.1.4 Local Stability Near the Disease-Free Equilibrium Point

asymptotically stable in the positive (% - plane region, then it is also globally

#₍(()

#₍((^∗)

T heorem 6. Local asymptotic stability of the infection-free ﬁxed point

asymptotically stable in the same region if

≥

.

ꢆ³((^∗, 0, %^∗) of system (3) is ensured if the following criteria are satisﬁed:

∗

(

1. ' < 1,

0

2. #₍⁰((^∗) > 0,

3. #₍⁰((^∗)((^∗)²− 2ꢈ + #₍⁰((^∗) (^∗+ ꢈ < 0.

Proof. Consequently, the system can be reduced to the following

two-dimensional subsystem

ꢄ

ꢅ

3(

= ((1 − () − ꢆ # (()%,

(17)

(18)

1

(

3)

3%

3)

Proof.

= #₍(()% − ꢈ%.

5₁

−ꢃ(^∗− ꢄ(^∗−ꢆ # ((^∗)

1

(

©

-

ª

®

1

(%

3

Consider ℎ((, %)

=

as a Dulac positive function in the positive

0

5₂

#⁰_ꢃ(0)%^∗

0

ꢇ(ꢆ ) =

(16)

-

®

quadrant. Also, deﬁne the following functions

#⁰((^∗)%^∗

5₃

«

¬

(

ℎ ((, %) = ((1 − () − ꢆ # (()%,

(19)

(20)

1

(

ℎ ((, %) = # (()% − ꢈ%.

2

where, 5₁= 1 − 2(^∗− ꢆ #⁰((^∗)%^∗, 5₂= ꢄ(^∗− ꢇ − ꢆ #⁰(0)%^∗and

(

1

2

(

ꢃ

5₃= #₍((^∗) − ꢈ = 0.

Then,

The associated auxiliary equation of ꢇ(ꢆ³) is

ꢂ

ꢃ

(#⁰(() − #₍(()

ꢊ

%

1

%

(

∗

ꢊ

)((, %) =

(ℎ ℎ ) +

(ℎ ℎ ) = −

+ ꢆ

.

5 − ꢊ

−ꢃ( − ꢄ(

−ꢆ # (( )

1

2

1

(

ꢊ

2

%(

%%

(

ꢊ

0

5 − ꢊ

2

0

= 0

Hence, )((, %) is a negative function of its arguments if (#₍⁰(() −

ꢊ

#₍(() ≥ 0. Note that by mean value thorem (#⁰(() − #₍(() ≥ 0 and

0

∗

# (( )%

(

∗

# (0)%

ꢃ

5 − ꢊ

3

(

#₍(()

#₍((^∗)

≥

are equivalent. Since )((, %) does not change sign and

∗

(

is not identically zero in the positive quadrant of the (%-plane, by the

Bendixson - Dulac criterion the infection free equilibrium point is globally

asymptotically stable in the region

Since the matrix is block triangular with respect to the infected variable,

one eigenvalue is

ꢊ

= 5₂= ꢄ(^∗− ꢇ − ꢆ #⁰(0)%^∗.

1

2

ꢃ

ꢆ

ꢇ

#₍(()

#₍((^∗)

ꢅ = ((, %) ∈ '²

:

≥

, % > 0

.

Using the basic reproduction number

+

∗

(

ꢄ(^∗₀

'

=

,

0

ꢇ + ꢆ # (0)%^∗

Moreover, the system has no limit cycle in the region.

7

2

ꢃ

Alemu et.al (2026)

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

4.1.6 System stability conditions near endemic equilibrium point

Proof. (1). Consider the community matrix of model (3) evaluated at

ꢆ¹(1, 0, 0):

−1 −(ꢃ + ꢄ) −ꢆ # (1)

T heorem 8. Local asymptotic stability of the endemic equilibrium point

©

-

ª

®

ꢇ(ꢆ¹) =

.

0

ꢄ − ꢇ

¹0⁽

ꢆ^∗((^∗, ꢃ^∗, %^∗) holds provided that the following conditions are met:

0

#₍(1) − ꢈ

«

¬

The eigenvalues of ꢇ(ꢆ¹) are ꢊ = −1 < 0, ꢊ = ꢄ − ꢇ, and ꢊ = # (1)− ꢈ.

(8) ꢅ + ꢈ < 0,

1

2

3

(

Therefore, ꢆ¹is locally asymptotically stable provided that ꢄ < ꢇ and

#₍(1) < ꢈ hold. Substituting either ꢇ = ꢄ or ꢈ = #₍(1) into ꢇ(ꢆ¹) yields a

zero eigenvalue in the characteristic equation.

(88) @ @₄− ꢆ@₅+ ꢈ@₂+ ꢅ@ < 0,

1

3

(888) (ꢅ + ꢈ)(ꢅꢈ + @₂+ @₃+ ꢆ@ ) + @ @₄− ꢆ@₅+ ꢈ@₂+ ꢅ@ > 0,

4

1

3

where the parameters ꢅ, ꢆ, ꢄ, and ꢈ are deﬁned in the proof.

With ꢇ

=

ꢇ^[1], the eigenvectors + and ,, associated with the zero

eigenvalue of the Jacobian ꢇ^[1](ꢆ¹, ꢇ^[1]) and its transpose, respectively,

are

Proof. The positive ﬁxed point ꢆ^∗((^∗, ꢃ^∗, %^∗) of the dynamics (3) is locally

asymptotically stable if all the characteristic roots of the Jacobian matrix,

ꢇ, has negative real parts, where

−(ꢃ + ꢄ)0

0

©

-

ª

®

© ª

0

1

+ =

,

, =

,

- ®

0

«

¬

« ¬

)

where 0, 1 ≠ 0. The derivative of the vector ﬁeld ꢄ((, ꢃ, %) = (ꢄ , ꢄ , ꢄ )

1

2

3

ꢅ

ꢆ

−ꢃ(^∗− ꢄ(^∗ꢄ²−ꢆ #₍((^∗)

1

with respect to ꢇ is

©

-

ª

®

ꢇ(ꢆ^∗) =

.

ꢈ

−ꢆ #_ꢃ(ꢃ^∗)

2

0

©

ª

®

© ª

ꢄ (-, ꢇ) = −ꢃ

,

ꢄ_ꢇ(ꢆ¹, ꢇ^[1]) = 0 ,

#⁰((^∗)%^∗

#⁰_ꢃ(ꢃ^∗)%^∗

0

-

- ®

ꢇ

«

¬

(

0

«

¬

« ¬

implying

1

,⁾ꢄ_ꢇ(ꢆ¹, ꢇ^[1]) = 0.

where, ꢅ = 1 − 2(^∗− ꢃꢃ^∗− ꢄꢃ^∗ꢄ − ꢆ #⁰((^∗)%^∗, ꢆ = ꢄꢃ^∗ꢄ, ꢄ =

and

1

1 + ꢅꢃ^∗

(

ꢈ = ꢄ(^∗ꢄ²− ꢇ − ꢆ #⁰(ꢃ^∗)%^∗.

Hence, the ﬁrst condition of Sotomayor’s theorem Pirayesh et al., 2016 for

a transcritical bifurcation is met.

2

ꢃ

The characteristic equation is

Next, we compute

ꢊ

ꢅ − ꢊ

−ꢃ(^∗− ꢄ(^∗ꢄ²^{−ꢆ #}₍⁽⁽^∗⁾ꢊ

0

−1

0

−0

0

1

ꢊ

©

-

ª

®

©

-

ª

®

ꢅꢄ_ꢇ(ꢆ¹, ꢇ^[1]) =

,

ꢅꢄ_ꢇ(ꢆ¹, ꢇ^[1])+ =

,

∗

ꢊ

ꢆ

ꢈ − ꢊ

#⁰_ꢃ(ꢃ^∗)%^∗

−ꢆ #_ꢃ(ꢃ )

2

= 0

ꢊ

«

¬

«

¬

0

#₍((^∗)%^∗

−ꢊ

so that

,⁾[ꢅꢄ_ꢇ(ꢆ¹, ꢇ^[1])+] = −01 ≠ 0.

This implies

Finally, the second derivative of ꢄ along + is

ꢊ³+ : ꢊ²+ : ꢊ + : = 0

(21)

1

2

3

−2E₁²+ 2ꢄꢅE₂²+ 2(−ꢃ − ꢄ)E E

1

2

©

-

ª

®

ꢅ²ꢄ(ꢆ¹, ꢇ^[1])(+, +) =

−2ꢄꢅE₂²+ 2ꢄE E

,

1

2

0

where, : = −(ꢅ + ꢈ), : = ꢅꢈ + @ + @ + ꢆ@ , : = −(@ @

−

1

2

3

4

3

4

ꢆ@ + ꢈ@ + ꢅ@ ), @ = ꢆ # (ꢃ^∗)#⁰((^∗)%^∗, @ = ꢆ # ((^∗)#⁰(¹(^∗)%^∗,

«

¬

5

2

3

1

2

1

ꢃ

(

giving

@

= ꢆ # (ꢃ^∗)#⁰(ꢃ^∗)%^∗, @ = ꢃ(^∗+ ꢄ(^∗ꢄ²and @ = ꢆ # ((^∗)#⁰(ꢃ^∗⁽)%^∗.

3

2

4

5

1

ꢃ

(

ꢃ

,⁾[ꢅ²ꢄ(ꢆ¹, ꢇ^[1])(+, +)] = −21ꢇ(ꢅE₂²− E E ) ≠ 0 whenever ꢅE ≠ E .

Consequentially, if ꢅ +ꢈ < 0, then : > 0. If @ @ −ꢆ@ +ꢈ@ +ꢅ@ < 0,

then : > 0. Moreover, : : − : > 0 if (ꢅ +¹ꢈ)(ꢅꢈ + @ +²@ + ꢆ@ ) +

1

4

5

3

1

2

1

3

1

2

3

2

3

@ @₄− ꢆ@₅+ ꢈ@₂+ ꢅ@ > 0. According to the Routh–Hurwitz criter⁴ion,

the model system 3 is l³ocally asymptotically stable at the endemic ﬁxed

1

Therefore, by Sotomayor’s theorem Pirayesh et al., 2016; Yu et al., 2020,

point ꢆ^∗= ((^∗, ꢃ^∗, %^∗) if the corresponding conditions are satisﬁed.

the model (3) demonstrate transcritical bifurcation at ꢇ = ꢇ^[1]near the

axial ﬁxed point ꢆ¹(1, 0, 0) provided that ꢅE ≠ E .

2

1

4.2 Local bifurcation analysis

Now, let us examine the bifurcation at ꢈ = ꢈ^[2]= #₍(1). The eigenvectors

+ and ,, associated to the zero eigenvalues of the matrices ꢇ^[2](ꢆ¹, ꢈ^[2]

)

and its transpose respectively, can be written as

Bifurcations analysis helps to predict and understand transitions in

the behavior of dynamical system as parameters value change. Local

bifurcation refers, change in the qualitative behavior of dynamical system

near ﬁxed point as a system’s parameters are varied.

ꢄ

ꢅ

)

+ = (E

E

E ) = −ꢆ # (1)

0

2

and , = (0

0

3)

1

2

3

1

(

where 2 and 3 are nonzero real numbers. It follows that

T heorem 9. (Transcritical Bifurcation)

0

−%

0

©

-

ª

®

© ª

0

ꢄ (-, ꢈ) =

,

ꢄ_ꢇ(ꢆ¹, ꢈ^[2]) =

.

- ®

ꢈ

0

1. The diseased model (3) demonstrate a transcritical bifurcation at parameter

«

¬

« ¬

values ꢇ = ꢇ^[1]= ꢄ or ꢈ = ꢈ^[2]= #₍(1) in the neighborhood of the

This implies that, ,⁾ꢄ_ꢈ(ꢆ¹, ꢈ^[2]) = 0. Moreover,

equilibrium point ꢆ¹(1, 0, 0).

0

∗

ˆ

©

-

0⁽

2

ª

®

2. When the parameter ꢈ attains the bifurcation threshold ꢈ = #₍(( ) +

ꢅꢄ_ꢈ(ꢆ¹, ꢈ^[2]) =

,

#_ꢃ(ꢃ^∗), the system (3) near the equilibrium ꢆ²((^∗, ꢃ^∗, 0) exhibits

−1

«

¬

(i). No saddle-node bifurcation occurs but

(ii). A transcritical bifurcation is observed.

0

−ꢆ # (1)2

0

1

©

-

ª ©

® -

ª

®

©

ª

®

DF_ꢈ(ꢆ¹, ꢈ^[2])+ =

=

.

-

−1

−2

«

¬ «

¬

«

¬

Alemu et.al (2026)

8

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

¤

Hence, ,⁾[ꢅꢄ_ꢈ(ꢆ¹, ꢈ^[2])+] = −23 ≠ 0. Furthermore,

3ꢊ

3ꢅ

ꢅ₂+ ꢊꢅ

¹. Thus, at ꢅ = ꢅ^∗it is reduced to

This implies that,

= −

ꢅ₁+ 2ꢊ

−2E₁²− 2ꢆ #⁰(1)E E

1

0

1

3

(

ꢈ

ꢉ

ꢈ

ꢉ

√

2

©

-

ª

®

ꢅ²ꢄ(ꢆ¹, ꢈ^[2])(+, +) =

¤

3ꢊ

3ꢅ

ꢅ₂+ ꢊꢅ

−ꢅ

ꢅ

1

2

2#⁰(1)E E

= −

=

+ 8

2ꢅ

.

√

ꢅ

2

ꢅ₁+ 2ꢊ

2ꢅ

1

3

«

¬

∗

2

(

ꢅ=ꢅ

ꢊ=8

which implies that,

ꢉ

√

¤

3ꢊ₉

ꢅ

1

ꢅ

2

Hence, '4

= −

≠ 0.

,⁾[ꢅ²ꢄ(ꢆ¹, ꢈ^[2])(+, +)] = 2#⁰(1)E E = −232²ꢆ ꢈ#⁰(1) ≠ 0.

3ꢅ

2ꢅ

∗

2

ꢅ=ꢅ

1

3

1

(

Hence, the transversality condition

Hence, based on the Sotomayor’s theorem as in Pirayesh et al. (2016) and

ꢂ

ꢃ

Yu et al. (2020) the model exhibit transcritical bifurcation at ꢈ = ꢈ^[2]

#₍(1) near to the axial equilibria ꢆ¹(1, 0, 0).

=

3ꢊ₉

<

≠ 0, 9 = 2, 3,

3ꢅ

∗

ꢅ=ꢅ

See the proof of (2) in the appendix A.

is satisﬁed, which conﬁrms the occurrence of a Hopf bifurcation at ꢅ = ꢅ^∗.

Furthermore, it can be demonstrated that there exists a threshold value of

the parameter ꢄ at which the present model also demonstrate a stability

switch via Hopf bifurcation.

T heorem 10 (Hopf Bifurcation). The system undergoes a Hopf bifurcation

near to the equilibrium point ꢆ²((^∗, ꢃ^∗, 0) at the parameter value ꢅ = ꢅ^∗,

provided the following criteria are hold:

1. At ꢅ = ꢅ^∗, we have ꢅ = 0 and ꢅ > 0, which guarantees the existence

1

2

T heorem 11. If the bifurcation parameter is given by

ꢆ ꢈꢄ ( − ꢆ #⁰(0)((1 − ()

of a pair of purely imaginary eigenvalues, and

2. The transversality condition is satisﬁed, i.e.,

ꢂ

ꢃ

3ꢊ₉

1

2

<

≠ 0, 9 = 2, 3,

ꢃ

3ꢅ

ꢇ^[0]

=

,

∗

ꢅ=ꢅ

ꢆ ꢈ

1

where ꢊ₉denote the eigenvalues of the auxiliary equation

then the model system (3) at the infection-free ﬁxed point ꢆ³((^∗, 0, %^∗) does not

exhibit a saddle-node bifurcation. Instead, the transcritical bifurcation of the

system is observed. See the proof in the appendix B.

ꢊ²+ ꢅ ꢊ + ꢅ = 0

1

2

associated with ꢆ². Here, ꢅ = −(ꢋ + ꢋ ) and ꢅ = ꢋ ꢋ − ꢌ ꢌ ,

2

1

2

1

2

with ꢋ₈and ꢌ₈(8 = 1, 2) deﬁ¹ned in th¹e proof part.

5 Computational Analysis

Proof. The auxiliary equation of the model (3) at ꢆ²is

In order to validate the theoretical results, the researchers numerically

explore the dynamic behavior of their model using the ode45 solver in

MATLAB. Owing to the unavailability of empirical data, a biologically

feasible and representative set of parameter values is adopted for the

ꢊ

∗

ꢊ

ꢋ − ꢊ

ꢌ

−ꢆ # (( )

1

(

ꢊ

∗

ꢌ

ꢋ − ꢊ −ꢆ # (ꢃ ) = 0,

2

ꢃ

ꢊ

purpose of numerical simulations. - = {ꢁ = 0.001, ꢁ = 0.007, ꢄ = 2.5,

1

2

ꢊ

ꢆ

1

= 0.04, ꢆ = 0.02, ꢇ = 0.04, ꢃ = 1.3, ꢈ = 0.01, ꢅ = 0.03, ℎ = 0.064,

2

0

ꢋ − ꢊ

3

and 4 = 0.05}. Moreover, for simulation purposes, we consider four

representative models selected from the sixteen possible combinations of

Holling-type functional responses. Speciﬁcally:

ꢄꢃ^∗

ꢄ(^∗

where, ꢋ = 1 − 2(^∗− ꢃꢃ^∗−

, ꢋ

2

=

− ꢇ, ꢋ = #₍((^∗) +

1

3

1 + ꢅꢃ^∗

(1 + ꢅꢃ^∗)²

ꢄ(^∗

ꢄꢃ^∗

# (ꢃ^∗) − ꢈ, ꢌ = −ꢃ(^∗

−

< 0 and ꢌ

=

> 0. After

1 + ꢅꢃ^∗

1

2

ꢃ

(1 + ꢅꢃ^∗)²

• Model 1: Represents the disease model with(HT-I–HT-II).

• Model 2: Corresponds to the combination (HT-II–HT-III).

• Model 3: Deﬁned by the combination (HT-III–HT-II).

• Model 4: Represents the combination (HT-IV–HT-III).

simpliﬁcation

(ꢋ₃− ꢊ)(ꢊ²+ ꢅ ꢊ + ꢅ ) = 0

(22)

1

2

where

ꢂ

ꢃ

ꢄꢃ^∗

ꢄ(^∗

ꢅ

= −(ꢋ₁+ ꢋ ) = − 1 − 2(^∗− ꢃꢃ^∗−

+

− ꢇ

,

1

2

1 + ꢅꢃ^∗

(1 + ꢅꢃ^∗)²

ꢂ

ꢃ ꢂ

ꢃ

ꢄꢃ^∗

ꢄ(^∗

= ꢋ ꢋ₂− ꢌ ꢌ

=

1 − 2(^∗− ꢃꢃ^∗−

− ꢇ +

1

2

1 + ꢅꢃ^∗

(1 + ꢅꢃ^∗)²

These four sample models are selected as representative cases among the

sixteen Holling-type response function combinations to capture diﬀerent

nonlinear interaction mechanisms and to compare their qualitative

impacts on the disease transmission dynamics. When the consumer

species is absent in the dynamics 3, then the model is dominated by prey

population. The time-series plots in Figure (1) show that both ((C) and

ꢃ(C) converge smoothly to their ﬁxed point values ꢆ²((^∗, ꢃ^∗, 0), indicating

local asymptotic stability. The infected prey population increases initially

due to infection transmission, then stabilizes at a moderate level, while

the susceptible population decreases and reaches a steady state.

ꢂ

ꢃ

ꢄ(^∗

ꢄꢃ^∗

ꢃ(^∗+

.

1 + ꢅꢃ^∗

(1 + ꢅꢃ^∗)²

Eq’n (22) has pure imaginary roots if ꢅ = 0 and ꢅ > 0 from which the

1

2

√

threshold value ꢅ = ꢅ^∗. Thus, when ꢅ = ꢅ^∗, ꢊ = ꢋ , ꢊ = 8 ꢅ , and

1

3

2

√

ꢊ

= −8 ꢅ . Diﬀerentiating equation (22) with respect to ꢅ, we have

3

2

3ꢊ

3ꢅ

3ꢊ 3ꢅ

1

2

3ꢅ

2ꢊ

+ ꢊ

+ ꢅ

+

= 0.

3ꢅ

¹3ꢅ

Alemu et.al (2026)

9

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

(a)

(b)

(c)

(d)

Figure 1: Time series plot of the system (3), where the parametric values ꢃ = 0.3; ꢄ = 0.8; ꢅ = 0.2; ꢆ = 0.4; ꢆ = 0.3; ꢇ = 0.25; G8 = 0.5; 0₁= 0.6; 0₂= 0.5; 1₁= 0.4;

1

2

1₂= 0.3; 2₁= 0.2; and initial condition (0.70, 0.120, 0.3).

When host is absent in the model system (3), the dynamics reduce to

a predator-prey subsystem involving ( and %. In Figure (2) the time

series plots show both populations converging to the host-free ﬁxed point

ꢆ³((^∗, 0, %^∗). Consumers grow up initially fueled by prey availability,

then stabilize as prey density decreases. Additionally, the phase diagrams

conﬁrm that trajectories approach the subspace ꢃ = 0. The infection-free

ﬁxed point is locally asymptotically stable under the parameter sets

considered, consistent with Theorem 6.

predation in prey dynamics and the stabilizing inﬂuence of functional

response saturation.

The bifurcation diagrams in Figures 4 conﬁrm the predicted transcritical

bifurcations in system (3) (Theorem 9). As the parameter ꢇ cross its

thresholds, equilibria exchange stability, with the infected equilibrium ꢃ^∗

smoothly transitioning from stable to unstable.

Figure 3 shows the system dynamics begin with periodic oscillation and

through time it goes to a locally asymptotically stable endemic ﬁxed point,

where the computational laboratory is performed for some possible

Holling Type response function combinations of the mathematical

Eco-Epidemiology model for the diseased-model (3).

Ecologically, small changes in disease-induced mortality ꢇ can shift

the system between disease-free and infected states, or from predator

extinction to coexistence, reﬂecting the inﬂuence of nonlinear functional

responses #₍and #_ꢃ.

Models 1 and 2 show bifurcations at lower thresholds, indicating higher

sensitivity and faster infection spread under simpler responses. In

contrast, Models 3 and 4, show stronger saturation or complex predation

terms, display delayed transitions, highlighting greater ecological

resilience.

Overall, the simulations show that both predator-free and disease-free

equilibria are stable, while nonlinear functional responses mainly aﬀect

the rate and amplitude of transient dynamics rather than the ﬁnal steady

state. These results emphasize the regulatory roles of infection and

Alemu et.al (2026)

10

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Figure 2: Time series plot of the system (3), where the parametric values ꢃ = 0.10; ꢄ = 0.4; ꢅ = 0.10; ꢇ = 0.30; ꢆ = 0.50; ꢆ = 0.70; ꢈ = 0.40; 0 = 0.65; 1 = 0.3;

1

2

ℎ = 0.20; ℎ₁= 0.30; ℎ₂= 0.25; and initial condition (0.80, 0.1, 0.3).

(a)

(b)

Figure 3: Time seires plot of the model system 3, where the parameter values ꢄ = 1.2, ꢆ = 0.5, ꢆ = 0.3, ꢇ = 0.3, ꢃ = 0.1, ꢈ = 0.4, 0 = 0.6, 1 = 0.3, 2 = 0.2, ꢅ = 0.03,

1

2

ℎ₁= 0.3, ℎ₂= 0.25 and for the, initial condition (0.8, 0.1, 0.3).

Alemu et.al (2026)

11

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Figure 4: Bifurcation diagrams for Models 1–4 with respect to ꢇ = ꢄ near the equilibrium ꢆ¹(1, 0, 0). The horizontal axis represents the bifurcation parameter ꢇ, and the

vertical axis represents the infected equilibrium ꢃ^∗. Solid lines denote stable equilibria, while dashed lines denote unstable equilibria and the bifurication value is ꢇ = 0.4.

5.1 Impact of the inhibition rate, ꢅ

such as crowding, limited contact, or behavioral avoidance. Ecologically,

ꢅ represents density-dependent inhibition in the infection process due

to immunity, crowding, or behavioural avoidance among prey. As ꢅ

increases, the eﬀective contact rate between susceptible and infected

prey decreases, reducing infection pressure. This reduction weakens

the oscillatory feedback between prey and predator populations, thereby

promoting stability in the coexistence equilibrium.

The inhibition rate ꢅ appears in the infection term which regulates the rate

at which susceptible prey become infected. From Figure(5), the parameter

ꢅ controls the saturation level of the infection process for small values of ꢅ,

the incidence rate is almost linear in ꢃ, leading to rapid spread of infection;

for large ꢅ, the infection saturates quickly, representing inhibitory eﬀects

Figure 5: Time series showing the impact of the inhibition rate ꢅ on system stability for HT-II–HT-III, where parameters value ꢅ = 0.1, 0.3, 0.5, ꢃ = 0.5; ꢄ = 0.8; ꢆ = 0.2;

1

ꢆ

2

= 0.15; ꢇ = 0.6; ꢈ = 0.5; 0 = 0.6; 1 = 0.3; 2 = 0.2; and initial condition (0.7, 0.2, 0.1).

5.2 Impact of the transmission rate, ꢄ

transmission rates ꢄ increasing the infected prey population rise up, and

making the nature of stability of the coexistence equilibrium is becomes

more periodic and take long time to stable. Moreover, for diﬀerent

Holling Type response functions combination, the patterns of stability of

endemic equilibrium point are identical.

As illustrated in Figure 6, the transmission parameter ꢄ has also have a

signiﬁcation eﬀect on the dynamical behaviour of population ꢃ , and as

Alemu et.al (2026)

12

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Figure 6: Time series plot of the dynamical system ( 3) for diﬀerent values of transmission rate ꢄ = 0.4, 0.6, 0.8 for Mode 1-4, where other parameter values ꢃ = 0.5,

ꢅ = 0.3, ꢆ = 0.2, ꢆ = 0.15, ꢇ = 0.6, ꢈ = 0.5; 0 = 0.6; 1 = 0.3; 2 = 0.2 and initial condition (0.7, 0.2, 0.1).

1

2

(a)

(b)

Figure 7: Time series plot(periodic solution) and phase diagram (limit cycle) of the model system HT-II–HT-III, where the parameter values ꢃ = 0.4, ꢄ = 5.0, ꢇ = 0.6,

= 1.2, ꢆ = 1.0, 0₁= 1.5, 0₂= 1.3, 1 = 0.4, ℎ = 0.2, ꢈ = 0.3, ꢅ = 2.0, and initial condition (0.6, 0.3, 0.1).

ꢆ

1

2

Figure 7 indicates existence of Hopf bifurcation which veriﬁes Theorem

4.6. Ecologically, ꢅ measures the strength of inhibitory (saturation)

eﬀects regulating predation or disease transmission. For ꢅ < ꢅ^∗, the

populations coexist at a stable steady state. When ꢅ > ꢅ^∗, the equilibrium

loses stability and sustained oscillations emerge due to feedback between

infection spread and predation pressure. Increased infection enhances

predator growth, which subsequently suppresses the host population,

leading to predator decline and eventual host recovery. This recurring

mechanism generates persistent population cycles, reﬂecting realistic

eco-epidemiological ﬂuctuations observed in natural ecosystems.

6 Result

In this section, we concisely summarize the main analytical and numerical

ﬁndings obtained in Sections 4 and 5 for the eco-epidemiological

model (3). The local and global stability conditions of all equilibrium

points, corresponding to diﬀerent combinations of Holling-Type

functional responses, are presented in Table (3). These results establish

the parametric regimes under which the system exhibits disease-free,

endemic, predator-free, or coexistence states.

The bifurcation analysis results of the eco-epidemiological dynamics 3

near equilibrium points are summarized in table (4).

Alemu et.al (2026)

13

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Table 3: Stability analysis result of Equilibrium points of the model system 3: Note; LAS ≡ locally asymptotically stable, GAS ≡ globally asymptotically stable

Equilibria

Stability conditions

Stability status

ꢆ⁰

ꢆ¹

ꢆ²

ꢆ³

ꢆ^∗

Always

Unstable

LAS

ꢄ < ꢇ and #₍(1) = ꢇ

Conditions stated in Theorem (5) (i)–(iii)

LAS

Conditions stated in Theorem (6) (1)–(3)

#₍(() ≥ exp

LAS

∫

3D

(

)

0

GAS

Conditions stated in Theorem (8)(i)–(iii)

LAS

The collective results depicted in Figures 4 clearly demonstrate how

variations in the key bifurcation parameter ꢇ regulate the coexistence

and persistence of prey, infected prey, and predator populations. The

transcritical bifurcation marks a critical threshold where the system shifts

from a disease-free to an endemic equilibrium, reﬂecting a change in

ecological stability and disease prevalence.

From an ecological perspective, increasing the recovery rate (ꢈ) helps

drive the system toward a disease-free state. Comparing Models 1–4,

introducing nonlinear saturation in infection and predation stabilizes

the system by postponing bifurcations. This highlights the importance

of using realistic functional responses in eco-epidemiological models

to capture key biological feedbacks and better understand ecosystem

resilience under disease pressure.

Table 4: Local bifurcation analysis result of Equilibrium points of the model 3: TB ≡ Transcritical bifurcation, HB ≡ Hopf bifurcation

E.P

ꢆ¹

ꢆ²

T hreshold value

ꢄ = ꢇ

Stability condition

Bifurcation

ꢅE ≠ E

2

TB

HB

1

#₍(1) = ꢈ

always

ꢈ = #₍(( ) + #_ꢃ(ꢃ )

#⁰((^∗)+ + #⁰(ꢃ^∗)+ ≠ 0

∗

ˆ

1

2

(

ꢃ

ꢈ

ꢉ

3ꢊ₉

ꢅ = ꢅ^∗

ꢅ

1

= 0, ꢅ > 0 and '4

≠ 0

2

3ꢅ

∗

ꢅ=ꢅ

0

ꢆ

1

ꢈꢄ (−ꢆ # (0)((1−()

ꢆ³

ꢇ^[0]

=

−(2ꢄꢅ(^∗+ ꢆ #⁰⁰(0)%^∗)+²Υ + 2(ꢄ+ + Υ − ꢆ #⁰(0)+ + )Υ ≠ 0

TB

2

ꢃ

2

1

2

3

ꢆ

1

ꢈ

2

ꢃ

The inhibition (saturation) parameter ꢅ plays a critical regulatory

role. When ꢅ is below the critical threshold (ꢅ^∗

2), the system

results.

=

settles into a stable coexistence of susceptible prey, infected prey, and

predators. However, once ꢅ exceeds this value, a Hopf bifurcation occurs:

the equilibrium becomes unstable and a stable limit cycle emerges.

Biologically, this leads to recurring oscillations driven by feedback

between disease transmission and predation. Increased susceptible

prey boosts infection and predator growth; predators then reduce prey

populations, which in turn lowers predator numbers, allowing prey to

recover and restarting the cycle.

The system exhibits oscillatory behavior for lower values of the inhibition

rate (ꢅ), whereas higher inhibition rates promote stability. Hopf

bifurcation analysis, taking ꢅ as the bifurcation parameter, revealed

that increasing inhibition enhances system stability. Furthermore,

when the predation rates ($ , $ ) exceed a critical threshold, the

predator-free equilibrium becomes²unstable, and a stable disease-free

coexistence of prey and predator emerges. The bifurcation analysis

indicates that disease transmission and predator–prey interactions jointly

determine ecosystem stability. Managing infection parameters such

as the transmission rate can prevent oscillatory outbreaks and species

extinction. Hence, controlling ecological feedbacks through parameter

tuning plays a crucial role in maintaining biodiversity and long-term

coexistence within predator–prey systems. Overall, the theoretical and

numerical investigations are carried out under saturating incidence rates

demonstrate the biological consistency of the proposed model. The

results provide valuable insights into the interplay between infection,

predation eﬃciency, and population stability in eco-epidemiological

systems. The primary contribution of this work lies in providing a

comprehensive bifurcation analysis under these combined nonlinear

mechanisms. We rigorously establish threshold dynamics through

the basic reproduction number and employ bifurcation theory to

demonstrate the occurrence of transcritical and Hopf bifurcations.

The results reveal how inhibition and transmission parameters govern

transitions between disease-free equilibria, endemic coexistence, and

sustained oscillatory outbreaks.

1

Overall, the Hopf bifurcation shows how changes in inhibitory eﬀects

can shift the ecosystem from stable coexistence to sustained oscillations,

underscoring the delicate balance between disease dynamics and

predator–prey interactions.

7 Conclusion

In this study, the researchers have investigated an eco-epidemiological

mathematical model in which

a

prey population is infected by

microparasites, while predators feed on both susceptible and infected

prey following a general Holling-type functional response. The model

was developed to explore how disease transmission and predation

eﬃciency aﬀect the overall community structure and population

dynamics. An emergent carrying capacity was introduced to reﬂect the

fact that infected prey, having reduced ﬁtness, are more easily captured

by predators. The stability and bifurcation conditions were derived

for diﬀerent equilibrium points, including trivial, axial, predator-free,

disease-free, and endemic states. Analytical and numerical analyses

showed strong agreement between theoretical predictions and simulation

Future studies can extend this work by incorporating time delays

representing disease incubation or predator gestation periods, which

may lead to more complex dynamical behaviors such as chaos or

multiple attractors. Furthermore, integrating stochastic eﬀects, seasonal

Alemu et.al (2026)

14

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

variations, or optimal control strategies may enhance the model

applicability to real-world ecological management and conservation

policies.

Haldar, S., Khatua, A., Das, K., & Kar, T. K. (2021). Modeling and analysis

of a predator–prey type eco-epidemic system with time delay.

Modeling Earth Systems and Environment, 7, 1753–1768.

Hale, J. K. (2009). Ordinary diﬀerential equations. Courier Corporation.

Haque, M., & Venturino, E. (2007). An ecoepidemiological model with

disease in predator: The ratio‐dependent case. Mathematical

methods in the Applied Sciences, 30(14), 1791–1809.

Data Availability Statement

Holling, C. S. (1959). The components of predation as revealed by a study

of small-mammal predation of the european pine sawﬂy1. The

canadian entomologist, 91(5), 293–320.

The data supporting the ﬁndings of this study are available from the

authors upon reasonable request.

Hu, Z., Teng, Z., Zhang, T., Zhou, Q.,

&

Chen, X. (2017).

Globally asymptotically stable analysis in a discrete time

eco-epidemiological system. Chaos, Solitons & Fractals, 99,

20–31.

Conflicts of interest

Hugo, A., & Simanjilo, E. (2019). Analysis of an eco-epidemiological

model under optimal control measures for infected prey.

Applications and Applied Mathematics: An International Journal

(AAM), 14(1), 8.

The authors declare that they have no conﬂicts of interest relevant to this

study.

Kooi, B. W., van Voorn, G. A., & pada Das, K. (2011). Stabilization and

complex dynamics in a predator–prey model with predator

suﬀering from an infectious disease. Ecological Complexity, 8(1),

113–122.

Layek, G. C. (2015). An introduction to dynamical systems and chaos

(Vol. 449). Springer.

Author Contributions

All have equal contribution.

Li, H., & Takeuchi, Y. (2011). Dynamics of the density dependent

predator–prey system with beddington–deangelis functional

response. Journal of Mathematical Analysis and Applications,

374(2), 644–654.

Liu, W. M., Hethcote, H. W., & Levin, S. A. (1987). Dynamical behavior of

epidemiological models with nonlinear incidence rates. Journal

of mathematical biology, 25, 359–380.

Funding

This research received no speciﬁc grant from any funding agency.

Lotka, A. J. (1925). Elements of physical biology. Williams; Wilkins.

Maiti, A. P., Jana, C., & Maiti, D. K. (2019). A delayed eco-epidemiological

model with nonlinear incidence rate and crowley–martin

functional response for infected prey and predator. Nonlinear

Dynamics, 98, 1137–1167.

Omar, F. M., Sohaly, M. A., & El-Metwally, H. (2024). Lyapunov

functions and global stability analysis for epidemic model with

n-infectious. Indian Journal of Physics, 98(5), 1913–1922.

Panja, P. (2020). Prey–predator–scavenger model with monod–haldane

type functional response. Rendiconti del Circolo Matematico di

Palermo Series 2, 69(3), 1205–1219.

Pirayesh, B., Pazirandeh, A., & Akbari, M. (2016). Local bifurcation

analysis in nuclear reactor dynamics by sotomayor’s theorem.

Annals of Nuclear Energy, 94, 716–731.

Ruan, S., & Wang, W. (2003). Dynamical behavior of an epidemic model

with a nonlinear incidence rate. Journal of diﬀerential equations,

188(1), 135–163.

References

Allen, L. J., Bolker, B. M., Lou, Y., & Nevai, A. L. (2007). Asymptotic

proﬁles of the steady states for an sis epidemic patch model.

SIAM Journal on Applied Mathematics, 67(5), 1283–1309.

Anderson, R. M., & May, R. M. (1986). The invasion, persistence

and spread of infectious diseases within animal and plant

communities. Philosophical Transactions of the Royal Society of

London. B, Biological Sciences, 314(1167), 533–570.

Bezabih, A. F., Edessa, G. K., & Rao, K. P. (2021). Mathematical modeling

the impact of infectious diseases in prey-predator interactions.

Biswas, S., Samanta, S., & Chattopadhyay, J. (2015). A model based

theoretical study on cannibalistic prey–predator system

with disease in both populations. Diﬀerential Equations and

Dynamical Systems, 23, 327–370.

Brauer, F. (2005). The kermack–mckendrick epidemic model revisited.

Mathematical biosciences, 198(2), 119–131.

Saifuddin, M., Biswas, S., Samanta, S., Sarkar, S., & Chattopadhyay, J.

(2016). Complex dynamics of an eco-epidemiological model

with diﬀerent competition coeﬃcients and weak allee in the

predator. Chaos, Solitons & Fractals, 91, 270–285.

Cai, L. M., & Li, X. Z. (2010). Global analysis of a vector-host epidemic

model with nonlinear incidences. Applied Mathematics and

Computation, 217(7), 3531–3541.

Sieber, M., Malchow, H., & Hilker, F. M. (2014). Disease-induced

modiﬁcation of prey competition in eco-epidemiological

models. Ecological complexity, 18, 74–82.

Capasso, V.,

kermack-mckendrick

Mathematical biosciences, 42(1-2), 43–61.

&

Serio, G. (1978).

A

generalization of the

deterministic

epidemic model.

Volterra, V. (1927). Fluctuations in the abundance of a species considered

mathematically. Nature, 119(2983), 12–13.

Das, K. P. (2016). A study of harvesting in a predator–prey model with

disease in both populations. Mathematical Methods in the Applied

Sciences, 39(11), 2853–2870.

Dawed, M. Y., Tchepmo Djomegni, P. M., & Krogstad, H. E. (2020).

Complex dynamics in a tritrophic food chain model with

general functional response. Natural Resource Modeling, 33(2),

e12260.

Demir, M. (2019). Optimal control strategies in ecosystem-based ﬁshery

models.

Wang, N., Zhao, M., Yu, H., Dai, C., Wang, B., & Wang, P. (2016).

Bifurcation behavior analysis in

a

predator‐prey model.

Discrete Dynamics in Nature and Society, 2016(1), 3565316.

Wayesa, N. N., Obsu, L. L., & Dawed, M. Y. (2024). Analysis of

predator–prey model with inclusion of temperature variability

in prey refugees. Journal of Applied Mathematics, 2024(1),

5138320.

Wayesa, N. N., Obsu, L. L., & Dawed, M. Y. (2025). Predator–prey

population dynamics with time delay and prey refuge eﬀects.

Modeling Earth Systems and Environment, 11(2), 142.

Ghanbari, B. (2021). On the modeling of an eco-epidemiological model

using a new fractional operator. Results in physics, 21, 103799.

Gumel, A. B., & Moghadas, S. M. (2003). A qualitative study of

Yu, X., Zhu, Z., Lai, L., & Chen, F. (2020). Stability and bifurcation

a

vaccination model with non-linear incidence. Applied

analysis in

a

single-species stage structure system with

mathematics and computation, 143(2-3), 409–419.

michaelis–menten-type harvesting. Advances in Diﬀerence

Equations, 2020, 1–18.

Alemu et.al (2026)

15

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

Appendix A

2

∗

2

ˆ

let the Jacobian matrix the model(3) at the predator free equilibrium point ꢆ (( , ꢃ , 0) denote by ꢇ(ꢆ ) = (ꢀ₈₉)_3×3

ꢄꢃ^∗

1 + ꢅꢃ^∗

ꢄ(^∗

(1 + ꢅꢃ^∗)²

1 − 2(^∗− ꢃꢃ^∗−

−ꢃ(^∗−

−ꢆ # ((^∗)

1

(

©

ª

-

®

-

®

ꢄꢃ^∗

1 + ꢅꢃ^∗

ꢄ(^∗

2

ꢇ(ꢆ ) =

.

− ꢇ

−ꢆ # (ꢃ^∗)

2

ꢃ

(1 + ꢅꢃ^∗)²

-

®

-

®

0

#₍((^∗) + #_ꢃ(ꢃ^∗) − ꢈ

«

¬

From the condition at which ꢇ(ꢆ²) has zero eigenvalue, that is, ꢊ = ꢋ = # ((^∗) + # (ꢃ^∗) − ꢈ = 0 the bifurcation value is

1

3

(

ꢃ

∗

ˆ

ꢈ = #₍(( ) + #_ꢃ(ꢃ ).

2

)

ˆ

2

ˆ

Now we compute the Jacobian matrix ꢇ(ꢆ ) = (ꢀ )_3×3at ꢈ = ꢈ which is same as above ꢇ(ꢆ ) except ꢀ = 0. The eigenvectors of ꢇ(ꢆ , ꢈ) and ꢇ (ꢆ , ꢈ),

33

89

corresponding to the zero eigenvalue are, respectively

Ψ e

+

0

1

e

1

©

-

ª

®

©

-

ª

®

© ª

+

0

+ =

=

and , =

- ®

2

Ψ e

+

3

f

2

«

¬

«

¬

« ¬

ˆ

21

ˆ

21

ˆ

ꢀ

ꢀ₂₃− ꢀ

ꢀ

ˆ

ꢀ

ˆ

ꢀ

ꢀ₂₂− ꢀ

ꢀ

12

ˆ

22

ˆ

11

12

ˆ

where, Ψ

=

¹³, Ψ

=

²¹, moreover e and f are nonzero real numbers. From our model system (3), we have:

1

2

ˆ

ꢀ₁₃− ꢀ

11 23

0

©

-

ª

®

© ª

2

ˆ

0

ꢄ (-, ꢈ) =

ꢈ

=⇒ ꢄ (ꢆ , ꢈ) =

- ®

ꢈ

−%

0

«

¬

« ¬

)

2

ˆ

Thus, , ꢄ_ꢈ(ꢆ , ꢈ) = 0. Applying Sotomayor’s theorem (Pirayesh et al., 2016) for local bifurcation, the saddle node bifurcation does not occur near to

the equilibrium point ꢆ²((^∗, ꢃ^∗, 0). For Bogdanov– Takens bifurcation, there must be two equilibria : saddle and non-saddle. Therefore, BT bifurcation

cannot appear here also.

We noted that the ﬁrst condition , ꢄ_ꢈ(ꢆ , ꢈ) = 0 of Sotomayor’s theorem for the existence of transcritical bifurcation is satisﬁed. Now,

)

2

ˆ

0

Ψ e

0

1

e

©

-

ª

®

©

-

ª ©

® -

ª

®

©

-

ª

®

2

ˆ

ꢅꢄ (ꢆ , ꢈ) =

=⇒ ꢅꢄ (ꢆ , ꢈ)+ =

=

ꢈ

−1

−1 Ψ e

−Ψ e

2

«

¬

«

¬ «

¬

«

¬

)

2

ˆ

So, we have , [ꢅꢄ (ꢆ , ꢈ)+] = −Ψ ef ≠ 0. Moreover,

2

ꢈ

ꢂ

ꢃ

ꢀ

ꢁ

2ꢄꢅ(^∗

ꢄ

−2+²

+

− 2 ꢃ +

+ + + ꢆ #⁰((^∗)+ +

2

©

ª

®

1

2

1

3

2

1

2

(

(1 + ꢅꢃ^∗)³

(1 + ꢅꢃ^∗)²

-

®

ꢂ

ꢃ

−2ꢄꢅ(^∗

ꢄ

2

ˆ

+ + − ꢆ #⁰(ꢃ^∗)+ +

ꢅ ꢄ(ꢆ , ꢈ)(+, +) =

2

+

+ 2

1

2

3

2

ꢃ

-

(1 + ꢅꢃ^∗)³

(1 + ꢅꢃ^∗)²

®

-

®

2+ (#⁰((^∗)+ + #⁰(ꢃ^∗)+ )

3

1

2

«

¬

(

ꢃ

0

(

0

ꢃ

)

2

∗

ˆ

Thus, we have , [ꢅ ꢄ(ꢆ , ꢈ)(+, +)] = 2+ f(# (( )+ + # (ꢃ )+ ) ≠ 0. Therefore, by Sotomayor’s theorem, transcritical bifurcation occurs near to the

3

1

2

predator-free stationary point ꢆ²((^∗, ꢃ^∗, 0).

Appendix B

The Jacobian matrix of the system (3) at the infection free equilibrium point ꢆ³((^∗, 0, %^∗) denote by ꢇ(ꢆ³) = (ꢅ₈₉)_3×3as

1 − 2(^∗− ꢆ #⁰((^∗)%^∗

−ꢃ(^∗− ꢄ(^∗

−ꢆ # ((^∗)

1

(

©

«

ª

-

®

ꢇ(ꢆ³) =

0

ꢄ(^∗− ꢇ − ꢆ #⁰(0)%^∗

0

.

2

ꢃ

-

®

#₍⁰((^∗)%^∗

#⁰_ꢃ(0)%^∗

#₍((^∗) − ꢈ

¬

Thus, ꢇ(ꢆ³) has zero eigenvalue, while, ꢊ = 5₂= ꢄ (^∗− ꢇ − ꢆ #⁰(0)%^∗= 0. and the model bifurcate when

1

2

ꢃ

ꢆ ꢈꢄ(^∗− ꢆ #⁰(0)(^∗(1 − (^∗)

1

2

ꢃ

ꢇ^[0]

=

.

ꢆ ꢈ

1

Alemu et.al (2026)

16

East Afr. J. Biophys. Comput. Sci. (2026), Vol. 7, Issue. 1, 1-17

To perform the Jacobian matrix ꢇ^[0](ꢆ³) = (ꢅ )_3×3at ꢇ = ꢇ^[0]which is same as above ꢇ(ꢆ³) except ꢅ = 0. The eigenvectors of ꢇ^[0](ꢆ³, ꢇ^[0]) and

22

89

(ꢇ^[0](ꢆ³, ꢇ^[0])) , corresponding to the zero eigenvalue are, respectively

)

Υ

©

ª

-

®

ꢅ

ꢅ₃₁− ꢅ

ꢅ

32 13

+

0

13

12

11 33

1

©

-

ª

®

© ª

Υ

+

Γ

+ =

=

and , =

- ®

2

ꢅ₃₃− ꢅ

-

®

+

3

0

«

¬

« ¬

ꢅ

ꢅ₃₂− ꢅ

ꢅ

11

12 13

Υ

ꢅ₁₂ꢅ₃₃− ꢅ

ꢅ

32 13

)

«

¬

where,Υ and Γ are nonzero real numbers. From our model system (3), use derivative we get:

0

−ꢃ

0

©

-

ª

®

© ª

=⇒ ꢄ (ꢆ³, ꢇ^[0]) =

0

ꢄ (-, ꢇ) =

- ®

ꢇ

0

«

¬

« ¬

Thus, applying Sotomayor’s theorem ﬁrst condition ,⁾ꢄ_ꢇ(ꢆ³, ꢇ^[0]

)

=

0. Hence, the dynamical system (3) saddle node bifurcation does not

demonstrate at the disease free equilibrium point ꢆ³((^∗, 0, %^∗).

Now we try to perform the other conditions of Sotomayor’s theorem for the existence of transcritical bifurcation. Thus,

0

−1

0

©

-

ª

®

ꢅꢄ (-, ꢇ) =

ꢇ

«

¬

Υ

0

©

ª

®

0

−1

0

ꢅ

ꢅ₃₁− ꢅ

ꢅ

13

11 33

-

©

-

ª

®

ꢅ

ꢅ₃₃− ꢅ

ꢅ

32 13

©

-

ª

®

Υ

11

12

13 31

=⇒ ꢅꢄ_ꢇ(ꢆ³, ꢇ^[0])+ =

=

Υ

ꢅ₃₃− ꢅ

- ꢅ¹²ꢅ − ꢅ³²ꢅ¹³

ꢅ₃₃− ꢅ

11 32

12 13

«

¬

Υ

0

«

¬

ꢅ

− ꢅ

ꢅ

)

«

12 33

32 13

¬

ꢅ

ꢅ₃₃− ꢅ

ꢅ

11

12

13 31

Hence, we arrived that ,⁾[ꢅꢄ_ꢇ(ꢆ³, ꢇ^[0])+] =

ΥΓ ≠ 0. In addition,

ꢅ₃₃− ꢅ

32 13

ꢀ

ꢁ

(−2 − ꢆ #⁰⁰((^∗)%^∗)+²+ 2ꢄꢅ(^∗+₂²− 2 (ꢃ + ꢄ1)+ + + ꢆ #⁰((^∗)+ +

1

2

1

3

2

1

©

(

ª

-

®

-

®

ꢅ²ꢄ(ꢆ³, ꢇ^[0])(+, +) =

(−2ꢄꢅ(^∗− ꢆ #⁰⁰(0)%^∗)+²+ 2(ꢄ+ + − ꢆ #⁰(0)+ + )

2

1

2

3

2

ꢃ

-

®

#⁰⁰((^∗)%^∗+₁²+ #⁰⁰(0)%^∗+₂²+ 2+ (#⁰((^∗)+ + #⁰(0)+ )

3

1

2

«

¬

(

ꢃ

(

ꢃ

Thus, we have

,⁾[ꢅ²ꢄ(ꢆ³, ꢇ^[0])(+, +)] = −(2ꢄꢅ(^∗+ ꢆ #⁰⁰(0)%^∗)+²Υ + 2(ꢄ+ + Υ − ꢆ #⁰(0)+ + )Υ ≠ 0.

2

1

2

3

2

ꢃ

Therefore, according to Sotomayor’s theorem (Pirayesh et al., 2016), our dynamical system (3) experience transcritical bifurcation at the disease-free

ﬁxed point ꢆ³((^∗, 0, %^∗).

Alemu et.al (2026)

17